Genetic Sharing with Cardiovascular Disease Risk Factors and Diabetes Reveals Novel Bone Mineral Density Loci.

Bone Mineral Density (BMD) is a highly heritable trait, but genome-wide association studies have identified few genetic risk factors. Epidemiological studies suggest associations between BMD and several traits and diseases, but the nature of the suggestive comorbidity is still unknown. We used a novel genetic pleiotropy-informed conditional False Discovery Rate (FDR) method to identify single nucleotide polymorphisms (SNPs) associated with BMD by leveraging cardiovascular disease (CVD) associated disorders and metabolic traits. By conditioning on SNPs associated with the CVD-related phenotypes, type 1 diabetes, type 2 diabetes, systolic blood pressure, diastolic blood pressure, high density lipoprotein, low density lipoprotein, triglycerides and waist hip ratio, we identified 65 novel independent BMD loci (26 with femoral neck BMD and 47 with lumbar spine BMD) at conditional FDR < 0.01. Many of the loci were confirmed in genetic expression studies. Genes validated at the mRNA levels were characteristic for the osteoblast/osteocyte lineage, Wnt signaling pathway and bone metabolism. The results provide new insight into genetic mechanisms of variability in BMD, and a better understanding of the genetic underpinnings of clinical comorbidity.

Introduction

Low bone mineral density (BMD) is an important human phenotype predisposing for bone fractures [1]. Primary and secondary osteoporosis, (defined as BMD less than 2.5 SD of young controls) occur frequently in all populations and lead to high risk for fractures and lasting functional impairment, resulting in long term personal suffering and high social costs [2]. Several lines of evidence suggest an overlap between BMD/osteoporosis and several traits related to metabolism and cardiovascular disease (CVD): -presence of osteoporosis is associated with a ~4-fold increase in risk for an acute cardiovascular event [3].—BMD loss is associated with increased mortality from coronary heart disease and pulmonary diseases [4]- an inverse relationship is found between high-density lipoprotein (HDL) cholesterol and BMD [5-9]. The relationship between low-density lipoprotein (LDL) cholesterol and BMD appears to be less profound, but a positive association has been found in some studies [5,10]. While not all studies have identified a relationship between Triglycerides (TG) and BMD, a few larger studies have shown an inverse relationship [7,8,10]. Furthermore, statins are widely used as cholesterol-lowering drugs, and a recent meta-analysis indicates that statins may help improve and maintain BMD at the lumbar spine, hip and femoral neck, especially in Caucasians and Asians [11].

Blood pressure and anthropometric measures have also been found to be associated with BMD in epidemiological studies. Lee et al. [12]. found that both high systolic blood pressure (SBP) and high diastolic blood pressure (DBP) were associated with low femoral BMD, but not with lumbar BMD in a total study sample consisting of 8439 men and postmenopausal women aged 50 years and older. A study of 586 postmenopausal Turkish women also showed a significant correlation between SBP and femur BMD [13]. It should be noted that several studies also failed to find a link between blood pressure and osteoporosis, e.g. [14].

There is also clinical and epidemiological evidence for association between BMD and metabolic traits. As reviewed [15-17], it is well documented that Type 1 Diabetes (T1D) and Type 2 Diabetes (T2D) increase risk of fracture. Also, it is well established that a major part of the increased fracture risk in T1D is caused by reduced BMD, due to defects in osteoblast differentiation and activity as well as contributing factors including accumulation of advanced glycation end products (AGEs)[18]. Thus, it is plausible that the microenvironment in which B cells develop, the bone marrow including osteoblasts, is influenced by genetic factors that affect both an autoimmune disease like T1D and osteoporosis.

The relationship between T2D and BMD or fracture is more complicated, since the effect on bone microstructure appears to be more important. However, Sayers et al. [19] found an inverse association between insulin and both periosteal circumference and cortical BMD in adolescents after adjusting for all body composition variables, indicating that insulin levels and diabetes have effects on bone metabolism. In adults T2D has been associated with high BMD [16,17] and Billings et al. [20] identified Integrin, Alpha 1 (ITGA1) as a new locus candidate, capable of influencing both fasting glucose and BMD, thus pointing to a possible explanation for the epidemiological observations linking T2D diabetes and BMD/osteoporosis. The previous concept, that obesity is protective for osteoporosis is weakened since several studies have shown a negative correlation between WHR and BMD [21-23]. Many of the previous studies did not take into consideration that DXA measurements are falsely elevated by increased body fat and that the associated increase in bone marrow adiposity occurs at the expense of bone [23].

The co-morbidity between BMD and CVD risk factors or metabolic traits have been postulated to be, at least partly, caused by overlapping genes (pleiotropy) [24]. GWAS have identified several genes and single nucleotide polymorphisms (SNPs) associated with BMD [25], and CVD risk factors or metabolic traits, including HDL [26], LDL [26], TG [26], T1D [27], T2D [28], SBP [29], DBP [29] and WHR [29]. Despite the strong heritable component of BMD, the genes identified in GWAS so far explain only a small proportion of the variance (‘missing heritability‘) [25]. Due to the polygenic architecture of BMD, a large number of SNPs have associations too weak to be identified in the currently available cohorts. Thus, pleiotropic enrichment together with cost-effective analytical methods may identify a larger proportion of SNPs associated with BMD.

