Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture.

Bone mineral density (BMD) is the most widely used predictor of fracture risk. We performed the largest meta-analysis to date on lumbar spine and femoral neck BMD, including 17 genome-wide association studies and 32,961 individuals of European and east Asian ancestry. We tested the top BMD-associated markers for replication in 50,933 independent subjects and for association with risk of low-trauma fracture in 31,016 individuals with a history of fracture (cases) and 102,444 controls. We identified 56 loci (32 new) associated with BMD at genome-wide significance (P < 5 × 10(-8)). Several of these factors cluster within the RANK-RANKL-OPG, mesenchymal stem cell differentiation, endochondral ossification and Wnt signaling pathways. However, we also discovered loci that were localized to genes not known to have a role in bone biology. Fourteen BMD-associated loci were also associated with fracture risk (P < 5 × 10(-4), Bonferroni corrected), of which six reached P < 5 × 10(-8), including at 18p11.21 (FAM210A), 7q21.3 (SLC25A13), 11q13.2 (LRP5), 4q22.1 (MEPE), 2p16.2 (SPTBN1) and 10q21.1 (DKK1). These findings shed light on the genetic architecture and pathophysiological mechanisms underlying BMD variation and fracture susceptibility.

Osteoporosis is a disease characterized by low bone mass and microarchitectural deterioration of bone tissue leading to increased risk of fracture. The disease accounts for approximately 1.5 million new fracture cases each year representing a huge economic burden on health care systems, with annual costs estimated to be $17 billion in the USA alone and expected to rise 50% by the year 2025.[1] Osteoporosis is defined clinically through the measurement of bone mineral density (BMD), which remains the single best predictor of fracture.[2,3]

Twin and family studies have shown that 50%–85% of the variance in BMD is genetically determined.[4] Osteoporotic fractures are also heritable by mechanisms that are partly independent of BMD.[5] Over the past 5 years, genome-wide association studies (GWAS) have revolutionized the understanding of the genetic architecture of common, complex diseases.[6] This strategy is providing key insights into the mechanisms of disease with prospects of designing effective strategies for risk assessment and development of new interventions.[7]

Previous GWAS have identified to-date 24 loci which influence BMD variation.[8-14] While several variants in these BMD loci have also been nominally-associated with fracture risk[15,16], none have shown robust association at genome-wide significant levels (P<5×10−8). We report here the results of the largest effort to date searching for BMD loci in >80,000 subjects and testing them for association with fracture in >130,000 cases and controls. In addition, we employed bioinformatics tools and gene expression analyses to place the identified variants in the context of pathways relevant to bone biology.

Results

This study was performed across three main stages (Fig. 1): 1) discovery of BMD loci, 2) follow-up replication and 3) association of the BMD loci with fracture.

Figure 1

Description of Study design

Stage 1: Meta-analysis of 17 Genome-wide association studies for BMD. Stage 2: 96 top independent SNPs (82 autosomal SNPs with P<5×10−6 + 5 X-chromosome SNPs + 9 SNPs from conditional analysis) were followed-up for de-novo and in-silico replication of the BMD association in 34 studies. Stage 3: the same 96 SNPs were tested for association with fracture in 50 studies with de-novo and in-silico data.

Discovery of BMD Loci (Stage 1)

We first performed a meta-analysis of GWAS for BMD of the femoral neck (FN-BMD; n=32,961) and lumbar spine (LS-BMD; n=31,800) including ~2.5 million autosomal genotyped or imputed SNPs from 17 studies from populations across North America, Europe, East Asia and Australia, with a variety of epidemiological designs and patient characteristics (see Online Methods). We also performed meta-analysis in men and women separately to identify sex-specific associations. The quantile-quantile (Q-Q) plots of the discovery meta-analysis displayed strong (and not early) deviation of the observed statistics from the null distribution for both BMD traits (Supplementary Fig. 1). After double Genomic Control (GC) correction of the overall (λFN_POOLED=1.112; λLS_POOLED=1.127) and sex-stratified analyses (λFN_WOMEN=1.091; λFN_MEN=1.059; λLS_WOMEN=1.086; λLS_MEN=1.061), SNPs in 34 loci surpassed GWS level while a total of 82 loci were associated at P<5×10−6 (Supplementary Fig. 2–3). Thirty eight loci were associated with FN-BMD, 25 with LS-BMD and 19 with both. The overlap reflects the correlation between the femoral neck and lumbar spine measurements (Pearson correlation = 0.53). Of these 82 loci, 59, 18 and 5 were prioritized from the analysis in the sex-combined, women and men sample sets, respectively (Supplementary Table 1). The meta-analysis was extended to include the evaluation of 76,253 X-chromosome imputed markers across 14 of the discovery GWAS including 31,801 participants (see Online Methods). Five X-chromosome loci were associated at P<5×10−5 four of which were derived from the sex-combined analysis and one from the analysis in men only (Supplementary Table 1). We further performed genome-wide conditional analyses in all sex-combined stage 1 studies. Each study repeated the GWAS analysis but additionally adjusted for 82 SNPs representing the autosomal loci associated at P<5×10−6 (see Online Methods). We then meta-analyzed these studies in the same way as for the primary GWA study meta-analysis. Nine loci showed at least two independent association signals arising from this conditional analysis (Supplementary Fig. 4 and Supplementary Table 2) suggesting that allelic heterogeneity is underlying BMD variation. We also assessed all possible pairwise interactions of the 82 SNPs, but none were significant after adjusting for the number of tests (Supplementary Fig. 5 and Supplementary Table 3). A total of 96 independent SNPs (82 autosomal SNPs with P<5×10−6 + 9 SNPs from conditional analysis and 5 X-chromosome SNPS) from 87 genomic loci were selected for further replication (Fig. 1).

Follow-up replication (Stage 2)

We de-novo genotyped these 96 SNPs and tested them for association with BMD in up to 50,933 additional participants from 34 studies (see Online Methods). The meta-analysis of the 96 SNPs in the discovery and replication studies (n=83,894) yielded 64 replicating SNPs from 56 associated loci. Of these loci, 32 were novel (Table 1 and Supplementary Table 4A) and 24 were reported previously[8-14] (Supplementary Table 4B). Thirty two SNPs did not reach genome-wide significance after replication (Supplementary Table 4C) and included 10 markers remaining associated at a suggestive level. Of all analyzed SNPs only one (rs9533090 mapping to 13q14.11 near RANKL) displayed high degree of heterogeneity of effects (I2>50%) across studies, despite being the marker with highest significance (P=4.82×10−68) in the fixed-effect meta-analysis (Supplementary Table 4B). After applying random effects meta-analysis, this marker was still genome-wide significant (P=3.98×10−13).

