A comprehensive 1,000 Genomes-based genome-wide association meta-analysis of coronary artery disease.

Existing knowledge of genetic variants affecting risk of coronary artery disease (CAD) is largely based on genome-wide association study (GWAS) analysis of common SNPs. Leveraging phased haplotypes from the 1000 Genomes Project, we report a GWAS meta-analysis of ∼185,000 CAD cases and controls, interrogating 6.7 million common (minor allele frequency (MAF) > 0.05) and 2.7 million low-frequency (0.005 < MAF < 0.05) variants. In addition to confirming most known CAD-associated loci, we identified ten new loci (eight additive and two recessive) that contain candidate causal genes newly implicating biological processes in vessel walls. We observed intralocus allelic heterogeneity but little evidence of low-frequency variants with larger effects and no evidence of synthetic association. Our analysis provides a comprehensive survey of the fine genetic architecture of CAD, showing that genetic susceptibility to this common disease is largely determined by common SNPs of small effect size.

Coronary artery disease (CAD) is the main cause of death and disability worldwide and represents an archetypal common complex disease with both genetic and environmental determinants[1,2]. To date, 48 genomic loci have been found to harbour common SNPs in genome-wide significant association with the disease. Previous GWAS of CAD have tested the common disease/common variant hypothesis with meta-analyses typically based on HapMap imputation training sets or tagging SNP arrays with up to 2.5 million SNPs (85% with MAF > 0.05)[3,4]. The 1000 Genomes Project[5] has considerably expanded the coverage of human genetic variation especially for lower frequency and insertion/deletion variants (indels). We assembled 60,801 cases and 123,504 controls from 48 studies for a GWAS meta-analysis of CAD; 34,997 (57.5%) of the cases and 49,512 (40.1%) of the controls had been previously included in our Metabochip-based CAD meta-analysis (Supplementary Fig. 1)[3]. Imputation was based on the 1000 Genomes phase 1 version 3 training set with 38 million variants of which over half are low frequency (MAF < 0.005) and one-fifth are common (MAF > 0.05) variants. The majority (77%) of the participants were of European ancestry; 13% and 6% were of south (India and Pakistan) and east (China and Korea) Asian ancestry with smaller samples of Hispanic and African Americans (Supplementary Table 1). Case status was defined by an inclusive CAD diagnosis (e.g. myocardial infarction (MI), acute coronary syndrome, chronic stable angina, or coronary stenosis >50%). After selecting variants that surpassed allele frequency (MAF > 0.005) and imputation quality control criteria in at least 29 (>60%) of the studies, 8.6 million SNPs and 836K (9%) indels were included in the meta-analysis (Fig. 1); of these, 2.7 million (29%) were low frequency variants (0.005 < MAF < 0.05).

Figure 1

Spectrum of minor allele frequencies (MAF) and median imputation quality (MEDIAN INFO) showing the number (N) of variants in each bin (a) shows the distribution for the 9.4M 1000 Genomes phase1v3 variants (b) shows the distribution for 2.5M HapMap2 SNPs. Imputation quality was calculated as the median of the respective values in up to 48 contributing studies; the imputation quality for genotyped variants was set equal to 1.0. The 1000 Genomes training set included more low frequency variants, many of which have imputation qualities > 0.9.

Scanning for additive associations

The results of an additive genetic model meta-analysis are summarized in Manhattan plots (Fig. 2 and Supplementary Fig. 2). 2,213 variants (7.6% indels) showed significant associations (P < 5 × 10−8) with CAD with a low false discovery rate (FDR q-value < 2.1 × 10−4). When these 2,213 variants are grouped into loci, eight represent regions not previously reported at genome-wide levels of significance (Fig. 2; Table 1). Of 48 previously reported loci at genome-wide levels of significance, 47 showed nominally significant associations (Supplementary Table 2). The exception was rs6903956, the lead SNP for the ADTRP–C6orf105 locus detected in Han Chinese[6], which previously showed no association in the Metabochip meta-analysis of Europeans and South Asians[3]. Thirty-six previously reported loci showed genome-wide significance (Supplementary Table 2). Monte Carlo simulations, guided by published effect-sizes, suggest that our study was powered to detect 34 of the previously reported loci (95%CI 31 – 41 loci) at genome-wide significance. Hence, our findings are fully consistent with the previously identified CAD loci. The majority of the loci showing GWAS significance in the present analysis were well imputed (82% with imputation quality > 0.9) (Fig. 3a) and had small effect sizes (odds ratio (OR) < 1.25) (Fig. 3b). An exception was the lead SNP in the novel chromosome 7q36.1 (NOS3) locus, rs3918226, which was only moderately well imputed (quality 0.78) but the validity of this association is supported by existing genotype data as it was present on the HumanCVD BeadChip available for some of the cohorts used in the present analysis and therefore directly measured genotypes could be compared with imputed genotypes (Supplementary Table 3)[7]. Three additional lower frequency and moderately well imputed SNPs in LPA and APOE (Fig. 3a), which were not previously reported in CAD GWAS[3,4], also showed strong associations (LPA: rs10455872 P = 5.7 × 10−39, rs3798220 P = 4.7 × 10−9; APOE: rs7412 P = 8.2 × 10−11). The LPA SNPs have been previously shown to be strongly associated with CAD in candidate gene studies based on experimental genotype data[7,8]. SNP rs7412 encodes the epsilon 2 allele of APOE and it has been well documented that ε2 carriers have lower cholesterol levels and significant protection from CAD was confirmed in a large meta-analysis[9] and in the Metabochip study (P = 0.0009)[3]. However, rs7412 is not present on most commercially available genome-wide genotyping arrays and cannot be imputed using HapMap reference panels, supporting the value of the expanded coverage of the 1000 Genomes reference panels. Finally, SNP rs11591147, which encodes the low frequency (MAF = 0.01) R46L variant in PCSK9 that has been associated with low LDL cholesterol levels and cardioprotection[10-13], was imperfectly imputed (quality = 0.61). Nonetheless these data provide the strongest evidence to date for a protective effect of this variant for CAD (P = 7.5 × 10−6).