Standard methods to assess genetic pleiotropy have not taken full advantage of the existing GWAS data and the majority of these studies have focused on the subset of SNPs exceeding a Bonferroni-corrected threshold of significance for each trait or disorder [30,31]. However, this Bonferroni–based approach cannot detect SNPs that only reach genome-wide significance in the combined analysis but do not meet significance cutoffs in the individual phenotype. In the current study, we applied a recently developed genetic pleiotropy-informed approach for GWAS to leverage the power of multiple large GWAS of CVD risk factors blood lipids (HDL, LDL, TG), metabolic disorders (T1D, T2D), blood pressure (SBP, DBP), and waist-hip ratio (WHR) to identify susceptibility SNPs, and capture more of the polygenic effects in BMD [32-34]. This novel genetic epidemiological approach is able to take advantage of polygenic pleiotropy among several types of diseases to identify genetic variants with smaller effect sizes, and thus elucidate the mechanism of variability in BMD. We used summary statistics (p-values and allele frequencies) from the analysis data (up to 32,961 individuals) in the primary study of BMD [25] for both femoral neck (FN) and lumbar spine (LS) BMD phenotypes.

Materials and Methods

Participant Samples and Statistical Strategy

The study was approved by the Norwegian Regional Ethical Committee (REK no: 2010/2539) and conducted according to the Declaration of Helsinki (2000). Written informed consent was given by participants for their clinical records to be used in this study. We obtained complete stage 1 GWAS results in the form of summary statistics p-values from public access websites or through collaboration with investigators (T1D cases and controls from The Type 1 Diabetes Genetics Consortium, BMD cases and controls from the GEFOS Consortium). There was some overlap among several of the participants in the anthropometric GWAS and the BMD GWAS sample (for further details, see S1 Table).

Statistical Analyses

Overall Approach

After applying genomic inflation control, we compute the stratified empirical cumulative distribution functions (cdfs) of the nominal p-values. Strata are determined by relative enrichment of pleiotropic SNPs in BMD as a function of increased nominal p-values in the different associated traits and disorders. Using this stratified methodology, we construct two-dimensional FDR “look-up” tables (S1 and S2 Figs), with FDR in BMD SNPs computed conditional on nominal associated phenotypes p-values (conditional FDR). Using this table we identify loci that are significantly associated with BMD at a conditional FDR level of 0.01. All p-values were corrected for inflation using the genomic control procedure [35], and for overlap in samples [36] as previously described [37]. Finally, the SNP gene associations were validated using information from global transcriptional mapping of bone biopsies from postmenopausal women [38,39].

Genomic Control

The empirical null distribution in GWAS is affected by global variance inflation due to population stratification and cryptic relatedness and deflation due to over-correction of test statistics for polygenic traits by standard genomic control methods. We used the same formulism as in Schork et al. [35]. The genomic inflation factor λGC for each phenotype were estimated based on intergenic SNPs as the median z-score squared divided by the expected median of a chi-square distribution with one degree of freedom and divided all test statistics by λGC. We have previously reported that intergenic SNPs, as defined in our annotation protocol (Schork et al, 2013) are deplete of association with >30 complex traits/diseases, and it seems that this is a generic feature for SNPs in this category. Furthermore, intergenic SNPs do not show skewed distribution towards small minor allele frequency (MAF) based on the 1000 Genomes Project (1KGP) [32,33,37].

Conditional Q-Q Plots for Pleiotropic Enrichment

To assess pleiotropic enrichment, we used Q-Q plot conditional by ‘pleiotropic’ effects as described in detail earlier (Fig 1) [33,34,37]. For a given associated phenotype, enrichment for pleiotropic signals is present if the degree of deflection from the expected null line is dependent on SNP associations with the second phenotype. Specifically, we computed the empirical cumulative distribution of nominal p-values for a given phenotype for all SNPs and for SNPs with significance levels below the indicated cut-offs for the other phenotype (–log10(p) ≥ 0,–log10(p) ≥ 1,–log10(p) ≥2,–log10(p) ≥3 corresponding to p < 1, p < 0.1, p < 0.01, p < 0.001, respectively). The nominal p-values (–log10(p)) are plotted on the y-axis, and the empirical quantiles (–log10(q), where q = 1-cdf(p)) are plotted on the x-axis. To assess for polygenic effects below the standard GWAS significance threshold, we focused the conditional Q-Q plots on SNPs with nominal–log10(p) < 7.3 (corresponding to p > 5x10-8).

Fig 1

Genetic enrichment.

Conditional Q-Q plot of nominal versus empirical -log10 p-values (corrected for inflation) in bone mineral density (BMD, femoral neck) below the standard GWAS threshold of p < 5x10-8 as a function of significance of association with CVD risk factors, including systolic blood pressure (SBP), diastolic blood pressure (DBP), high density lipoproteins (HDL) and triglycerides (TG) at the level of -log10(p) ≥ 0 (all SNPs),–log10(p) ≥ 1,–log10(p) ≥ 2,–log10(p) ≥ 3 corresponding to p ≤ 1, p ≤ 0.1, p ≤ 0.01, p ≤ 0.001, respectively. Dotted lines indicate the null-hypothesis.

Conditional Statistics–Test of Association with BMD