Table 1

Estimated effects of novel genome-wide significant SNPs on bone mineral density of femoral neck (FN-) BMD and lumbar spine (LS-) BMD across stages.

SNP	Locus	Closest Gene/Candidate	Functional evidence						Aa	Freqb	FN-BMD				LSBMD				Phet Sited
			eQTL	KO mouse	OMIM	Tags funct.	GRAIL	Pathway				STAGE 1 (32,961)	STAGE 2 (50,933)	STAGE 1 + STAGE 2 (83,894)		STAGE 1 (31,800)	STAGE 2 (45,708)	STAGE 1 + STAGE 2 (77,508)
											Betac	P	P	P	Betac	P	P	P
rs479336	1q24.3	DNM3							T	0.74	−0.04	1.1×10−7	1.3×10−8	8.5×10−15	−0.03	0.01	5.0×10−4	2.1×10−5	0.05
rs7584262	2p21	PKDCC		•					T	0.23	0.03	1.4×10−7	3.4×10−4	1.3×10−9	0.01	0.13	0.28	0.07	0.01
rs17040773	2q13	ANAPC1	•						A	0.76	0.03	4.3×10−6	6.1×10−5	1.5×10−9	0.01	0.61	0.21	0.19	5.2×10−3
rs1878526	2q14.2	INSIG2							A	0.22	0.00	0.70	0.97	0.79	0.04	7.3×10−6	3.4×10−6	1.2×10−10	8.6×10−5
rs1026364	3q13.2	KIAA2018				•			T	0.37	0.03	2.0×10−6	2.5×10−5	4.1×10−10	0.02	0.04	7.3×10−3	7.6×10−4	0.11
rs344081	3q25.31	LEKR1							T	0.87	0.03	1.1×10−4	2.5×10−3	2.2×10−6	0.06	2.8×10−5	3.5×10−8	4.5×10−12	0.12
rs3755955	4p16.3	IDUA		•		•			A	0.16	−0.05	3.9×10−7	6.1×10−9	1.5×10−14	−0.05	1.4×10−7	5.5×10−9	5.2×10−15	0.80
rs11755164	6p21.1	SUPT3H/RUNX2		•			•	•	T	0.40	−0.01	0.23	0.12	0.05	−0.03	3.5×10−7	9.2×10−6	5.6×10−11	2.1×10−3
rs9466056	6p22.3	CDKAL1/SOX4							A	0.38	−0.03	1.8×10−8	1.6×10−6	2.7×10−13	−0.03	6.5×10−5	1.1×10−4	3.6×10−8	0.34
rs3801387	7q31.31	WNT16				•	•	•	A	0.74	−0.08	4.2×10−14	2.0×10−27	5.0×10−40	−0.10	1.4×10−16	1.5×10−36	3.2×10−51	0.09
rs13245690 g	7q31.31	C7orf58							A	0.62	0.00	8.6×10−5	0.69	8.2×10−4	0.03	1.1×10−9	1.3×10−3	6.0×10−11	0.05
rs7812088	7q36.1	ABCF2							A	0.13	0.04	1.2×10−6	4.4×10−4	7.3×10−9	0.04	2.9×10−5	1.1×10−3	2.2×10−7	0.86
rs7017914e	8q13.3	XKR9/LACTB2	•			•			A	0.49	0.02	4.7×10−8	7.1×10−3	1.9×10−8	−0.01	0.35	0.41	0.98	9.1×10−5
rs7851693	9q34.11	FUBP3							C	0.64	0.05	3.1×10−8	1.4×10−15	3.4×10−22	0.04	0.06	6.7×10−8	6.1×10−8	0.02
rs3905706	10p11.23	MPP7							T	0.22	−0.02	0.63	1.7×10−3	0.03	0.05	2.9×10−9	6.7×10−9	2.4×10−16	5.8×10−11
rs1373004	10q21.1	MBL2/DKK1		•				•	T	0.13	−0.04	1.4×10−5	1.5×10−4	1.5×10−8	−0.05	5.4×10−8	2.2×10−6	1.6×10−12	0.28
rs7071206	10q22.3	KCNMA1							T	0.78	0.01	0.29	0.26	0.81	−0.05	1.5×10−12	6.2×10−9	5.0×10−19	5.9×10−9
rs7084921	10q24.2	CPN1							T	0.39	0.03	1.4×10−4	1.6×10−6	9.0×10−10	0.03	0.01	1.9×10−5	9.2×10−7	0.58
rs10835187	11p14.1	LIN7C	•						T	0.55	−0.01	0.17	0.08	0.03	−0.02	3.0×10−5	2.4×10−4	4.9×10−8	0.03
rs7953528	12p11.22	KLHDC5/PTHLH		•	•		•	•	A	0.18	0.04	5.8×10−8	2.4×10−6	1.9×10−12	−0.02	0.94	0.05	0.13	2.3×10−7
rs2887571	12p13.33	ERC1/WNT5B					•	•	A	0.76	−0.03	1.1×10−4	1.6×10−5	6.5×10−9	−0.04	2.2×10−7	2.9×10−6	5.6×10−12	0.37
rs12821008	12q13.12	DHH						•	T	0.39	0.03	1.9×10−4	5.2×10−4	3.3×10−7	0.05	1.5×10−7	1.9×10−9	1.2×10−15	0.06
rs1053051	12q23.3	C12orf23	•						T	0.52	−0.03	1.4×10−5	1.8×10−5	9.6×10−10	−0.02	2.5×10−6	2.4×10−3	7.9×10−8	0.76
rs1286083	14q32.12	RPS6KA5							T	0.81	−0.05	2.9×10−8	9.3×10−9	2.0×10−15	−0.04	1.7×10−11	7.1×10−6	1.8×10−14	0.92
rs4985155	16p13.11	NTAN1	•			•			A	0.67	−0.03	3.5×10−4	1.4×10−7	1.7×10−10	−0.03	8.7×10−7	1.8×10−4	2.2×10−9	0.98
rs9921222	16p13.3	AXIN1		•			•	•	T	0.48	−0.03	2.5×10−7	2.4×10−6	5.2×10−12	−0.04	2.2×10−8	8.3×10−10	1.0×10−16	0.26
rs13336428	16p13.3	C16orf38/CLCN7		•	•				A	0.43	−0.04	2.9×10−7	1.1×10−10	1.5×10−16	−0.04	5.9×10−5	5.8×10−10	1.7×10−13	0.80
rs1566045	16q12.1	SALL1/CYLD		•					T	0.80	−0.06	5.0×10−12	3.0×10−12	1.9×10−22	−0.01	7.8×10−3	0.55	0.04	7.5×10−6
rs1564981 g	16q12.1	CYLD		•					A	0.47	−0.02	1.1×10−3	0.01	4.4×10−5	−0.03	6.2×10−8	5.4×10−4	2.0×10−10	0.50
rs4790881	17p13.3	SMG6							A	0.69	0.05	1.7×10−8	1.2×10−11	9.8×10−19	0.04	6.0×10−4	1.7×10−6	3.4×10−9	0.13
rs7217932	17q24.3	SOX9		•	•		•	•	A	0.46	0.03	3.7×10−8	2.7×10−5	1.9×10−11	0.01	0.31	0.15	0.08	3.8×10−3
rs4796995	18p11.21	C18orf19							A	0.63	0.02	3.2×10−6	1.1×10−3	4.9×10−8	0.01	5.2×10−4	0.11	6.7×10−4	0.29
rs10416218	19q13.11	GPATCH1				•			T	0.73	−0.02	5.7×10−6	7.1×10−4	5.5×10−8	−0.03	9.2×10−9	1.2×10−4	6.6×10−11	0.38
rs5934507f	Xp22.31	FAM9B/KAL1							A	0.74	−0.08	0.01	8.3×10−4	1.6×10−4	−0.09	5.7×10−6	3.2×10−4	1.2×10−8	0.17