Figure 2

A circular Manhattan plot summarizing the 1000 Genomes CAD association results. The meta-analysis statistics have been adjusted for over-dispersion (before double genomic control, lambda = 1.18); over-dispersion is predicted to be a regular feature in GWAS under the polygenic inheritance model[60]. The association statistics have been capped to P = 1 × 10−20. Genome-wide significant variants (P < 5 × 10−8) are indicated by red triangles. Novel CAD loci are named in red (Table 1). Previously reported loci showing genome-wide significance are shown in black and those showing nominal significance (P < 0.05) in our meta-analysis in blue (Supplementary Table 2). The inner track (see inset) shows the imputation quality score of the lead variants of the novel loci. The middle track shows numbered chromosome ideograms with the centromeres indicated by pink bars.

Table 1

Ten new CAD loci

									Association model
									Additive		Recessive
Lead variant	Locus name	Chr.	A1/A2	Effect allele (A1) freq.	Imputation quality	I2	Heterogeneity P	n studies	OR (95% CI)	P	OR (95% CI)	P
rs17087335	REST-NOA1	4	T/G	0.21	0.99	0.20	0.11	48	1.06 (1.04-1.09)	4.60E-08	1.11 (1.05-1.17)	3.30E-04
rs3918226	NOS3	7	T/C	0.06	0.78	0.15	0.19	45	1.14 (1.09-1.19)	1.70E-09	1.26 (0.99-1.60)	5.96E-02
rs10840293	SWAP70	11	A/G	0.55	0.94	0.17	0.16	47	1.06 (1.04-1.08)	1.30E-08	1.05 (1.02-1.09)	1.51E-03
rs56062135	SMAD3	15	C/T	0.79	0.98	0.00	0.67	46	1.07 (1.05-1.10)	4.50E-09	1.17 (1.10-1.25)	8.88E-07
rs8042271	MFGE8-ABHD2	15	G/A	0.9	0.93	0.16	0.19	46	1.10 (1.06-1.14)	3.70E-08	1.25 (1.13-1.37)	7.27E-06
rs7212798	BCAS3	17	C/T	0.15	0.95	0.14	0.21	48	1.08 (1.05-1.11)	1.90E-08	1.17 (1.07-1.28)	6.12E-04
rs663129	PMAIP1-MC4R	18	A/G	0.26	1.00	0.00	0.6	47	1.06 (1.04-1.08)	3.20E-08	1.11 (1.06-1.17)	7.15E-06
rs180803	POM121L9P-ADORA2A	22	G/T	0.97	0.86	0.00	0.67	41	1.20 (1.13-1.27)	1.60E-10	N/A	N/A
rs11830157	KSR2	12	G/T	0.36	0.99	0.14	0.22	42	1.04 (1.02-1.06)	3.90E-04	1.12 (1.08-1.16)	2.12E-09
rs12976411	ZNF507-LOC400684	19	T/A	0.09	0.93	0.50	5.09E-04	34	0.95 (0.92-0.99)	9.10E-03	0.67 (0.60-0.74)	1.18E-14

Association results are presented for two inheritance models; results in bold text indicate the discovery association model. P values have been adjusted for over-dispersion following meta-analysis. Heterogeneity P values are for the respective discovery association model. Chr., chromosome; A1, effect allele; A2, non-effect allele; freq., frequency; I, heterogeneity inconsistency index; OR, odds ratio; CI, confidence interval; N/A, not available due to insufficient numbers (<60%) of studies for reliable results. The n studies column shows the number of studies participated in the discovery result where up to 48 studies participated in the additive model and up to 43 studies in the recessive model meta-analysis.