To improve detection of SNPs associated with BMD, we used a conditional False discovery rate (FDR) approach, leveraging pleiotropic phenotypes [32-34,37]. Specifically, the conditional FDR of a trait (e.g. BMD) for a SNP with p-value P 1 on a second pleiotropic trait with p-value P 2, is computed as the posterior probability that the SNP is null for the first trait given that the p-values for both phenotypes are as small as or smaller than the observed p-values, FDR(P 1│P 2) = π0(P 2)P 1/F(P 1│P 2), where F(P 1│P 2)is the conditional cdf and π0(P 2)the conditional proportional of null SNPs for the first phenotype given that p-value for the second phenotype are P 2 or smaller. The values of FDR(P 1│P 2) were conservatively estimated by setting π0(P 2) equal one and replacing F(P 1│P 2) by empirical conditional cdf. The conditional FDR values for BMD on second pleiotropic traits (denoted by FDR , where the dot denotes a second phenotype) were assigned, based on the combination of p-value for the SNP correlated to BMD and the associated trait, by interpolation into a 2-D look-up table (S1 and S2 Figs). All SNPs with FDR < 0.01 (-log10(FDR) > 2) in BMD given the different associated phenotypes were identified. A significance threshold of FDR < 0.01 corresponds to 1 false positive per 100 reported associations.

Annotation of Novel Loci

Based on 1KGP linkage disequilibrium (LD) structure, significant SNPs identified by conditional FDR were clustered into LD blocks at the LD-r2 > 0.2 level. This threshold was chosen since it has been used in a large number of reported GWAS, thus making our result comparable to previous studies, e.g.[25,39,40]. The blocks were numbered as loci # in Table 1 and S2, S3 and S4 Tables and any one block may contain more than one SNPs. Genes close to each SNPs were obtained from NCBI gene database. Blocks that do not contain SNPs or close-by genes to SNPs from primary study were deemed as novel loci in current study (Table 1 and S3 Table). And, loci that contain SNPs or genes from primary study were considered as replication of primary findings (S2 and S4 Tables for FN and LS BMD, respectively). The same procedure was applied to both FN BMD and LS BMD phenotypes. To identify non-overlapping loci between FN BMD and LS BMD, the SNP rs-numbers and gene symbols for these two phenotypes were compared. Loci containing SNPs with same rs-number or same genes were considered overlapping.

Table 1

Novel femoral neck BMD associated genes at conditional FDR <0.01.

									Expressed QTL (Age and BMI adj.)
Loci #	SNP	Gene symbol	Map Loc.	BMD pvalue	BMD FDR	Min cond FDR	Waldstats	Drivingphenotype	Affymetrix ID	r
1	rs10779702	RERE	1p36.23	7,78E-08	3,06E-04	1,60E-04	-5.26	HDL	200940_s_at	-0.23
6	rs12137389	TESK2	1p32	1,88E-06	4,15E-03	4,01E-03	4.67	HDL	206758_at	-0.11
9	rs11809524	COL11A1	1p21	8,21E-07	2,03E-03	1,34E-03	-4.83	SBP	37892_at	0.25
11	rs9309664	PPP1CB	2p23	7,55E-06	1,20E-02	8,22E-03	4.39	HDL	228222_at	-0.30
15	rs11675051	TMEM194B	2q32.2	1,46E-06	3,47E-03	1,56E-03	-4.72	SBP	238014_at	0.09
15	rs13005335	NAB1	2q32.3-q33	1,54E-06	3,47E-03	1,56E-03	-4.71	SBP	209272_at	0.05
16	rs12995369	CDK15	2q33.2	1,07E-07	3,69E-04	2,80E-04	-5.2	SBP	1552559_a_at	0.16
17	rs7594560	METTL21A	2q33.3	3,42E-06	5,91E-03	3,74E-03	4.55	HDL	1553743_at	-0.11
23	rs4957742	RAB9BP1	5q21.2	2,98E-06	5,91E-03	6,27E-03	-4.58	DBP	NA	NA
27	rs6583337	FAM20C	7p22.3	3,30E-06	5,91E-03	3,38E-03	4.56	LDL	229438_at	0.18
29	rs2282930	GRB10	7p12.2	5,20E-06	8,40E-03	7,20E-03	4.46	TG	210999_s_at	-0.35
32	rs10953178	C7orf76	7q21.3	3,75E-11	6,36E-07	3,53E-07	-6.48	HDL	NA	NA
32	rs10464592	SHFM1	7q21.3	4,28E-10	2,35E-06	4,07E-06	6.11	SBP	202276_at	-0.05
35	rs1670346	PTPRN2/MIR595	7q36	1,73E-06	3,47E-03	1,80E-03	-4.68	SBP	203030_s_at	0.16
37	rs980299	EYA1	8q13.3	1,18E-07	4,45E-04	3,39E-04	5.19	HDL	214608_s_at	-0.03
38	rs13272568	PKIA	8q21.11	1,29E-06	2,90E-03	2,52E-03	4.74	SBP	204612_at	-0.35
40	rs665556	KLF4	9q31	6,68E-06	1,00E-02	5,84E-03	4.41	DBP	220266_s_at	-0.34
49	rs600231	MALAT1	11q13.1	7,75E-06	1,20E-02	7,60E-03	-4.38	SBP	231735_s_at	0.29
51	rs258415	KLHL42	12p11.22	3,55E-08	1,69E-04	1,43E-04	-5.4	SBP	NA	NA
53	rs11614913	MIR196A2	12q13.13	4,20E-08	1,69E-04	1,25E-04	5.37	SBP	NA	NA
54	rs10746070	RIC8B	12q23.3	2,14E-06	4,15E-03	3,13E-03	-4.64	HDL	229637_at	0.04
58	rs7175531	CYP19A1	15q21	2,30E-06	4,96E-03	4,82E-03	-4.63	HDL	240705_at	0.25
58	rs10851498	MIR4713	15q21	2,73E-06	4,96E-03	4,45E-03	-4.59	TG	NA	NA
61	rs3198697	PDXDC1	16p13.11	1,01E-05	1,44E-02	5,00E-03	4.32	HDL	212053_at	0.08
67	rs199529	NSF	17q21	2,39E-06	4,96E-03	2,55E-03	4.62	SBP	202395_at	-0.13
71	rs8090312	NFATC1	18q23	4,54E-06	8,40E-03	6,40E-03	-4.49	T1D	211105_s_at	-0.15
74	rs756632	RTDR1/GNAZ	22q11.2	3,39E-06	5,91E-03	4,75E-03	-4.55	HDL	220105_at	-0.01
74	rs4820539	RAB36	22q11.2	3,06E-06	5,91E-03	7,44E-03	4.57	HDL	211471_s_at	0.14