Boldface indicates P < 5 × 10−8 or Site-specificity P<5×10−4.

The effect estimates (Beta) are expressed as standardized values per copy of the SNP allele (A)from fixed effects meta-analysis.

Freq is the allele frequency of A.

Effect estimates were calculated in the Stage 2 samples.

Site-specificity null hypothesis, i.e., βLS = βFN.

rs7017914 was discovered in the meta-analysis of women-only, the effects and p-value for this marker are for the meta-analysis of women samples.

rs5934507 was discovered in the meta-analysis of men-only, the effects and p-value for this marker are for the meta-analysis of men samples.

The snps rs13245690 and rs1564981 were independently associated to their main signals in conditional analysis. Black dots in the six functional evidence columns indicate that: first: SNP is an expression quantitative trait loci (eQTL); second: there is a knockout mouse with skeletal phenotypes (MGI database 2011); third: candidate gene has been involved in a monogenic syndrome with skeletal phenotypes (OMIM database 2011); fourth: most significant SNP tags a SNP predicted to have impact on function of the candidate gene; fifth: gene is best candidate in GRAIL analysis; and sixth: candidate gene is part of a bone active pathway.

Two of the novel loci were discovered in the sex-stratified meta-analysis: 8q13.3 in women and Xp22.31 in men; however, only the association at Xp22.31 showed significant evidence for sex-specificity as reflected by significant heterogeneity of effects across sex strata (Phet=1.62×10−8). Yet, we acknowledge that the association at 8q13.3 in women may be driven by a lower number of men in the discovery and replication datasets (Table 1 and Supplementary Table 5). Furthermore, evidence for BMD site-specificity (Phet<5×10−4) was observed in a fraction of the loci including 6 of the 32 novel and 4 of the 24 known loci (Table 1 and Supplementary Fig. 6). Among the newly identified loci, 2q14 (INSIG2), 12p11.22 (PTHLH) and 16q12.1 (CYLD) displayed site-specificity with FN-BMD while 8q13.3 (LACTB2), 10p11.23 (MPP7) and 10q22.3 (KCNMA1) displayed site-specificity with LS-BMD.

After replication, the conditional analysis provided significant evidence of association (P<5×10−8) in 8 of the 9 loci containing secondary signals (Supplementary Fig. 4 and Supplementary Table 2). Three loci included variants localized less than 40 Kb from the initial main signal suggesting allelic heterogeneity and included the 1p31.3 (represented by rs17482952 near WLS), 6q25.1 (rs7751941 near ESR1) and the 16q12.1 (rs1564981 near CYLD) loci. The secondary signal in 16q12.1 (rs1564981) showed a strong association with LS-BMD, while the main signal in this locus (rs1566045) was only associated with FN-BMD. The other five secondary signals were represented by variants localized at more than 180kb away from the initial main signal and contained different candidate genes including the 1p36.12 (rs7521902 near WNT4), 7p14.1 (rs10226308 near SFRP4), 7q31.31 (rs13245690 near C7orf58), 12q13.13 (rs736825 near HOXC6) and the 17q21.31 (rs4792909 near SOST) loci. The secondary signal mapping to the 13q14.11 locus (rs7326472) did not achieve genome-wide significance after replication.

Association of the BMD loci with fracture (Stage 3)

We tested the 96 markers for association with fracture in 31,016 cases and 102,444 controls from 50 studies with fracture information. This collection included: 5,411 cases and 21,909 controls tested in the BMD GWAS discovery samples, 9,187 cases and 45,057 controls tested by in-silico replication and 16,418 cases and 35,478 controls tested by de-novo genotyping (Figure 1 and Online Methods). In this fracture meta-analysis fourteen loci were significantly associated with any type of fracture at a Bonferroni level (P=5×10−4), of which five included novel BMD loci. None of the markers displayed large estimates of heterogeneity (Table 2, Supplementary Table 6 and Supplementary Fig. 7). Markers at six of these loci reached P<5×10−8 including 18p11.21 (C18orf19), 7q21.3 (SLC25A13), 11q13.2 (LRP5), 4q22.1 (MEPE), 2p16.2 (SPTBN1) and 10q21.1 (DKK1). The proportion of the overall fracture risk explained by BMD ranged between 0.09 and 0.40 across markers (Supplementary Table 7) and was estimated in a subset of Stage 2 samples (including n=8,594 cases and 23,218 controls) by modeling the BMD SNP effect on fracture risk with and without the inclusion of BMD as covariate. In general, the effect of these SNPs on BMD was larger than on fracture risk (Fig. 2A) except for the most significantly associated fracture locus 18p11.21 (Fig. 2B). SNPs in genes of the RANK-RANKL-OPG pathway (TNFRSF11A-TNFSF11-TNFRSF11B) despite being the strongest-associated BMD loci were not significantly associated with fracture. All 31 BMD loci with nominal association with fracture risk (P<0.05) showed consistent direction (decreasing BMD allele increased risk of fracture). When we performed subgroup analyses using “cleaner” phenotype definitions such as limiting to clinically-validated fractures and stratifying by anatomical site (i.e. “non-vertebral” fractures and “vertebral” fractures), we did not gain any additional signals (Supplementary Table 8). At a nominally significant level (P<0.05) only three loci were associated with vertebral fracture and all 14 BMD loci were associated with non-vertebral fracture, but these difference in effects between fracture sites were not significant. Therefore, the power of our study did not benefit from improving phenotype definition at expense of (a lower) sample size.