Figure 3

Imputation quality and effect size of lead variants at 46 genome-wide significant loci. (a) Imputation quality and minor allele frequency (MAF) for lead variants at 46 genome-wide significant susceptibility loci. Blue circles indicate novel additive loci, red squares - novel recessive loci, black triangles - previously mapped additive loci, black diamonds - key SNPs in LPA and APOE. Imputation quality and MAF were calculated as the median of the respective values in up to 48 contributing studies; the imputation quality for studies with genotype data was fixed at 1.0. (b) Odds ratios and effect allele frequency (EAF) for lead variants at 46 genome-wide significant loci. Blue circles indicate novel additive loci; red squares - novel recessive loci, black triangles - previously mapped additive loci. SNPs rs55730499 and rs2891168 are lead variants in the LPA and chromosome 9p21 susceptibility loci. EAF was calculated as the median value in up to 48 contributing studies.

Scanning for non-additive associations

Few GWAS of CAD have systematically scanned for associations that include dominance effects and few truly recessive loci have been reported[14,15]. We used a recessive inheritance model to search for susceptibility effects conferred by homozygotes for the minor (i.e. less frequent) allele. Two novel recessive loci were identified with MAF 0.09 and 0.36 and genotypic OR 0.67 and 1.12 (Table 1, Fig. 2); these loci showed very little evidence of association under an additive model (Table 1). A supplementary dominant model analysis revealed multiple strong associations with variants that all overlapped with loci identified in the additive model analysis (Supplementary Table 4).

MI sub-phenotype analysis

Sub-group analysis in cases with a reported history of MI (~ 70% of total number of cases) did not identify any additional associations reaching genome-wide significance. The association results for the MI sub-phenotype for the 48 previously known CAD loci and the 8 novel additive CAD loci discovered in this study are shown in Supplementary Table 5. Supplementary Fig. 3 compares the OR for the lead SNPs at 56 loci for the broader CAD phenotype (full cohort) and the MI sub-phenotype. While, as expected, for most of the loci the ORs are very similar, for ABO and HDAC9 the ORs are sufficiently distinct for their 95% confidence intervals to lie away from the line of equality, suggesting that ABO preferentially associates with MI and HDAC9 with stable coronary disease but not infarction per se.

FDR and heritability analysis

We performed a joint association analysis to search for evidence of synthetic associations[16] where multiple low-frequency susceptibility variants at a locus might be in LD with a common variant that is discovered as the lead variant in a GWAS, and to compile an FDR-defined list of informative variants for annotation and heritability analysis[3]. Variants that showed suggestive additive associations (P < 5 × 10−5) were assigned to 214 putative susceptibility loci centred on a lead variant ±1 centiMorgan and all variants in these loci examined; consequently the search space for the joint analysis included 1,399,533 variants. Using the GCTA software[17] to perform an approximate joint association analysis (Online Methods), we identified 202 variants (q-value < 0.05) in 129 loci (Supplementary Table 6) with multiple (2 – 14) tightly linked variants in 57% of CAD loci. The 202 FDR variants were mostly common (median MAF 0.22) and well imputed (median imputation quality 0.97). 95 variants (explaining 13.3±0.4% of CAD heritability) mapped to 44 GWAS significant loci. 93 variants (explaining 12.9±0.4% of CAD heritability) mapped to loci that include a previously reported GWAS significant variant; 109 variants (explaining a further 9.3±0.3% of CAD heritability) mapped to other loci. Fifteen low frequency (MAF < 0.05) variants explained only 2.1±0.2% of CAD heritability noting that our study was ~90% powered to detect OR > 1.5 with low frequency variants (Supplementary Table 7).

Common variants showing typical GWAS signals might be coupled with one or more low frequency variants with relatively large effects[16]. We found no evidence for such synthetic associations in the joint association analysis i.e. all low frequency variants were either lead variants or jointly associated (q-value < 0.05) with a common variant. Twenty of the 202 FDR variants (9.9%) were indels (4–14 bp size) compared to 8.8% of all the variants in the meta-analysis (P = 0.60). Low frequency variants (MAF < 0.05) were strikingly underrepresented (6.9% vs. 29.0%; P = 4.9 × 10−12) which may reflect on the statistical power to detect modest effects associated with these variants (vide infra).

Annotation and ENCODE analysis