Independent complex or single gene loci (LD-r2 < 0.2) with SNP(s) with a conditional FDR (condFDR) < 0.01 in bone mineral density (BMD, Femoral neck) given the association in other phenotype. We defined the most significant BMD SNP in each LD block based on the minimum condFDR (min condFDR) for each phenotype. The most significant SNPs in each gene of the LD block are listed and the second phenotype which provides the minimal FDR signal (Driving phenotype). All loci with SNPs with condFDR < 0.01 were used to define the number of the loci. The following abbreviations were used: Type 1 diabetes (T1D), type 2 diabetes (T2D), systolic blood pressure (SBP), diastolic blood pressure (DBP), low-density lipoproteins (LDL) cholesterol and high-density lipoproteins (HDL) cholesterol, chromosome location (Map Loc.). BMD FDR values < 0.01 are in bold. Bold r values represent nominally significant (p<0.05) Pearson correlations. Gene titles and ontology function terms are presented in S5 Table. Wald stats: z-score transformed from p values NA: not applicable (undetected)

Conditional FDR Manhattan Plots

To illustrate the localization of the genetic markers associated with BMD given the CVD risk factor effect, we used a ‘Conditional FDR Manhattan plot’, plotting all SNPs within an LD block in relation to their chromosomal location. As illustrated in Fig 2 and S3 Fig, the large points represent the SNPs with FDR < 0.01, whereas the small points represent the non-significant SNPs. All SNPs without ‘pruning’ (removing all SNPs with LD-r2 > 0.2 based on 1KGP LD structure) are shown. The strongest signal in each LD block is marked by larger points with black edges. This was identified by ranking all SNPs in increasing order, based on the conditional FDR value for BMD, and then removing SNPs in LD-r2 > 0.2 with any higher ranked SNP. Thus, the selected locus was the most significantly associated with BMD in each LD block (Fig 2 and S3 Fig).

Fig 2

‘Conditional FDR Manhattan plot’ of conditional–log10 (FDR) values for bone mineral density (BMD, femoral neck) alone (small black dots) and BMD given the associated phenotypes type 1 diabetes (T1D; BMD|T1D), type 2 diabetes (T2D; BMD|T2D), waist hip ratio (WHR, BMD|WHR), systolic blood pressure (SBP, BMD|SBP), diastolic blood pressure (DBP, BMD|DBP), high density lipoproteins (HDL, BMD|HDL), low density lipoproteins (LDL, BMD|LDL) and triglycerides (TG, BMD|TG).

SNPs with conditional–log10 FDR > 2 (i.e. FDR < 0.01) are shown with large points. A black line around the large points indicates the most significant SNP in each LD block and this SNP was annotated with the closest gene which is listed above the symbols in each locus, except for the HLA region on chromosome 6. Gene symbols were obtained from NCBI gene databases and colored in line with the second phenotype which gives the minimal conditional FDR value. Details for the novel loci with–log10 FDR > 2 (i.e. FDR < 0.01) are shown in Table 1 and S1 Table. Genes previously reported by other studies were marked by stars (*).

Validation by Expression Genetics

We looked for expressional association between the SNP associated genes and BMD in bone biopsies from postmenopausal women (n = 84) [38,39]. The Iliac biopsies were analyzed with Affymetrix microchips and log2 transformed signal values were correlated to BMD levels (Table 1, S2 Table). The primary data have been submitted to the European Bioinformatics Institute (EMBL-EBI; ID: E-MEXP-1618).

Results

Pleiotropic Enrichment-Polygenic Overlap

Conditional Q-Q plots for FN BMD conditioned on nominal p-values of association with T1D, T2D, SBP, DBP, HDL, LDL, TG and WHR showed enrichment across different levels of significance (Fig 1 and S5 Fig). Similar plots for LS BMD are shown in S4 Fig. The earlier departure from the null line (leftward shift) suggests a greater proportion of true associations for a given nominal FN BMD p-value (See S1 File for detailed explanation). Successive leftward shifts for decreasing nominal p-values of a second phenotype indicate that the proportion of non-null effects varies considerably across different levels of association with the comorbidity trait or disease.

Loci Associated with BMD