Table 2

Association of identified BMD loci with risk for any type of low-trauma fracture

SNP	Locus	Closest Gene/Candidate	Functional evidence						Meta-analysis without studies included in BMD discovery				Combined meta-analysis results
			eQTL	KO mouse	OMIM	Tags funct.	GRAIL	Pathway	25,605 cases, 80,535 controls				31,016 cases, 102,444 controls
									Risk Allele	Freqb	OR (95% CI)	P	OR (95%CI)	P	Qhet P	I2
Loci significantly associated with fracture risk at P < 5 × 10−8
rs4233949	2p16.2	SPTBN1				•			G	0.63	1.07 (1.04–1.09)	1.4×10−7	1.06 (1.04–1.08)	2.6×10−8	0.36	6
rs6532023	4q22.1	MEPE/SPP1	•	•			•	•	G	0.67	1.06 (1.04–1.09)	8.8×10−7	1.06 (1.04–1.09)	1.7×10−8	1.00	0
rs4727338	7q21.3	SLC25A13							G	0.32	1.08 (1.05–1.10)	1.0×10−8	1.08 (1.05–1.10)	5.9×10−11	0.03	31
rs1373004	10q21.1	MBL2/DKK1		•				•	T	0.13	1.09 (1.06–1.13)	7.2×10−7	1.10(1.06–1.13)	9.0×10−9	0.64	0
rs3736228	11q13.2	LRP5		•	•	•	•	•	T	0.15	1.09 (1.05–1.12)	2.1×10−6	1.09 (1.06–1.13)	1.4×10−8	0.78	0
rs4796995	18p11.21	C18orf19							G	0.39	1.06 (1.04–1.09)	6.4×10−7	1.08 (1.06–1.10)	8.8×10−13	0.12	20
Other significant loci associated with fracture risk at P < 5 × 10−4 (Bonferroni)
rs6426749	1p36.12	ZBTB40							G	0.83	1.06 (1.03–1.09)	2.4×10−4	1.07 (1.04–1.10)	3.6×10−6	0.07	24
rs7521902	1p36.12a	WNT4	•				•	•	A	0.27	1.10(1.06–1.14)	3.5×10−6	1.09 (1.06–1.13)	1.4×10−7	0.87	0
rs430727	3p22.1	CTNNB1		•			•	•	T	0.47	1.05 (1.03–1.08)	2.4×10−5	1.06 (1.03–1.08)	2.9×10−7	0.93	0
rs6959212	7p14.1	STARD3NL							T	0.33	1.04 (1.02–1.07)	1.0×10−3	1.05 (1.02–1.07)	7.2×10−5	0.43	2
rs3801387	7q31.31	WNT16				•	•	•	A	0.74	1.08 (1.05–1.11)	4.9×10−9	1.06 (1.04–1.08)	2.7×10−7	0.69	0
rs7851693	9q34.11	FUBP3							G	0.37	1.04 (1.01–1.06)	1.9×10−3	1.05 (1.02–1.07)	3.5×10−5	0.65	0
rs163879	11p14.1	DCDC5							T	0.66	1.06(1.03–1.09)	6.4×10−6	1.05 (1.03–1.07)	3.3×10−5	0.05	28
rs1286083	14q32.12	RPS6KA5							T	0.81	1.05 (1.02–1.08)	9.8×10−4	1.05 (1.03–1.08)	7.2×10−5	0.01	34
rs4792909	17q21.31a	SOST		•	•			•	G	0.62	1.07 (1.04–1.11)	4.0×10−5	1.07 (1.04–1.10)	6.9×10−6	0.31	10
rs227584	17q21.31	C17orf53	•						A	0.67	1.05 (1.02–1.08)	2.2×10−4	1.05 (1.03–1.07)	4.1×10−5	0.49	0

Odds Ratios (OR)estimatedper risk allele copy for any low-trauma fracture among cases as compared with controls. Q is the Cochran’s Q statistic and I2 is the measure of heterogeneity.

rs7521902 and rs4792909 are secondary independent signals.

Freq is the frequency ofthe risk allele. Boldface indicates gene names from novel loci and/or P<5×10−8. Black dots in the six functional evidence columns indicate that: first: SNP is an expression quantitative trait loci (eQTL); second: there is a knockout mouse with skeletal phenotypes (MGI database 2011); third: candidate gene has been involved in a monogenic syndrome with skeletal phenotypes (OMIM database 2011); fourth: most significant SNP tags a SNP predicted to have impact on function of the candidate gene; fifth: gene is best candidate in GRAIL analysis; and sixth: candidate gene is part of a bone active pathway.

Figure 2

A) Phenotype-wide effect for the BMD loci associated with fracture and those part of the OPG/RANK/RANKL pathway. Genetic effect estimates (± 95%CI) are shown for fracture (blue circles), lumbar-spine BMD (yellow rectangles) and femoral-neck BMD (green diamonds) for the 14 loci associated with fracture risk. Horizontal lines represent 95% confidence limits. Effect estimates are shown after transformation of the standardized mean difference(SMD) in the BMD effect to odds ratio equivalents[34] (e.g. a 0.02 SMD in the BMD effect corresponds to an OR of 1.04). Secondary signals for markers rs227584 and rs6426749 are marked with an asterisk and the signals mapping to the OPG (rs2062377), RANK (rs884205), and RANKL (rs9533090) genes are marked with a hash. B) Regional association plot for the 18p11.21 locus displaying the P-values for the top SNP associated with fracture (rs4796995) together with P-values of the BMD discovery (Stage 1) and combined with the BMD replication (Stage 1 + 2). SNPs are plotted by position in a 500Kb window of chromosome 18 against association with FN-BMD (−log10 P). Estimated recombination rates (from HapMap) are plotted in cyan to reflect the local LD structure. SNPs surrounding most significant SNP are color-coded reflecting LD between markers (pairwise r2). Genes, exons and transcription direction are derived from the UCSC genome browser.