Functional annotations were assigned to the 9.4 million variants studied in the CAD additive meta-analysis using the ANNOVAR software[18] (Supplementary Table 8). The 202 FDR variants were depleted in intergenic regions (P = 2.5 × 10−7) and enriched in introns (P = 0.00035). Variants were also assigned to three sets of ENCODE (Encyclopedia of DNA Elements) features, namely histone/chromatin modifications (HM), DNase I hypersensitive sites (DHS) and transcription factor binding sites (TFBS) (Supplementary Table 9). The FDR variants showed independent enrichment across 11 cell types for the HM (P = 2.8 × 10−6) and DHS (P = 0.0003) ENCODE feature sets and with genic annotation status (P = 0.0013) (Supplementary Table 10 and Supplementary Table 11). These associations were also evident in three cell types selected with maximal CAD relevance with a 2.6-fold enrichment for DHS, 2.2-fold enrichment for HM and 1.6-fold enrichment for genic status (Supplementary Table 12 and Supplementary Table 13). These findings suggest that the 202 FDR variants are enriched for functional variants with potential relevance to CAD pathogenesis.

Post hoc power calculations

8.2M of 9.4M (87%) of the analysed variants were highly powered (> 90%) to detect an OR ≥ 1.3 (Supplementary Table 7). The number of variants with power ≥ 90% to detect associations varies systematically with allele frequency and imputation quality (Supplementary Fig. 4 for OR = 1.3); 1.5M of 2.7M (55%) low frequency variants (0.005 < MAF < 0.05) in the meta-analysis were adequately powered to detect an OR ≥ 1.3 as most of these variants were accurately imputed (median imputation quality 0.94, interquartile range 0.88, 0.98). With more common variants (MAF > 0.05) almost all (99.8%) were highly powered to detect OR ≥ 1.3. However, in terms of the total coverage of low frequency variation, only 15.3% of 9.3M low frequency (0.005 < MAF < 0.05) variants in the 1000 Genomes phase 1 v3 training set surpassed the allele frequency and imputation quality entry criteria in the 60% of the studies required for inclusion in the meta-analysis and were predicted to be adequately powered to detect significant associations; 100% of these variants were highly powered (>90%) to detect OR ≥ 3.15.

Interrogation of 10 novel additive and recessive loci

We examined whether there were any eQTLs, associations with known cardiovascular risk factors or prior evidence of involvement of genes with atherosclerotic processes in each of the newly identified loci to define putative mechanisms by which the loci may affect risk of CAD.

Chromosome 4q12 (REST – NOA1) locus: The lead SNP rs17087335 lies within an intron of the nitric oxide associated 1 (NOA1) gene; 23 SNPs in linkage disequilibrium (LD r2>0.8) show CAD associations (P < 1 × 10−6) across the NOA1 and the repressor element-1 silencing transcription factor (REST) genes (Fig. 4a). NOA1 is a GTP binding protein involved in the regulation of mitochondrial respiration and apoptosis[19]. REST is a transcription factor that suppresses the expression of voltage dependent sodium and potassium channels[20] and has been shown to maintain vascular smooth muscle cells in a quiescent non-proliferative state and is itself down-regulated in neointimal hyperplasia[21]. SNP rs17087335 shows a cis-eQTL signal with REST in lung[22] (Supplementary Table 14).

Figure 4

Regional association plots of the eight additive (a–h) and two recessive (i–j) novel CAD loci. The association statistics have been adjusted for over-dispersion following meta-analysis (genomic control parameter 1.18 for the additive and 1.05 for the recessive models). Linkage disequilibrium (r2) calculations were based on the combined 1000 Genomes phase 1 v3 training dataset. Genomic coordinates in mega-base pairs (Mb) refer to the hg19 sequence assembly.

Chromosome 7q36.1 (NOS3) locus: The lead SNP rs3918226 (MAF = 0.07) lies in the first intron of the nitric oxide synthase 3 (NOS3) gene (Fig. 4b). This SNP has been tentatively associated with CAD (OR = 1.14, P = 1.4 × 10−4) in a candidate gene meta-analysis based on 15.6K cases and 35K controls genotyped with the HumanCVD BeadChip[7] and firmly associated with essential hypertension (OR = 1.34, P = 1.0 × 10−14)[23]. NOS3 is involved in production of nitric oxide (NO), a potent vascular smooth muscle relaxant and a well-studied candidate gene for CAD. Indeed, proteins forming the NO receptor (soluble guanylyl cyclase) display both linkage as well as genome-wide association with CAD[3,24]. There are several overlapping ENCODE features in the NOS3 intron 1 suggesting a functional role for rs3918226. However, there are 30 genes neighbouring NOS3 within a ±1 centiMorgan window and the current data does not allow candidacy of these genes to be excluded. A non-synonymous SNP, rs1799983, in NOS3 previously associated with cardiovascular phenotypes[25] is in weak LD with rs3918226 but did not achieve significance in an additive or joint association analysis.