To identify SNPs associated with FN BMD, we constructed a “conditional FDR” Manhattan plot showing the FDR conditional on each of the risk factors (Fig 2). We identified significant loci associated with FN BMD leveraging the reduced FDR obtained by the associated phenotype. To estimate the number of independent loci, we pruned the associated SNPs (removed SNP with LD-r2 > 0.2), and identified a total of 74 independent loci with a conditional FDR < 0.01 of which 15 were complex loci and 59 single gene loci (marked in Fig 2 by points with black edges). The 74 loci encompassed 84 different genes. Using the FDR method in FN BMD alone, 70 loci were identified (bold values in the “BMD-FDR” column, Table 1 and S2 Table). The remaining 4 loci would not have been identified in the current sample without using the conditional FDR method. Similarly, the 95 independent loci for LS BMD encompasses 107 different genes, and the extra number of loci identified with our conditional FDR compared with FDR method is 21 (bold value in the “BMD-FDR” column, S3 and S4 Tables, marked by points with black edges in S3 Fig). Since there are overlaps in loci between the two phenotypes, we identified a total of 122 independent loci for FN BMD and/or LS BMD, representing 155 different genes in all.

SNP Detection and Verification

The previous study of BMD related SNPs by Estrada et al. [25] identified a total of 56 loci associated with FN BMD and/or LS BMD (49 loci with FN BMD and 48 loci with LS BMD). This was based on two-stage analysis (consisting up to 83,894 and 77,508 individuals), whereas in the stage-1 sample analysis (consisting 32,961 and 31,800 individuals), 20 and 26 loci were associated with FN BMD and LS BMD, respectively. Our analysis re-identified all (20 FN and 26 LS) loci reported in the primary study stage-1 analysis by Estrada et al. [25]. Also in the cross stage (I and II) analyses, all but 5 loci for FN BMD and 8 loci for LS were successfully re-identified (S2 and S4 Tables).

The FDR method identified 26 novel loci associated with FN BMD and 47 novel loci associated with LS BMD, not reported in the previous BMD GWAS [25].

Gene Expression Analysis

Global gene expression profiling in iliac bone biopsies from 84 postmenopausal women [38] permitted us to calculate the correlation values between BMD and the mRNA levels of all genes associated with the identified loci, as shown in the rightmost columns of Table 1 (novel genes) and S2 Table (genes identified also by Estrada et al. [25]). We found a similar fraction of transcripts that were significantly correlated with FN BMD among the novel BMD associated genes (8 out of 26 reaching detection level), very similar to the Estrada study [25], 31% vs. 30%, respectively.

Functional Enrichment Analysis

The 155 genes encompassed by all loci at FDR < 0.01 for FN and LS BMD were analyzed with Ingenuity Pathway Analysis (IPA). The top-most significantly affected canonical pathway was “Role of Osteoblasts and Chondrocyte in Rheumatoid Arthritis” (p = 4.1x10-12), which includes Wnt signaling, and the function and interaction of many of the identified genes in bone related cells (Table 2).

Table 2

Top Canonical pathways and Top diseases and Bio Functions from Ingenuity Pathway Analysis