Allele Risk Modelling for Osteoporosis and Fracture

The combined effect of all significant autosomal SNPs on BMD, osteoporosis and any type of fracture was modelled in the PERF study (n=2,836), a prospective study in postmenopausal Danish women aged 55–86 years.[17] This study comprises an independent validation setting since it was excluded from the overall meta-analysis for this purpose (see Supplementary Note for details). Risk alleles in the score (i.e., BMD-decreasing alleles) were weighted by their individual effect on BMD and grouped in 5 bins (Supplementary Table 9). The difference in mean FN-BMD between individuals in the highest bin of the risk score (9% of the population; n=244) and those in the middle bin (34% of the population; n=978) was −0.33 SDs (Fig. 3A). This analysis was based on 63 SNPs and explained 5.8% (95%CI [4.0–7.6]) of the total genetic variance in FN-BMD.

Figure 3

Combined effect of BMD-decreasing alleles and fracture risk-increasing risk alleles modelled in the population-based Prospective Epidemiological Risk Factor (PERF) study (n=2,836 women) on A) Baseline FN-BMD standardized residuals (Z-scores), B) Risk for Osteoporosis and C) Risk for Any type of fracture. The genetic score of each individual for A) and B) was based on the 63 SNPs displaying genome-wide significant association with BMD (55 main and 8 secondary signals), and for C) was based on the 16 BMD SNPs associated with fracture. Both genetic scores are weighted for relative effect sizes estimated without the PERF study. Weighted allele counts summed for each individual were divided by the mean effect size making them equivalent to the percent of alleles carried by each individual and binned into 5 categories. Histograms describe counts of individuals in each genetic score category (left axis scale). Diamonds (right axis scale) represent A) mean FN-BMD standardized levels, risk estimates in the form of odds ratio for B) Osteoporosis (defined as NHANES T-score<−2.5) and for C) Any type of Fracture using the middle category as reference (OR=1). Vertical lines represent 95% confidence limits.

The ability of this genetic score to predict the risk for osteoporosis (defined as T-score<−2.5) and for fracture was modelled in the PERF study using the middle bin as reference (OR=1). Women in the highest bin had 1.56 (95%CI [1.12–2.18]) increased odds for osteoporosis (Fig. 3B), while women in the lowest bin were protected for both osteoporosis (OR=0.38 (95%CI [0.23–0.63])). A model based on the 16 BMD SNPs associated with fracture risk showed that women in the highest bin had 1.60 (95%CI [1.15–2.24]) increased odds for fracture, while women in the lowest bin had a decreased risk for fracture (OR=0.54 (95%CI [0.36–0.83])) (Fig. 3C). Despite serving as a robust proof of principle of the relation between the BMD-decreasing alleles and the risk of osteoporosis and fracture, prediction ability was modest. The ROC analysis showed a significant but relatively small discrimination ability of the genetic score alone with an area under the curve (AUC) of 0.59 (95%CI [0.56–0.62]) for osteoporosis (Supplementary Fig. 8). Adding this score to a model with age and weight alone (AUC 0.75 (95%CI [0.73–0.77])) did not substantially increase discrimination (AUC 0.76 (95%CI [0.74–0.78])). A similar pattern was observed for fracture discrimination with an AUC of 0.57 (95%CI [0.55–0.59]) in a model with the score alone and of 0.62 (95%CI [0.60–0.64]) in a model with age, weight and height. A model considering all 63 SNPs did not change the AUC for fracture risk prediction (0.57 (95%CI [0.54–0.59])).

Functional annotations and pathway analyses

1000 Genomes Project (1KGP)

For the purpose of fine-mapping and identifying additional SNPs with putative functional implication using linkage disequilibrium (LD), a subset of nine discovery studies (FN-BMD n=21,699; LS-BMD n=20,835) used 1KGP data (Release June/2010) to re-impute genotypes contained in the 55 autosomal BMD loci (see Supplementary note for details). In 13 of the 55 BMD loci (X-chromosome SNP not included) we identified markers in a surrounding 1 Mb region that were imputed from 1KGP, and that were more significant than the original HapMap signal (Supplementary Tables 10 and 11) highlighting the benefit of using a denser reference panel of markers. All HapMap markers in LD with variants with functional annotation and displaying higher significance in the 1KGP meta-analysis are shown on Supplementary Table 12. In 14 of the 56 discovered BMD loci a marker from the HapMap imputation was highly correlated (r2 >0.8) with at least one putative functional variant annotated in the 1KG reference. Three of the 14 BMD loci associated with fracture contained putative functional variants tagged by the top SNPs of the BMD meta-analysis. These included the known rs3736228 LRP5 (Ala→Val) functional marker,[16,18] the intronic marker rs3779381 within a promoter/regulatory region of WNT16, and one intronic marker (rs4305309) within a promoter/regulatory region of SPTBN1.

eQTL analyses

Expression profiles at the GWS BMD loci were analyzed within four datasets (see Supplementary Note). In trans-iliac bone biopsies, expression of five genes correlated with LS-BMD and/or FN-BMD of the donors with P < 0.001 including PSME4 (2p16.2), DKK1 (10q21.1), C17orf91 (17p13.3), SOST (17q21.31_1) and DUSP3 (17q21.31_1) (Supplementary Table 13). Among them DKK1 (10q21.1) was the most significantly correlated with FN-BMD (P=1.3×10−5) and LS-BMD (P=3.2×10−4). Variants in all these BMD loci (with exception of 17p13.3) were also associated with fractures.

The SNP-eQTL analyses were performed across diverse tissues examining the correlation between marker alleles and transcript levels at the associated BMD loci. Fourteen of the BMD-associated SNPs correlated with the expression of one or more of the nearby genes with P < 5×10−5 and were either the strongest cis-variants, or good surrogates thereof, for those genes (Supplementary Tables 14 and 15). The most significant BMD-SNP eQTL was observed for rs10835187[T] with reduced expression of the LIN7C gene at the 11p14.1 locus (P = 2.8×10−39 in adipose tissue). Of particular interest were BMD-SNP cis-variants at three loci that were also associated with fracture including: 1p36.12, 4q22.1 and 17q21.31. At 1p36.12, rs6426749[G] correlated with reducedWNT4 expression in fibroblast, osteoblast and adipose tissue; at 4q22.1 rs6532023[G] correlated with reduced SPP1 (osteopontin) expression in adipose tissue and at 17q21.31 rs227584[A] correlated with increased C17orf65 expression in monocytes, adipose tissue, whole blood and lymphoblasts.