Chromosome 11p15.4 (SWAP70) locus: SNP rs10840293 is intronic to the switch-associated protein-70 gene (SWAP70) (Fig. 4c). SWAP-70 is a signalling molecule involved in the regulation of filamentous-actin networks[26] in cell migration and adhesion. SNP rs10840293 with other SNPs in strong LD are cis-eQTLs in naïve and challenged monocytes[27], SNP rs93138 shows a strong association with CAD (P = 5.5 × 10−8) and is a cis-eQTL in naïve and challenged monocytes[28], fat[29], skin[29] and lung[22] (Supplementary Table 14); three of these SNPs (rs93138, rs173396 and rs472109) are intronic and lie within ENCODE regulatory functional elements. Although this CAD-associated locus includes 33 genes, the eQTL and ENCODE data implicate SWAP70 as a plausible causal gene and suggest putative causal SNPs.

Chromosome 15q22.33 (SMAD3) locus: The lead SNP rs56062135 is intronic to SMAD3 and the CAD association is tightly localized between two recombination hot spots (Fig. 4d). Mice lacking Smad3, a major downstream mediator of TGF-β, show enhanced neointimal hyperplasia with decreased matrix deposition in response to vascular injury[30]. SMAD3 had been tentatively associated in an earlier CAD GWAS[31], although that lead SNP (rs17228212) is in linkage equilibrium with rs56062135, showed modest association (P = 0.009) in the present GWAS and no evidence of joint association (Supplementary Table 6).

Chromosome 15q26.1 (MFGE8 – ABHD2) locus: The lead intergenic SNP rs8042271 maps 117kb upstream of milk fat globule-EGF factor 8 (MFGE8) and 57kb upstream from abhydrolase domain-containing protein 2 (ABHD2) (Fig. 4e). MFGE8 (lactadherin) has a crucial role in VEGF dependent neovascularization[32] and is secreted from activated macrophages and binds to apoptotic cells, facilitating phagocytic engulfment[33]. ABHD2[34] has been shown to be expressed in human atherosclerotic lesions with higher levels in patients with unstable angina. There were no overlapping risk factor QTL, eQTL or ENCODE features to guide the nomination of a putative causal gene.

Chromosome 17q23.2 (BCAS3) locus: The lead intronic SNP rs7212798 lies in the breast carcinoma amplified sequence 3 (BCAS3) gene (Fig. 4f). Multiple variants in LD with rs7212798 map to BCAS3 introns and show strong association with CAD. BCAS3 encodes the Rudhira protein that has been shown to activate Cdc42 to affect actin organization and control cell polarity and motility in endothelial cells thus contributing to angiogenesis[35].

Chromosome 18q21.32 (PMAIP1 – MC4R) locus: The lead intergenic SNP rs663129 lies 266 kb downstream of the phorbol-12-myristate-13-acetate-induced protein 1 (PMAIP1) gene and 200kb downstream of the melanocortin 4 receptor (MC4R) gene (Fig. 4g). PMAIP1 is a HIF1A-induced proapoptotic gene that mediates hypoxic cell death by the generation of reactive oxygen species[36]. MC4R is a well-studied obesity-related locus and the variant (and corresponding proxy variants) that were associated with higher CAD risk were also associated with BMI (P = 6×10−42) and obesity-associated risk factors including higher triglyceride and lower HDL concentrations and type 2 diabetes[37-41]. However, we found no eQTL data or ENCODE features for the lead or proxy SNPs to further implicate MC4R as the causal gene underlying the CAD susceptibility.

Chromosome 22q11.23 (POM121L9P – ADORA2A) locus: The lead SNP rs180803 lies in the non-coding RNA POM121 transmembrane nucleoporin-like 9, pseudogene (POM121L9P). A ±1 centiMorgan region spans 1.2Mb and includes 21 variants that are genome-wide significantly associated with CAD, most of them in LD (r2>0.6) with the lead SNP and mapping to intronic regions of the SPECC1L and ADORA2A genes (Fig. 4h).

Chromosome 12q24.23 (KSR2) locus: The lead SNP rs11830157 (MAF=0.36) associates with CAD risk in a recessive model (genotypic OR = 1.12) and is intronic to the kinase suppressor of ras 2 (KSR2) gene (Fig. 4i) and overlaps with ENCODE functional elements. KSR2 interacts with multiple proteins including AMPK and rare loss-of-function coding variants in KSR2 are associated with severe obesity, hyperphagia and insulin resistance, a phenotype recapitulated in Ksr2 null mice[42].

Chromosome 19q13.11 (ZNF507 – LOC400684) locus: The lead SNP rs12976411 (MAF=0.09) lies in an uncharacterized non-coding RNA (LOC400684) and is 3.4 kb 3’ of ZNF507 (Fig. 4j). The minor allele shows a protective CAD effect (genotypic OR = 0.69) in the recessive model. ENCODE analysis of this locus suggests that several SNPs including rs12981453 and rs71351160, which are in strong LD (r2>0.8) and intronic to ZNF507 overlap with ENCODE functional elements.