Canonical Pathways	Ratio(p-value)		Molecules
Role of Osteoblasts, Osteoclasts and Chondrocytes in Rheumatoid Arthritis	16/231 (4.10E-12)		SFRP4,RELA,LRP5,TNFSF11,SPP1,AXIN1,WNT2B,WNT16,SP7, TNFRSF11A, NFATC1, WNT4, BMP7, SOST, CTNNB1, TNFRSF11B
Role of Macrophages, Fibroblasts and Endothelial Cells in Rheumatoid Arthritis	12/329 (1.86E-6)		SFRP4, RELA, TNFSF11, LRP5, AXIN1, WNT2B, WNT16, WNT4, SOST,CTNNB1, NFATC1, TNFRSF11B
Wnt/Î2-catenin Signaling	9/174 (3.69E-6)		SFRP4, LRP5, SOX6, AXIN1, WNT2B, WNT16, WNT4, SOX9, CTNNB1
Basal Cell Carcinoma Signaling	6/75 (1.32E-5)		AXIN1, WNT2B, WNT16, WNT4, BMP7, CTNNB1
Role of NANOG in Mammalian Embryonic Stem Cell Pluripotency	7/117 (1.52E-5)		AXIN1, WNT2B, WNT16, WNT4, BMP7, CTNNB1, ZFP42
Human Embryonic Stem Cell Pluripotency	7/153 (5.06E-5)		AXIN1, SMAD3, WNT2B, WNT16, WNT4, BMP7, CTNNB1
Colorectal Cancer Metastasis Signaling	9/254 (5.87E-5)		RELA, LRP5, AXIN1, SMAD3, WNT2B, ADCY6, WNT16, WNT4,CTNNB1
Protein Kinase A Signaling	11/389 (9.02E-5)		DHH, RELA, PTPRD, SMAD3, ADCY6, PPP1CB, CTNNB1, EYA1,ANAPC1, NFATC1, AKAP11
Role of Wnt/GSK-3Î2 Signaling in the Pathogenesis of Influenza	5/82 (2.48E-4)		AXIN1, WNT2B, WNT16, WNT4, CTNNB1
Regulation of the Epithelial-Mesenchymal Transition Pathway	7/190 (3.13E-4)		RELA, AXIN1, SMAD3, WNT2B, WNT16, WNT4, JAG1
Categories in Top Diseases and Bio Functions	Diseases orFunctionsAnnotation	# Molecules(p-Value)	Molecules
Connective Tissue Development and Function, Embryonic Development, Organ Development, Organ Morphology, Organismal Development, Skeletal and Muscular System Development and Function, Tissue Development	abnormalmorphology ofbone	27 (1,86E-15)	ARHGAP1, BMP7, CYP19A1, ESR1, EYA1, FAM20C, GALNT3,HOXC4, HOXC5, HOXC6, IBSP, IDUA, LRP5, MEOX1, MEPE,NAB1, PKDCC, SALL1, SMAD3, SOST, SOX6, SOX9, SPP1,TNFRSF11A, TNFRSF11B, TNFSF11, ULK4
Organismal Development, Skeletal and Muscular System Development and Function	abnormalmorphology of limb	18 (2,48E-13)	BMP7, ESR1, EYA1, FAM20C, GALNT3, IBSP, IDUA, LRP4,LRP5, PKDCC, SALL1, SMAD3, SOST, SOX9, TNFRSF11A,TNFRSF11B, TNFSF11, WNT4
Skeletal and Muscular System Development and Function	abnormalmorphology ofskeleton	19 (1,16E-12)	ARHGAP1, BMP7, ESR1, EYA1, FAM20C, GALNT3, HOXC4,HOXC5, HOXC6, IBSP, IDUA, LRP5, MEOX1, PKDCC, SMAD3,SOST, SOX9, TNFRSF11B, TNFSF11
Connective Tissue Development and Function, Skeletal and Muscular System Development and Function	bone mineraldensity	15 (4,51E-12)	ARHGAP1, CYP19A1, ERCC1, ESR1, FAM20C, GALNT3, IBSP,LRP5, NAB1, SMAD3, SOST, SPP1, TNFRSF11A, TNFRSF11B,TNFSF11
Connective Tissue Development and Function, Embryonic Development, Organ Development, Organ Morphology, Organismal Development, Skeletal and Muscular System Development and Function, Tissue Development	morphology oflimb bone	13 (1,36E-11)	BMP7, ESR1, FAM20C, GALNT3, IBSP, IDUA, LRP5, PKDCC,SMAD3, SOST, SOX9, TNFRSF11B, TNFSF11
Connective Tissue Development and Function, Embryonic Development, Organ Development, Organismal Development, Skeletal and Muscular System Development and Function, Tissue Development	mineralization ofbone	13 (4,78E-11)	BMP7, ESR1, FAM20C, IBSP, LRP5, MEPE, PK DCC, SMAD3,SOST, SOX9, SPP1, TNFRSF11B, WNT4
Skeletal and Muscular System Development and Function	abnormalmorphology ofappendicularskeleton	13 (5,66E-11)	BMP7, ESR1, FAM20C, GALNT3, IBSP, IDUA, LRP5, PKDCC,SMAD3, SOST, SOX9, TNFRSF11B, TNFSF11
Cellular Development	differentiation ofconnective tissuecells	23 (1,94E-10)	AREG/AREGB, AXIN1, BMP7, CTNNB1, FAM20C, JAG1, KLF4,LGR4, LRP5, MEF2C, NFATC1, PKDCC, RELA, SFRP4, SMAD3,SOST, SOX9, SP7, SPP1, TNFRSF11A, TNFRSF11B, TNFSF11,WNT4
Organismal Injury and Abnormalities	calcinosis	9 (3,26E-10)	BMP7, CTNNB1, GALNT3, IBSP, LRP5, SOX9, SPP1,TNFRSF11B, TNFSF11
Cardiovascular Disease	degenerativemitral valvedisease	5 (3,48E-10)	CTNNB1, IBSP, LRP5, SOX9, SPP1

The genes associated with all identified loci (min Cond FDR < 0.01) were subjected to Ingenuity Pathway Analysis. The topmost significantly affected canonical pathways (upper panel) and Categories in Top Diseases and Bio Functions (lower panel) from the analysis are shown.

Out of all the loci at FDR < 0.01 (LS and FN BMD), 48 associated gene transcripts were significantly correlated to BMD in bone biopsies from postmenopausal women. This subset of genes was also analyzed by IPA, and a network of interacting genes including NFATC1, RELA, NFΚB and SMAD3 as central nuclear hubs were generated (Fig 3).

Fig 3

Network analysis IPA-generated network illustrating molecular interactions among the genes correlating inversely (green) or positively (red) to FN or LS BMD.