GRAIL analysis

We applied the GRAIL text-mining algorithm[19] to investigate connections between genes in the 55 autosomal BMD loci. This analysis revealed significant (GRAIL SNP P<0.01) connections between genes in 18 of the 55 input loci (Fig. 4 and Supplementary Table 16). The strongest connections were seen for members of three key biologic pathways: RANK-RANKL-OPG pathway (TNFRSF11A-TNFSF11-TNFRSF11B); mesenchymal stem cell differentiation (RUNX2, SP7, SOX9); and Wnt-signalling (LRP5, CTNNB1, SFRP4, WNT3, WNT4, WNT5B, WNT16, AXIN1) with the ten most frequently connecting terms being: ‘bone’, ‘catenin’, ‘signaling’, ‘differentiation’, ’rank’, ‘osteoblast’, ‘diacylglycerol’, ‘kappab’, ‘development’, and ‘osteoclast’. To assess the significance of this “biological” gene connection enrichment we applied GRAIL to 2000 random matched sets of 55 SNPs (See Supplementary note for details) and we did not observe any set with 15 or more loci with significant enriched connectivity (Supplementary Fig. 9) providing strong statistical evidence of the significant clustering of our BMD loci (P < 0.0005).

Figure 4

Graphic representation of GRAIL connections between SNPs and corresponding genes for the 18 SNPs as determined with GRAIL P<0.01. The top 10 keywords linking the genes were: ‘bone’, ‘catenin’, ‘signaling’, ‘differentiation’, ’rank’, ‘osteoblast’, ‘diacylglycerol’, ‘kappab’, ‘development’, and ‘osteoclast’. Thicker redder lines imply stronger literature-based connectivity. Blue and black boxes depict loci boundaries represented per top-associated marker (outer circle) and per gene in the region (inner circle).

Discussion

In this, the largest GWAS for osteoporosis traits to date, we identified 32 novel genomic loci bringing to 56 the number of loci robustly associated with BMD variation. Furthermore, we report for the first time that six of these BMD loci are associated with low-trauma fractures at P<5×10−8. As for other complex traits, our results indicate hundreds of variants with small effects may be contributing to the genetic architecture of BMD and fracture risk.[20] Our hypothesis-free assessment of common variants of the genome provides novel insights into biology, implicating several factors clustering in bone-active pathways.

Our results highlight the highly polygenic nature underlying BMD variation and the critical role of several biological pathways influencing osteoporosis and fracture susceptibility (Supplementary Fig. 10). On top of the Wnt factors known to be associated with BMD (CTNNB1, SOST, LRP4, LRP5, WLS, WNT4, MEF2C) several of the newly discovered loci also implicate additional Wnt signalling factors (including WNT5B, WNT16, DKK1, PTHLH, SFRP4 and AXIN1). Another clearly delineated pathway is that involved in mesenchymal cell differentiation, including the newly identified RUNX2, SOX4 and SOX9 BMD loci along with the previously known SP7. Another bone-relevant pathway includes that of “Endochondral Ossification” which involves essential processes during the fetal development of the mammalian skeleton and which implicated several of our identified BMD loci including: SPP1, MEF2C, RUNX2, SOX6, PTHLH, SP7 and SOX9. In addition, the biological relevance of our associations is accentuated by the identification of genes underlying rare monogenetic forms of osteoporosis and/or high bone mass such as SOST, CLCN7, LRP5 [21-23] (Supplementary Table 17) which also contain common variants involved in normal BMD variation at the population level.[11,14,16] This is supportive of a genetic architecture where both common and rare genetic variation may reside in the same locus.[24] Other genes have not been reported to be associated with monogenic forms for osteoporosis but have clear involvement in bone development in animal models. For example, SNPs in the 16q12.1 BMD locus map near CYLD. Human mutations in this gene have been described to cause familial cylindromatosis a condition without phenotypic skeletal manifestations. However, it has been shown that Cyld knock-out mice have significant bone loss leading to a severe osteoporosis phenotype[25] and also that CYLD regulates osteoclastogenesis.[26] Moreover, evidence from the GWAS and eQTL analyses also suggests some loci contain more than one common variant with independent effects on BMD and fracture risk. On the other hand, when no correlation is observed between gene expression and a particular SNP, it is difficult to draw conclusions. A correlation might be missed if the expression of the transcript was not measured in a relevant tissue or if the expression of a particular splice-variant was not measured.[27]

BMD and fracture genetic effects correlate to some extent, but some important fracture risk variants may have minimal impact on BMD and vice versa. This is the case for the 18p11.21 signal (Fig. 2B) mapping to a gene coding for a protein of unknown function, which despite a modest effect on BMD (0.02% variance explained) displayed the most significant association with fracture risk (OR=1.08, 95%CI[1.06–1.10], P=8.8×10−13). This is in contrast to variants with known stronger effects on BMD which were not significantly associated with fracture risk. For example, variants in the RANK-RANKL-OPG pathway, known to play a critical role in osteoclastogenesis, had clear associations with BMD but not fracture risk (Fig. 2A). Even though loci discovery was based on the BMD phenotype, these findings reflect the heterogeneous and complex nature of the mechanistic pathways leading to fracture. Therefore, given our study design, we cannot rule out the possibility that yet unidentified genetic loci are influencing risk of fracture independently of BMD. Future well-powered GWAS meta-analyses on fracture risk will address this question while corroborating the associations with fracture that we report for some of the BMD loci (particularly those not associated with fracture at P<5×10−8).