Discussion

We demonstrate that the ability of GWAS to investigate the genetic architecture of complex traits is enhanced by the 1000 Genomes Project. This has allowed us to conclude that low frequency variants of larger effect, synthetic associations and insertion/deletion polymorphisms are unlikely to explain a significant portion of the missing heritability of CAD. Rather, all ten novel CAD associations identified in the present analysis, as well as all but one of the previously identified loci, are represented by risk alleles with a frequency of greater than 5%. Thus, this comprehensive analysis strongly supports the common disease/common variant hypothesis[43] given that it was powered to detect MAF < 0.05 variants with OR > 1.5. Moreover, risk alleles are significantly clustered within or close to genes and enriched in regions with functional annotations. Finally, genes implicated by this unbiased approach suggest hypotheses that explore the biology of the arterial vessel wall as a critical component of CAD pathogenesis.

The success of the GWAS meta-analysis strategy to map common, small-effect susceptibility variants for complex diseases has leant heavily on genotype imputation with publically available training sets. The 1000 Genomes Project provides a substantial step-up from the HapMap era in terms of coverage of lower frequency variants and the integration of insertion/deletion polymorphism (Fig. 1). Lead SNPs for 4 of the 10 novel CAD loci were either absent or imperfectly tagged (r2 < 0.8) in the HapMap2 training set, which reduced the power of discovering these loci in previous GWAS meta-analyses. Although lower frequency variants often show geographical differentiation[5], the phase1 v3 training set includes numerous low MAF variants that are tractable to a global meta-analysis that includes multiple continental ancestry groups. Key SNPs in apolipoproteins E and (a) and PCSK9, which mediate their CAD-effects via LDL-cholesterol linked mechanisms showed strong associations and reinforces the sensitivity of our 1000 Genomes analysis to detect lower frequency, imperfectly imputed susceptibility variants which were missed in HapMap-based GWAS.

Association analysis under the customary additive inheritance model widely used in GWAS analyses is optimally powered to detect traits with no dominance variance but conveniently has adequate power to detect dominantly inherited traits as well[44]. However, the additive model is systematically underpowered to detect recessively inherited traits particularly for lower frequency alleles[44]. This motivated our recessive meta-analysis, which revealed two novel CAD loci, KSR2 and ZNF507-LOC400684 that escaped detection in a conventional additive association scan.

Our GWAS explores two potential sources of missing heritability for CAD as it includes indels and an extended panel of lower frequency variants. Although there was no evidence that indels were systematically enriched for CAD associations, they were present in 10% of the 202 variants with an FDR q-value < 5%. In terms of surveying the totality of human genetic variation, the 1.5M of 2.7M lower frequency variants that were included in the meta-analysis with power to detect moderate-penetrance alleles (OR > 1.3) might seem modest. Yet the relative paucity of significant associations, 15 variants with MAF < 0.05 that explain 2% of CAD heritability and provide no evidence of synthetic associations, will temper expectations for the role of low frequency variants underlying CAD susceptibility, specifically with respect to risk prediction in a population-based setting. It is though important to acknowledge that GWAS based on SNP array data has limited power to resolve genes carrying rare mutation burdens. For example, LDLR[45], APOA5[45], APOC3[46] and NPC1L1[47] are loaded with susceptibility or protective mutations for CAD. These discoveries were only revealed by whole-exome sequencing studies in large series of cases and controls yet explain less than 1% of the missing heritability of CAD[45].

Annotation analysis showed that the CAD associated variants were significantly clustered within or close to genes. Furthermore, there was a strong and independent enrichment for an overlap of the CAD associations with ENCODE features, particularly in cell types relevant to CAD pathogenesis. This phenomenon has been previously reported for other diseases and traits[48] and can guide candidate gene nomination and the design of future functional studies. We found few suggestions of overlap with risk factor QTL or eQTL in available datasets; this may in part reflect that the use of proxy variants can be limiting in cross-referencing 1000 Genomes and HapMap association databases.

Coronary atherosclerosis underlies the development of the vast majority of cases of MI; therefore the two are intimately related. However, additional factors such as plaque vulnerability or the extent of the thrombotic reaction to plaque disruption may predispose to MI in the presence of CAD[49]. We confirmed that ABO is particularly associated with risk of MI[50] suggesting that this locus may specifically increase the risk of plaque rupture and/or thrombosis. In contrast, HDAC9 showed a stronger association with CAD than with MI suggesting that it might predispose to atherosclerosis but not the precipitant events leading to an MI. However, HDAC9 shows even stronger associations with ischemic strokes involving thrombosis or embolism due to atherosclerosis of a large artery[51]. Although further epidemiological as well as experimental data are required to substantiate these findings, they suggest that certain loci may affect distinct mechanisms related to the development and progression of CAD.

