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 ã 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][6][7][8][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][16][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][22][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][33][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.

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  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 (-log 10 (p) ! 0,-log 10 (p) ! 1,-log 10 (p) !2,-log 10 (p) !3 corresponding to p < 1, p < 0.1, p < 0.01, p < 0.001, respectively). The nominal p-values (-log 10 (p)) are plotted on the y-axis, and the empirical quantiles (-log 10 (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-log 10 (p) < 7.3 (corresponding to p > 5x10 -8 ).
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][33][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 Conditional Q-Q plot of nominal versus empirical -log 10 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 -log 10 (p) ! 0 (all SNPs),-log 10 (p) ! 1,-log 10 (p) ! 2,-log 10 (p) ! 3 corresponding to p 1, p 0.1, p 0.01, p 0.001, respectively. Dotted lines indicate the null-hypothesis. are as small as or smaller than the observed p-values, FDR(P 1 jP 2 ) = π 0 (P 2 )P 1 /F(P 1 jP 2 ), where F (P 1 jP 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 jP 2 ) were conservatively estimated by setting π 0 (P 2 ) equal one and replacing F(P 1 jP 2 ) by empirical conditional cdf. The conditional FDR values for BMD on second pleiotropic traits (denoted by FDR BMD , 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 (-log 10 (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-r 2 > 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.
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-r 2 > 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-r 2 > 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).
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 log 2 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).

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 Detection and Verification
The previous study of BMD related SNPs by Estrada et al.   Table 1 and S1  [25]. Also in the cross stage (I and II) analyses, all but 5 loci for FN BMD and 8 loci for LS were successfully reidentified (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).
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, NFKB and SMAD3 as central nuclear hubs were generated (Fig 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  AREGB, AXIN1, BMP7, CTNNB1,  FAM20C, JAG1, KLF4,LGR4, LRP5, MEF2C,  NFATC1, PKDCC, RELA, SFRP4, SMAD3,SOST,  SOX9, SP7, SPP1, TNFRSF11A, TNFRSF11B  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  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. 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 -log 10 (p) ! 0 (all SNPs),log 10 (p) ! 1,-log 10 (p) ! 2,-log 10 (p) ! 3 corresponding to p 1, p 0.1, p 0.01, p 0.001, respectively. Dotted lines indicate the null-hypothesis. (TIF) S7 Fig. Conditional QQ plot for lumbar spine BMD on CAD. Conditional Q-Q plot of nominal versus empirical -log 10 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 -log 10 (p) ! 0 (all SNPs),log 10 (p) ! 1,-log 10 (p) ! 2,-log 10 (p) ! 3 corresponding to p 1, p 0.1, p 0.01, p 0.001, respectively. Dotted lines indicate the null-hypothesis. (TIF) S1 File. Details of Statistical Analysis (DOC) S1 Table. Summary data from all GWAS used in the current study (DOCX) S2 Table. All identified loci associated with femoral neck BMD (DOCX) S3 Table. Identified loci containing novel SNPs or genes associated with lumbar spine BMD (DOCX) S4 Table. Identified loci containing known SNPs or genes associated with lumbar spine BMD (DOCX) S5 Table. Gene titles and gene ontology function terms of genes associated with LS an FN BMD loci at FDR <0.01 (DOCX)