Our study also provides indication that there is sex- and site-specificity underlying BMD variation. One of the GWAS signals (Xp22.31) was only significant in the sex-stratified analysis in men and displayed significant sex heterogeneity (Phet=1.62×10−8). This is expected considering the sexual dimorphism of bone.[28,29] In fact, in a recent GWAS, the rs5934507 SNP mapping to Xp22.31, which is associated with BMD in the current study, has been previously associated with male serum testosterone levels.[30] Thus, it is likely that rs5934507 affects serum testosterone, which in turn regulates BMD. In line with the different types of bone composition at the different skeletal sites (predominantly trabecular at the lumbar spine while predominantly cortical at the femoral neck) we observed some indication of site specificity in 10 of the 56 BMD loci, suggesting differential genetic influences on BMD determination across skeletal sites. As has been previously shown[31], we did not find in our results major differences in effect sizes between individuals of European and East Asian ancestry (Supplementary Fig. 7). However, this may be due to reduced power given the smaller number of individuals of East Asian ancestry. We tested a genetic risk score to identify individuals at risk of osteoporosis and fracture and showed that cumulatively, the identified variants generate a gradient of risk. These gradients reach ORs of 1.56 for osteoporosis and 1.60 for fractures when comparing participants with the highest risk scores with those reflecting the mean score. Yet, at present there is limited clinical utility in using this score as evidenced by the non-significant contribution to case discrimination after considering clinical risk factors with strong effects on osteoporosis and fracture risk (like age and weight). This is not unexpected given the small fraction of genetic risk for either BMD or fracture that has been identified thus far.

Our study has limitations. The identified SNPs are probably not the causal variants; it is more likely that these markers are in LD with the underlying causal variants. Additional analyses on potential functional SNPs identified in this study will be required to determine if they are causal to these relationships with BMD. Moreover, the causal genes underlying the GWAS signals may be different from the candidate genes we describe, considering that our understanding of their role in bone biology is limited. Further exploration of these loci with more detailed sequencing, gene expression, and translational studies will be required. Such studies can also disentangle the diverse types of complex relationships we currently cannot distinguish in the BMD loci with secondary signals, i.e., if these are the result of true allelic heterogeneity or if they are driven by a second gene in the same region.[32] Similarly, despite our large sample size, power limitations still play a role for detecting additional associations with smaller effect sizes and/or arising from rarer variants. Finally, given the different levels of data availability and the difficulties for standardization across studies, we did not evaluate the effect of additional risk factors for osteoporosis, such as menopausal status and smoking, which can influence the genetic associations with BMD. Nonetheless, despite these limitations we have identified many novel and previously unsuspected associations with BMD variation and fracture risk.

Finally, the relatively weak effects of the variants discovered by GWAS do not undermine the biological relevance of the genes identified, as exemplified by the identification of genetic signals at the location of genes coding for proteins currently targeted by novel osteoporosis treatments (Supplementary Fig. 10). The novel genes identified in our study may represent new candidates to target for osteoporosis drug discovery. Most established treatments for osteoporosis currently focus on curtailing bone resorption (eg. bisphosphonates, RANKL inhibitors) while only few anabolic treatments are currently approved for the treatment of osteoporosis (i.e. recombinant truncated or altered PTH). Other anabolic compounds under Phase II development include PTHrP fragments and Wnt-signaling enhancers such as anti-Sclerostin antibodies.[33] Several of the variants robustly associated with BMD map in or close to genes of proteins involved in these pharmacologic pathways, namely osteoprotegerin (TNFRSF11B), RANK (TNFRSF11A), RANKL (TNFSF11), PTHrP (PTHLH), Low-density lipoprotein receptor-related protein 5 (LRP5), Sclerostin (SOST), and Dickkopf-1 (DKK1).

In conclusion, these findings highlight the highly polygenic and complex nature underlying BMD variation, shedding light on the pathophysiological mechanisms underlying fracture susceptibility and harbouring potential for the future identification of drug targets for the treatment of osteoporosis.

ONLINE METHODS

Study design

This study is part of the GEnetic Factors for OSteoporosis consortium (GEFOS), a coalition of teams of investigators dedicated to identify the genetic determinants of osteoporosis. The discovery samples comprised 17 GWA studies (n=32,961) from populations across North America, Europe, East Asia and Australia, with a variety of epidemiological designs (Supplementary Table 18A) and patient characteristics (Supplementary Table 18B); a subset of which had fracture information (Supplementary Table 18C). Subjects from 34 additional studies with BMD data (n=50,933) were used for replication while association with fracture was tested across 50 studies with fracture information, most of them also used for the BMD analysis (n=31,016 cases and 102,444 controls) (Figure 1 and Supplementary Tables 19A-C and 20A-C). All studies were approved by their institutional ethics review committees and all participants provided written informed consent.

BMD measurements and fracture definition

BMD of the lumbar spine (LS-BMD) and femoral neck (FN-BMD) was measured in all cohorts using dual-energy X-ray absorptiometry following standard manufacturer protocols (Supplementary Tables 18B, 19B and 20B). Three clinically-distinct fracture definitions were used: 1) Any type, consisting of low-trauma fractures at any skeletal site (except fingers, toes and skull) occurring after age 18 years assessed by X-ray, radiographic report, clinical record, clinical interview and/or questionnaire; 2) Validated non-vertebral, consisting of fractures occurring after age 50 years with diagnosis confirmed by hospital records and/or radiographs; and 3) Radiographic vertebral fractures, from lateral morphometry scored on X-rays. The first is most-inclusive, while the latter two are more stringent fracture definitions commonly used in randomized trials.[35,36] Controls were defined as individuals without a history of fracture using for each fracture type the same age limit categories of the cases.

Stage 1 Genome-wide association analysis

Genotyping and Imputation

GWAS genotyping was done by each study following standard manufacturer protocols followed by imputation to ~2.5 million SNPs from HapMap[37] Phase II release 22 using Genome Build 36. Quality control was performed independently for each study. To facilitate meta-analysis, each group performed genotype imputation with BIMBAM[38], IMPUTE[39], or MACH[40] software using genotypes from the HapMapPhase II release 22 (CEU or CHB/JPT as appropriate). HapMap release 21 was used as reference for SNPs residing on the X chromosome and IMPUTE software was used for imputation. Overall imputation quality scores for each SNP were obtained from IMPUTE (proper_info) and MACH (rsq_hat) statistics. Details on the genotyping platform used, genotype quality control procedures and software for imputation employed for each study are presented in the Supplementary Tables 18D and 19D.

Association analysis with BMD