Several of the genes implicated to date in large-scale analyses of CAD susceptibility encode proteins with a known role in the biology of circulating risk factors for CAD, notably lipid levels and the metabolism of lipoproteins; others relate to other known atherosclerosis risk factors, such as genes implicated in systemic inflammation and in hypertension. This is unsurprising, partly because of the undoubted importance of these known risk factors in the etiology of CAD but also because some of these prior analyses particularly targeted genes involved in risk factor traits, for example the HumanCVD BeadChip[52] was based on candidate genes and the Metabochip[3,53] drew on prior association data with risk factor traits as well as an earlier HapMap2 CAD GWAS meta-analysis[54]. The current experiment adopts a completely unbiased approach and is the first to do so at very large scale. In this respect, it is notable that for some of the novel loci where genomic data, biological precedent and eQTL associations suggest a plausible novel CAD gene, the genes so implicated have well documented roles in vessel wall biology. Their gene products are involved in diverse processes including cell adhesion and leukocyte and vascular smooth muscle cell (VSMC) migration (SWAP70[26] and ABHD2[55]), VSMC phenotypic switching (REST[20]), TGF-β signaling (SMAD3[56,57]), anti-inflammatory and infarct sparing effects (ADORA2A[58] and MFGE8[59]), angiogenesis (BCAS3[35]) and nitrous oxide signaling (NOS3[24]).

It is important to note that these putative new susceptibility genes require substantial further investigation and validation before firm vascular biology links can be established. A number of preventative strategies target the vessel wall (control of blood pressure, smoking cessation) but the large majority of existing drug treatments for lowering CAD risk operate through manipulating circulating lipids and few directly target vessel wall processes. Detailed investigation of new aspects of vessel wall biology that are implicated by genetic association, but have not previously been explored in atherosclerosis, may provide new insights into its complex aetiology and hence new targets.

Online Methods

Association Analysis

Three models of heritable disease susceptibility were analysed by logistic regression 1) an additive model where the log(genotype risk ratio) (log(GRR)) for a genotype was proportional to the number of risk alleles, 2) a recessive model where the log(GRR) for homozygotes for the minor allele was compared with a reference risk in pooled heterozygotes and homozygotes for the major allele and 3) a dominant model where the log(GRR) for homozygotes for the minor allele pooled with heterozygotes was compared with a reference risk in homozygotes for the major allele. Minor and major alleles were identified by reference to the allele frequencies in the pooled populations (i.e. all continents) of 1000 Genomes phase 1 v3 data. For the recessive and dominant analyses, genotype probabilities to allow for variable imputation quality were analysed by all contributing studies; for the additive analysis, genotype probabilities or allelic dosages were used (Supplementary Table 1).

Data Quality Control

Association data for each contributing study were individually filtered by MAF > 0.005 (estimated in combined cases and controls) and an imputation quality metric, rsq > 0.3 for minimac or info_proper > 0.4 for IMPUTE2[61]. Allele frequencies for each study were binned and compared with other studies to detect systematic flipping of alleles (Supplementary Fig. 5). Over-dispersion of association statistics was assessed by the genomic control method[62] (Supplementary Table 15) and adjusted values were submitted for meta-analysis. Variants that were retained in at least 60% of the studies were submitted for meta-analysis using the GWAMA program[63]. Following an inverse-variance weighted fixed-effects meta-analysis, heterogeneity was assessed by Cochran’s Q statistic[64] and the I2 inconsistency index[65] and variants showing marked heterogeneity were reanalysed using a random-effects model[66]. Over-dispersion in the resulting meta-analysis statistics was adjusted for by a second application of the genomic control procedure (Supplementary Fig. 6).

FDR estimation

The false discovery rate (FDR) was assessed using a step-up procedure encoded in the qqvalue Stata program[67]. This procedure has been reported to be well controlled under positive regression-dependency conditions[68]; simulations based on 1000 permuted replicates of the PROCARDIS imputed data demonstrated that the FDR was conservatively controlled (theoretical q-value = 0.05, empirical q-value = 0.026 95%CI 0.017 – 0.038) in the context of the linkage disequilibrium patterns prevalent in the 1000 Genomes phase 1 v3 training set.

GCTA and heritability analysis

Joint association analysis of the CAD additive meta-analysis results was performed using the GCTA software[17] that fits an approximate multiple regression model based on summary association statistics and linkage disequilibrium information derived from a reference genotype database (here the 1000 Genomes phase 1 v3 training set for all continents/populations that includes genotypes for 1,092 individuals). In this analysis, the lead variant is not necessarily retained in the final joint association model in situations where there might be multiple associated variants in strong LD. The accuracy of this analysis depends on appropriate ancestry matching as well as the sample size of the reference genotype panel to ensure that estimated LD correlations are unbiased and acceptably precise[69]. Simulations suggest that the expected correlation between p-values based on the GCTA method using a reference panel of 1,000 genotyped samples and p-values from an “exact” multiple regression based on experimental genotypes will be between 0.90 and 0.95[69]. We investigated the empirical accuracy of the GCTA joint association analysis by comparing GCTA joint association results with a standard multiple logistic regression analysis in 4 contributing studies (Supplementary Fig. 7). This showed that 95% of the betas and SE were accurately approximated. The −log10 p-values from the two analyses were positively correlated (0.86 < rho < 0.93) with the GCTA method showing an insignificant trend (P > 0.20) to yield slightly inflated values.