All genes associated with FN or LS SNPs were analyzed for over-representation in KEGG pathways using Gene Codis (http://genecodis2.dacya.ucm.es/). “Wnt signaling pathway” ranked 1st with 9 genes and corrected chi square p = 8.4x10-21. Other highly ranked pathways included “Hedgehog signaling pathway”, “Osteoclast differentiation”, “Focal adhesion” and “Endocrine and other factor-regulated calcium reabsorption”. Interestingly, the pathway “Vascular Smooth Muscle Contraction” also emerged as significant (corrected chi square p = 3.9x10-3).

Discussion

The current analyses of combined GWAS data from more than 250,000 individuals demonstrated genetic overlap between BMD and associated CVD risk factor phenotypes. This indicated that some of the co-morbidity observed in epidemiological and clinical studies may be due to shared risk gene variants. Based on the polygenic enrichment we identified 65 novel BMD loci (26 for FN BMD and 47 for LS BMD) not previously reported. Many of these loci are associated with genes that were validated in our expression assay. The high confirmation rate of the current FDR approach and the association to gene expression assay suggest these loci for follow-up analysis.

By comparing GWAS and gene expression profiling of bone, we can suggest which transcriptional regulators drive the expression of the suggested genes identified in this study. Bone remodeling continues throughout life and involves the fine balance between bone building osteoblasts and resorbing osteoclasts. The complexes NFATc1 and NFkB (including p65/RelA) can function as heterodimers and DNA binding transcriptional activators [41] and are central to osteoclast development and differentiation. They do, however, also have an important function in osteoblasts. Strontium ranelate was shown to increase NFATc1 transactivation in osteoblasts promoting increased expression of WNT3A and WNT5A as well as beta-catenin transcription in osteoblasts [42,43]. This positions NFATc1 activation upstream of canonical and non-canonical Wnt signaling pathways, networks whose interactions and strong associations to bone and metabolism are clearly underscored in the present work. NFATc1 activation is also pathogenetically associated with blood pressure via binding to promoter elements on endothelin-1 (ET-1) thereby regulating its expression [44]. ET-1 regulates salt excretion in the kidney collecting duct [45]. Through regulation of salt excretion, NFATc1 also has a role in mineral metabolism, and thus possibly also affecting the body’s Ca++ balance and metabolism. NFATc1 blockade has been shown to completely prevent oxidized LDL-induced osteogenic transformation of human coronary artery smooth muscle cells as well as oxidized LDL-induced stimulation of osteoblast differentiation [46]. NFATc1 may therefore be a master regulator contributing to predisposition in several of these conditions. Interestingly, the application of this approach has uncovered a uniquely rich and coherent gene network which fully reflects the biological relationship between NFATc1 and the Wnt signaling pathways governing osteoclast/osteoblast activity and engagement in metabolism. Future work should focus on the identification of surrogate markers (transcripts and proteins) of aberrant NFATc1 activation, which in combination with genotyping could provide more accurate risk predictors for the range of conditions affected by this important transcription factor. Vascular smooth muscle contraction was identified as significantly affected among the BMD associated genes. This process is relevant to bone because the contractile elements used in muscle are also a characteristic feature of the osteocytes which constitute 90–95% of bone cells [47], and are dynamic star shaped cells with stretching and contracting protrusions [48]. It is not known if the mechanisms for osteocyte motility are more characteristic to smooth or striated muscle. However, both smooth and striated muscle share common features with osteocytes[49], and muscle-related gene expression in bone has been shown to be affected in postmenopausal osteoporotic women [39] as well as in human iliac bone with reduced BMD due to primary hyperparathyroidism [50].

T1D and T2D are complex metabolic disorders with multiple possible interactions with BMD. However, our results are only to a minor degree influenced by these disorders, since only one of the 26 novel FN BMD associated SNPs has diabetes (T1D) as the driving phenotype (Table 1 and Fig 2) and only 9 (~10%) of the novel the LS BMD associated SNPs has T1D or T2D as the driving phenotype (S3 Table and S3 Fig).

Our results confirm the feasibility of using a genetic epidemiology framework that leverages overlap in genetic signal from independent GWAS of correlated phenotypes for revealing genetic basis of complex phenotypes/diseases. The increased power provided by additional GWAS of associated phenotypes together with the FDR method, roughly doubled the previous number of BMD associated loci [25]. Using the same methods for functional validation of the current findings obtained with our statistical approach, we report a similar rate of significantly expressed genes as in the original BMD report [25]. Furthermore, “Role of Osteoblasts and Chondrocyte in Rheumatoid Arthritis” was the top-most significantly affected canonical pathway when subjecting the 155 genes encompassed by all loci at FDR < 0.01 for FN and LS BMD to IPA. This pathway was also among the most significantly affected in a study by Gupta et al.[51], using a Bayesian block-clustering algorithm to analyze GWAS of multiple phenotypes related to bone, thus supporting our results. It should be noted that, when analyzing BMD associated genes by IPA and similar methods, intergenic, and also intragenic SNPs, not necessarily affects transcription of the closest gene. Gene polymorphisms have been shown to affect more distant genes located several Mbp away [52,53]. More detailed experimental validation of the current findings is warranted. Our method for correction of the overlap in some of the GWAS cohorts examined, should exclude contribution from environmental factors. We also controlled for inflation using genomic control correction of each primary single phenotype GWAS. Further, the overlapping loci were spread over all autosomes in the different phenotypes. If a single control group used in several samples were driving the findings, it would be expected that the same region would have been significant across different phenotypes. This is particularly evident in the GWAS of blood lipids, where the same sample was used to discover new genes for three different phenotypes [26], but the pattern of loci was quite different across the different traits. This suggests that the findings are not due to common genetic variation in potentially overlapping control groups.

In conclusion, we identified 26 and 47 novel genomic loci associated with BMD in FN and LS, respectively, by leveraging genetic pleiotropy with several CVD-related traits, including T1D, T2D, SBP, DBP, LDL, TG, WHR and HDL. Association analyses point to genes involved in metabolism and activated immunological pathways. The results warrant further experimental investigations to clarify the clinical implications, and could lead to improved screening programs and prevention strategies.

Conditional FDR 2-D lookup table for femoral neck BMD.

Based on the combination of p-value for the SNPs in femoral neck BMD (PBMD) and that of the pleiotropic trait: A. type 1 diabetes (T1D), B. type 2 diabetes (T2D), C. systolic blood pressure (SBP), D. diastolic blood pressure (DBP), E. high density lipoprotein (HDL), F. low density lipoprotein (LDL), G. triglycerides (TG), and H. waist hip ratio (WHR) we assigned a conditional FDR value to each SNP associated with femoral neck BMD, by interpolation into a 2-D look-up table. Color scale refers to the conditional FDR values.

(TIF)

Click here for additional data file.

Conditional FDR 2-D lookup table for Lumbar Spine BMD.

Based on the combination of p-value for the SNPs in lumbar spine BMD (PBMD) and that of the pleiotropic trait: A. type 1 diabetes (T1D), B. type 2 diabetes (T2D), C. systolic blood pressure (SBP), D. diastolic blood pressure (DBP), E. high density lipoprotein (HDL), F. low density lipoprotein (LDL), G. triglycerides (TG), and H. waist hip ratio (WHR), we assigned a conditional FDR value to each SNP associated with lumbar spine BMD, by interpolation into a 2-D look-up table. Color scale refers to the conditional FDR values.

(TIF)

Click here for additional data file.

Conditional FDR Manhattan plots for lumbar spine BMD.

‘Conditional Manhattan plot’ of conditional–log10 (FDR) values for bone mineral density (BMD, lumbar spine) alone (small black dots) and BMD given the associated phenotypes type 1 diabetes (T1D; BMD|T1D), type 2 diabetes (T2D; BMD|T2D), systolic blood pressure (SBP; BMD|SBP), diastolic blood pressure (DBP; BMD|DBP), high density lipoprotein (HDL; BMD|HDL), low density lipoprotein (LDL; BMD|LDL), triglycerides (TG; BMD|TG), and waist hip ratio (WHR; BMD|WHR). SNPs with conditional–log10 FDR > 2 (i.e. FDR < 0.01) are shown with large points. A black line around the large points indicates the most significant SNP in each LD block and this SNP was annotated with the closest gene, which is listed above the symbols in each locus. Gene symbols were obtained from HGNC gene databases and colored in line with the second phenotype, which gives the minimal conditional FDR value. Genes previously reported by other studies were marked by stars (*).

(TIF)

Click here for additional data file.

Genetic pleiotropy enrichment.

Conditional Q-Q plot of nominal versus empirical -log10 p-values (corrected for inflation) in bone mineral density (BMD, femoral neck) below the standard GWAS threshold of p < 5x10-8 as a function of significance of association with CVD risk factors, including type 1 diabetes (T1D), type 2 diabetes (T2D), low density lipoprotein (LDL) and waist hip ratio (WHR) at the level of -log10(p) ≥ 0 (all SNPs),–log10(p) ≥ 1,–log10(p) ≥ 2,–log10(p) ≥ 3 corresponding to p ≤ 1, p ≤ 0.1, p ≤ 0.01, p ≤ 0.001, respectively. Dotted lines indicate the null-hypothesis.

(TIF)

Click here for additional data file.

QQ plots for Lumbar Spine-BMD.

Conditional Q-Q plot of nominal versus empirical -log10 p-values (corrected for inflation) in bone mineral density (BMD, lumbar spine) below the standard GWAS threshold of p < 5x10-8 as a function of significance of association with A. type 1 diabetes (T1D), B. type 2 diabetes (T2D), C. systolic blood pressure (SBP), D. diastolic blood pressure (DBP), E. high density lipoprotein (HDL), F. low density lipoprotein (LDL), G. triglycerides (TG), and H. waist hip ratio (WHR) at the level of -log10(p) ≥ 0 (all SNPs),–log10(p) ≥ 1,–log10(p) ≥ 2,–log10(p) ≥ 3 corresponding to p ≤ 1, p ≤ 0.1, p ≤ 0.01, p ≤ 0.001, respectively. Dotted lines indicate the null-hypothesis.

(TIF)

Click here for additional data file.

Conditional QQ plot for Femoral neck BMD on CAD.

Conditional Q-Q plot of nominal versus empirical -log10 p-values (corrected for inflation) in bone mineral density (BMD, femoral neck) below the standard GWAS threshold of p < 5x10-8 as a function of significance of association with Coronary Artery Disease (CAD) at the level of -log10(p) ≥ 0 (all SNPs),–log10(p) ≥ 1,–log10(p) ≥ 2,–log10(p) ≥ 3 corresponding to p ≤ 1, p ≤ 0.1, p ≤ 0.01, p ≤ 0.001, respectively. Dotted lines indicate the null-hypothesis.

(TIF)

Click here for additional data file.

Conditional QQ plot for lumbar spine BMD on CAD.

Conditional Q-Q plot of nominal versus empirical -log10 p-values (corrected for inflation) in bone mineral density (BMD, lumbar spine) below the standard GWAS threshold of p < 5x10-8 as a function of significance of association with Coronary Artery Disease (CAD) at the level of -log10(p) ≥ 0 (all SNPs),–log10(p) ≥ 1,–log10(p) ≥ 2,–log10(p) ≥ 3 corresponding to p ≤ 1, p ≤ 0.1, p ≤ 0.01, p ≤ 0.001, respectively. Dotted lines indicate the null-hypothesis.

(TIF)

Click here for additional data file.

Details of Statistical Analysis

(DOC)

Click here for additional data file.

Summary data from all GWAS used in the current study

(DOCX)

Click here for additional data file.

All identified loci associated with femoral neck BMD

(DOCX)

Click here for additional data file.

Identified loci containing novel SNPs or genes associated with lumbar spine BMD

(DOCX)

Click here for additional data file.

Identified loci containing known SNPs or genes associated with lumbar spine BMD

(DOCX)

Click here for additional data file.

Gene titles and gene ontology function terms of genes associated with LS an FN BMD loci at FDR <0.01

(DOCX)

Click here for additional data file.