each study performed genome-wide association analysis for FN-BMD and LS-BMD using sex-specific, age- weight- and principal components-adjusted standardized residuals analyzed under an additive (per allele) genetic model. Analyses of autosomal and chromosome X markers were done separately. Analysis of imputed genotype data accounted for uncertainty in each genotype prediction by using either the dosage information from MACH or the genotype probabilities from IMPUTE and BIM-BAM. Studies used MACH2QTL[40] directly or via GRIMP[41] (which uses genotype dosage value as a predictor in a linear regression framework), SNPTEST[39], Merlin[42], BIM-BAM or the linear mixed effects model of the Kinship and ProbABEL[43] (Supplementary Tables 18D and 19D). For analysis of the X-chromosome either SNPTEST or R package was used in each participating study. We coded “effect allele homozygous genotype” as “2” and “other allele homozygous genotype” as “0” in the genotyped SNPs in men on the X chromosome. The imputed genotypes were coded as continuous variables from 0 to 2 to take into account imputation uncertainty. The genomic control method[44] was used to correct the standard error (SE) by the square root of the genomic inflation factor (lambda): SEcorrected = SE × √lambda.

Meta-analysis of the GWA studies

before performing meta-analysis on the genome-wide association data, SNPs with poor imputation quality scores (rsq_hat < 0.3 in MACH, proper_info < 0.4 in IMPUTE or the ratio of observed to expected dosage variance < 0.3 in BIMBAM) and markers with a minor allele frequency < 1% were excluded for each study. All individual GWAS were genomic control corrected before meta-analysis.[44] Individual study-specific genomic control values ranged from 0.98 to 1.08. (Supplementary Table 18D). A total of 2,483,766 autosomal SNPs were meta-analyzed across: 17, 16 and 13 studies for FN-BMD (Pooled, women-only, and men-only analyses, respectively) and 16, 13 and 12 studies for LS-BMD (Pooled, women-only, and men-only analyses, respectively). A total of 76,253 X-linked SNPs were meta-analyzed across: 14, 13 and 10 studies for LS and FN-BMD (Pooled, women-only, and men-only analyses, respectively). In our discovery analysis, we chose to implement a fixed effects models approach as it is generally preferable for the purposes of initial discovery, where the aim is to screen and identify as many of the true variants as possible.[45,46] SNPs present in less than three studies were removed from the meta-analysis yielding ~ 2.2 million SNPs in the final results. Genomic inflation factors (λ) were 1.11, 1.09, 1.06 for FN-BMD BMD (Pooled, women-only, and men-only analyses, respectively) and 1.13, 1.09, 1.06 for LS-BMD (Pooled, women-only, and men-only analyses, respectively). A second GC correction was applied to the overall meta-analysis results, although such second correction is considered overly conservative.[47] Significance for BMD association was set at P<5×10−8 while a Bonferroni correction was used for the association with fracture.[48]

Selection of SNPs for follow up

we took forward the most significant 96 SNPs for replication. Based on power estimations, after adding 30,000 samples in stage 2 these variants had a priori Power >=85% to reach P= 5 × 10−8 in the meta-analysis. Loci were considered independent when separated by at least 1 Mb down and upstream of the top GWAS signal. The 96 variants included the 82 index SNPs representing each of the 82 loci reaching P<5×10−6 in Stage 1, 9 SNPs that lie within the same 2Mb windows as the 82 but which were independent from the main signal (secondary signals), and the top-five most associated SNPs of the X-chromosome (with P <5×10−5).

Association analyses with fracture risk

effect estimates (odds ratio) for association of allele dosage of the top hits with fracture risk were obtained from logistic regression models adjusted for age, age[2] weight, sex, height and four principal components. The proportion of the fracture risk explained by FN-BMD was calculated from the regression coefficients as (βunadjusted - βBMDadjusted)/βunadjusted in a subset of replication samples for which both FN-BMD and complete fracture information was available.

Stage 2 follow up

Samples and genotyping

fracture association results were also obtained for the 82 most significant SNPs from 54,244 individuals of European ancestry from 7 GWAS (in-silico genotyping) that had not been included in the stage 1 analyses (Supplementary Tables 19A, 19B and 19C). Subjects from 34 studies of the GENOMOS consortium with BMD and/or fracture information were studied for replication (Supplementary Tables 3A, 3B and 3C). De-novo replication genotyping was done in the UK (Kbiosciences), Iceland (deCODE Genetics), Australia (University of Queensland Diamantina Institute) and the USA (WHI GeCHIP) using KASPar, Centaurus, OpenArray and iSelect assays respectively (Supplementary Note). Minimum genotyping quality control criteria were defined as: Sample call rate > 80%, SNP call rate > 90%, HWE P > 1×10−4, MAF > 1%.

Association analyses and meta-analysis

We tested the association between the 96 SNPs and BMD and fracture risk in each in-silico and de-novo “Stage 2” study separately as described for the “Stage 1” studies. We subsequently meta-analyzed effects and standard errors from the “Stage 2” studies, followed by a meta-analysis of the summary statistics of both “Stage 1” and “Stage 2” using the inverse-variance method in METAL. At this replication stage, where more than 30 studies were synthesized, we chose to first assess the underlying heterogeneity considering both the Cochran’s Q statistic and the I2 metric. If the heterogeneity was not significant fixed effects models were applied. If the Cochrane Q P-value<0.0005 and the I2 was > 50% we used the more conservative random effects models.

Additional analyses

Further analyses were performed for the SNPs carried forward for replication. Each of these analyses is described in detail in the “Supplementary Note”. In brief, we performed: 1) a conditional genome-wide association analysis to examine whether any of the 82 BMD loci harbored additional independent signals; 2) tested gene-by-gene pair-wise interactions between these BMD loci; 3) assessed within the independent setting of the PERF study (for details on study design see Supplementary Tables 20A, 20B & 20C) the predictive ability derived from the cumulative effect of the 63 genome-wide significant autosomal BMD SNPs in relation to BMD levels and osteoporosis risk; and that of the 16 BMD SNPs also associated with fracture risk in relation to fracture risk; 4) identified SNPs having r2 ≥ 0.80 with the lead SNP that were potentially functional (nonsense, nonconservative non-synonymous, synonymous, exonic splicing, transcription factor binding sites, etc) using regional imputation with the 1000 Genomes data (June 2010 release); 5) tested the relationship between gene expression profiles from a) trans-iliacal bone biopsies and BMD in 84 unrelated postmenopausal women[49] and b) also examined cis- associations between each of the 55 significant BMD SNPs and expression of nearby genes in different tissues including lymphoblastoid cell lines[50-52], primary human fibroblasts and osteoblasts[53], adipose tissue[54], whole blood[54] and circulating monocytes[55]; and finally 6) evaluated the connectivity and relationships between identified loci using the literature-based annotation with Gene Relationships across Implicated Loci (GRAIL[19]) statistical strategy.