Heritability calculations were based on a multifactorial liability threshold model[70] assuming that the disease prevalence was 5% and that the total heritability of CAD was 40%[3]; multiple regression estimates of allelic effect sizes were used following the GCTA joint association analysis. Standard errors of the heritability estimates were estimated by Monte Carlo sampling with 1,000 replicates (i.e. for each variant, draw effect sizes (beta) randomly from the variants’ beta±SE estimate, calculate heritability for each beta-by-replicate draw, sum heritability across N variants within each replicate and finally calculate the standard error of the heritability estimates across the 1,000 replicates).

Power Calculations

The power to detect genetic associations depends on the magnitude of the genetic risk (i.e. effect size), the type 1 error rate, the risk allele frequency and imputation quality and the sample size. Non-centrality parameter calculations were based on double genomic controlled standard error estimates from the additive model meta-analysis; these estimates integrate information on allele frequency, imputation quality and sample size, which typically vary across studies. The type 1 error was set at 5 × 10−8 and an additive risk model was assumed.

Risk Factor QTL survey

The 10 novel CAD associated loci were scanned for associations with CAD heritable risk factors using publically available resources including large-scale GWAS consortia data downloads[37-41,71-73] and the NHGRI GWAS catalogue accessed May 2014[74]. As previous risk factor GWAS were mainly based on HapMap2 imputed datasets, all SNPs in LD (r2>0.8 based on 1000 Genomes phase 1 v3 ALL reference panel) with the novel variants were examined for risk factor associations. The novel associated loci were cross-referenced with known cis and trans eQTL associations from the University of Chicago eQTL browser (accessed July 2014), the GTEx Portal (accessed June 2014), the Geuvadis Data Browser (accessed June 2014) and other published data[22,28,29,75-79].

Annotation and ENCODE analysis

Variants were annotated using the ANNOVAR[18] (version Aug 2013) software based on a GRCh/hg19 gene annotation database. Upstream/downstream status was assigned to variants that mapped ≤1kb from the transcript start/end. Variants without intergenic annotation were assigned a genic annotation status (42%). Supplementary Table 8 shows the annotation status of 9.4M variants included in the CAD additive meta-analysis; 86% of the genic variants map to introns.

ENCODE features were downloaded from the Ensembl database using the Funcgen Perl API module release 75. The list of the ENCODE experiments stored in the Ensembl database can be browsed at http://Feb2014.archive.ensembl.org/Homo_sapiens/Experiment/Sources?db=core;ex=project-ENCODE-. This summarized 100 different types of functional evidence in 11 different cell types, a total of 379 ENCODE experiments that revealed 6,099,034 features. Variants that overlaid one or more of these features were cross-tabulated with their ANNOVAR annotation status (Supplementary Table 10); 50% of variants mapped to one or more ENCODE features and variants in ENCODE features were strongly enriched for genic annotation status. Variants were grouped into three functional sets, histone/chromatin modifications (HM), DNase I hypersensitive sites (DHS) and transcription factor binding sites (TFBS) (Supplementary Table 9). Cell types were grouped into CAD relevant and others (Supplementary Table 12) based on their potential roles in CAD pathophysiology; hepatocytes (e.g. lipid metabolism[80]), vascular endothelial cells (atherosclerosis[81]) and myoblasts (injury/repair[82]) were selected as being most relevant to the CAD phenotype. Multi-way contingency tables reporting ENCODE feature and ANNOVAR annotation status with inclusion in the FDR < 5% variant list (FDR202 status) are summarized for 11 ENCODE cell types in Supplementary Table 11 and for the three CAD relevant cell types in Supplementary Table 13. Contingency table counts were modelled by a logistic multiple regression model predicting FDR202 status with independent explanatory variables HM, DHS, TFBS and genic/intergenic status. The ENCODE[83] project has previously mapped 4,492 GWAS significant SNPs from the NHGRI (June 2011) catalogue[74] to TF (12%) and DHS (34%) features in an extended dataset of 1,640 experiments. The 202 FDR variants were slightly less prevalent in these feature groups (10.4% TF and 19.8% DHS) which could reflect a CAD-specific issue or a more general consequence of our analysis being based on a subset of the ENCODE data retrieved from the Ensembl database.