Osteoporosis is a skeletal disorder characterized by compromised bone strength leading to an increased risk of fracture 1. In the United States alone, there are over 1.5 million fractures each year, including 280,000 hip fractures and 500,000 vertebral fractures. According to the recent U.S. Surgeon General's Report on Skeletal Health, fractures remain a large and growing public health concern 2. Presently, the gold standard for assessment of fracture risk is measurement of bone mineral density (BMD, g/cm²) by dual-energy X-ray absorptiometry (DXA). Whereas low BMD is among the strongest risk factors for fracture 34, a number of clinical studies have demonstrated that other measurements, such as quantitative ultrasound (QUS) and bone geometry, are important for fracture prediction and osteoporosis treatment monitoring 567. Thus, QUS of the calcaneus is associated with hip fracture, largely independent of BMD 89. A growing body of evidence in recent years indicates that femoral geometry also contributes importantly to hip fracture risk 1011.

BMD, QUS, and femoral bone geometry are approximately normally distributed, complex traits. A wealth of studies have documented BMD to be under strong genetic control with 50–70% heritability 1213. Similarly, QUS 81415 and hip geometry 16171819 are probably regulated by additive genetic factors. However, despite years of research in the field of osteoporosis genetics, progress to date has been modest in successfully identifying major genes determining BMD, QUS, and bone geometry in the general population. The numbers of quantitative trait loci (QTLs) and genes linked and/or associated with osteoporosis-related phenotypes continue to expand and the list has become considerably more detailed and complex. More than 20 genome linkage scans to date have revealed multiple QTLs covering all chromosomes but the Y chromosome 20. Moreover, the results from one study have inconsistently been replicated in other samples. To overcome these obstacles, a collaborative meta-analysis of 9 genome-wide linkage searches of BMD was recently conducted, including data from 11,842 subjects, members of 3,045 families 13. The meta-analysis suggested a number of specific QTLs to be pursued further (1p13.3-q23.3 and 1q32-q42.3, 3p25.3-p22.1, 11p12-q13.3, 12q24.31-qter, and 18p11-q12.3). An additional factor that adds to the complexity of finding genes for osteoporosis-related traits is the notion (including our own 2122) that QTLs for bone density and geometry 23 are skeletal site-specific, age-group-, and sex-specific 2425.

Numerous biological candidate genes for bone phenotypes have also been proposed, but few have been validated with large-scale evidence 13. Recent meta-analyses for estrogen receptor alpha 26, collagen type I alpha 1 27, methylenetetrahydrofolate reductase, 28 and vitamin D receptor 29 genes demonstrated that the polymorphisms in these genes each explain a small percentage of the variation in BMD or fracture. With the near completion of the International HapMap Project and the rapid improvements in high throughput genotyping technology, the ultimate understanding of the genetic basis of osteoporosis may come from a genome-wide association approach in which the whole human genome is surveyed for common genetic variation in osteoporosis-related heritable quantitative traits such as BMD 19. To further advance the field of skeletal genetics, efforts are now turning toward studies that are able to cover the vast majority of the genome with a dense set of SNPs, in order to identify the genes linked and/or associated with osteoporosis-related traits. We thus performed genome-wide linkage and association analyses using a dense genome scan in members of extended pedigrees from the Framingham Heart Study (FHS).

Bone-related quantitative traits and derived phenotypes analyzed in the FHS 100K SNP resource are listed in Table 1. Each trait was analyzed as multivariable-adjusted residuals from cohort- and sex-specific models; in addition, geometric traits were adjusted for age, age² only. Among these, the following ten multivariable-adjusted residuals are the primary phenotypes presented here: BMD of the femoral neck (FNBMD; n with data = 1141), trochanter (TRBMD; n = 1141) and the second to fourth vertebral bodies of the lumbar spine (LSBMD; n = 1127); a quantitative ultrasound measure of the calcaneus, broadband ultrasound attenuation (BUA; n = 1105); femoral geometry measures at the "narrow neck" region of the hip and the femoral shaft including narrow neck section modulus (NeckZr; n = 1107) and width (NeckWr; n = 1106), femoral shaft section modulus (ShaftZr; n = 1037) and width (ShaftW; n = 1037), femoral neck-shaft angle (NSA; n = 1107), and femoral neck length (NeckLeng; n = 1101). Additional bone-related quantitative phenotypes are available in the FHS 100K SNP resource (web posted at http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?id=phs000007) including calcaneal speed of sound; several hip geometry measures including neck cross-sectional moment of inertia; neck, shaft and intertrochanteric buckling ratio, and shaft cross sectional moment of inertia from these same examinations.

Heritability estimates for all bone phenotypes ranged from 30–66%. We consider here only multivariable-adjusted, residual phenotypes derived from ten primary traits defined above. There were 12 associations with 100K SNPs in GEE models at p < 0.000001 and 2 associations in FBAT models with p < 0.000001. The 25 SNPs with the lowest p-values in GEE (p < 3.5 × 10^-6) and FBAT models (p < 2.5 × 10^-5), ordered by chromosomal position, are shown in Tables 2a and 2b, respectively. Of these SNPs, 10 (GEE) and 19 (FBAT) were within 60 kb of a known gene. Several SNPs associated with phenotypes were in high pairwise LD (r² ≥ 0.8). When this occurred we displayed results for the higher ranking SNP based on a combination of FBAT and GEE p-values. Of note, some SNPs were associated with more than one trait: thus, rs10514345 on chromosome 5 was associated (GEE) with NeckZ1 (combined sample) and LSBMD (females) and rs1209921 on chromosome 21 was associated (FBAT) with TRBMD (combined sample) and FNBMD (females) (Tables 2a and 2b, respectively).

In Table 2c we present results with LOD scores above 2.0 from genome-wide linkage analyses with their corresponding 1.5-LOD confidence interval. There were linkage peaks observed for the majority of chromosomes, and in general there were a greater number of linkages with bone geometry traits than with either BMD or ultrasound measures. LOD scores ≥3.0 were found on chromosomes 15 (1.5 LOD support interval, 51,336,679–58,934,236 bp) and 22 (35,890,398–48,603,847 bp), with femoral shaft section modulus. There were no linkage regions of overlap with markers that were found to be associated with bone phenotypes in the FBAT and GEE analyses. Notably, when we restricted our sample to members of the Offspring Cohort, LOD scores for BMD phenotypes increased; for example, for LSBMD on chr. 7 (63,211,216 bp) from 1.25 to 2.43; for FNBMD on chr. 9 (104,660,212 bp) from 1.79 to 2.65, and on chr. 17 (63,038,023 bp) from 1.83 to 2.43 (not shown in the Table 2c).

In order to prioritize SNPs potentially associated with multiple phenotypes, we evaluated four phenotypic subgroups by both population- and family-based tests. The first subgroup was BMD (femoral neck, trochanter and spine), the second was BMD + bone ultrasound (BUA), the third was bone geometry (femoral neck-shaft angle, femoral neck length, femoral neck width, and shaft width), and the fourth, hip BMDs with hip geometry. For each SNP we calculated the number of traits in the subgroup that were significantly associated with the SNP at p < 0.01 in both GEE and FBAT, and then ranked SNPs with highest numbers of nominally significantly associated traits in which the direction of association was the same. In cases where SNPs had the same proportions of nominally significant traits, SNPs were additionally ranked by the mean of GEE p-values across traits.

In Table 3a we present only the SNPs selected based on associations with p-value of both GEE and FBAT < 0.01 with the first subgroup (BMD phenotypes only), for a combined sample of men and women. Using the annotation by UCSC genome browser tables http://genome.ucsc.edu/, of the top 20 SNPs for the BMD subgroup, 9 (45%) were identified in genes or within 60 kb of a known gene, and none were coding (Table 3a). The nominally significant intragenic SNPs were found in several genes, including cadherin 9 type 2 (CDH9) on 5p14 and gene deleted in colorectal carcinoma (DCC) on 18q21. It is also noteworthy that the association of TRBMD with rs1209926 on chr. 21 (39,078,994 bp) occurred in the region identified in our original genome-wide linkage study in the FHS 22 SNP rs2938392 in the peroxisome proliferator-activated receptor gamma (PPARG) was also associated with the BMD traits in GEE tests and less strongly in FBAT analyses (not shown in Table 3a).

We also examined men and women separately within the first phenotypic subgroup and confirmed the associations of BMD phenotypes, observed in the combined sample, with rs922028 (chromosome 4, no known gene), rs1479679 (in CDH9, chr. 5), and rs768207 (in DCC, chr. 18), all in women (Table 3b). There were also associations in one sex that were neither observed in the other sex nor in the combined sample, and men and women did not share many polymorphisms associated with BMD.

We also looked for evidence of pleiotropy by evaluating association results within the other three pre-specified phenotype subgroups, BMD + BUA and hip geometry, as well as combinations of hip BMD with hip geometry (fourth subgroup). Notably, pleiotropic associations between BMD and other traits were uncommon (not shown in a table). When considering BMD and BUA together (phenotype subgroup 2), there were nominally significant associations (p < 0.0005) with a SNP rs953934 in the interleukin 1 receptor-like 1 gene (IL1RL1) on 2q12 (not shown). Of the top SNPs for the Hip Geometry subgroup (NSA, NeckLeng, NeckW1r and ShaftW1), 50% were identified in genes or within 60 kb of a known gene; none were coding. The nominally significant intragenic SNPs were found in several genes, including rs7151976 in heat shock 70 kD protein 2 (HSPA2) on 14q24 and rs2834719 in runt-related transcription factor 1 (RUNX1, or CBFA2) on 21q22. Significant associations of SNPs with BMD phenotypes did not overlap with geometric phenotypes, and within the hip geometry subgroup, SNPs significant for one geometric phenotype did not overlap with other phenotypes within this subgroup (data not shown).

We further queried the FHS 100K array results, to verify, in silico, whether previously reported candidate genes for osteoporosis were associated with BMD, QUS and bone geometry in our sample. Thus Table 4 shows SNPs in these pre-defined candidate genes that were associated at p ≤ 0.05 in either GEE or FBAT models. We selected the following six previously reported genes that had been extensively studied: COL1A1, CYP19, ESR1, LRP5, MTHFR, and VDR. For the MTHFR gene, we found 1 non-synonymous coding SNP (rs1801133) that was associated with femoral shaft section modulus and neck-shaft angle by GEE and FBAT, respectively. We 38 have previously reported an association of this SNP with BMD in unrelated participants in the Framingham Osteoporosis Study and others have confirmed this association 2839. None of the three 100K SNPs in the ESR1 gene (rs1884052, rs3778099 and rs3866461) with nominally significant associations with bone phenotypes were in LD with widely published SNPs (rs2234693 and rs9340799). A single intronic 100K SNP from LRP5 (rs4988300) was associated with shaft section modulus in males and with femoral neck BMD in females. This 100K SNP was not in LD with published rs3736228 in LRP5 40. There were no LD data available on the HapMap website for the SNPs widely studied by others with respect to BMD (rs4988321 and rs627174). For the VDR gene there was a single 100K SNP (rs2189480) nominally significantly associated with femoral neck section modulus and spine BMD; however, there were no data available on LD between this 100K SNP and the 3 out of 5 SNPs reported in a recent meta-analysis of association studies for osteoporosis 29. Two SNPs reported in that meta-analysis (rs1544410 and rs11568820) were not in LD with the one 100K SNP studied here. There were two 100K SNPs in the CYP19 gene (rs10519297 and rs2008691) associated with femoral neck-shaft angle that were in LD with at least one of the previously reported SNPs studied with respect to bone phenotypes (rs700518, rs4775936, rs10046, and rs11575899). Also, there was a single 100K SNP in the COL1A1 gene (rs2075555) that was associated with femoral neck width in women and shaft width in men; however, there were no data at the HapMap website regarding LD with the three other SNPs in this gene reported by other osteoporosis research groups (rs1107946, rs11327935, rs1800012).

To our knowledge, this is the first genome-wide association study of many of the most commonly measured osteoporosis-related phenotypes. This unbiased approach allowed simultaneous exploration of multiple genetic loci rather than focusing on a single-gene association. Using a strategy that combined GEE and FBAT association tests within the subgroup of BMD traits, 40 SNPs were each found to be associated with several bone phenotypes. More than half of these were in or near several genes, including cadherin 9 type 2 (CDH9) on 5p14, nuclear receptor subfamily 5, group A gene (NR5A2) at 1q32, gene deleted in colorectal carcinoma (DCC) on 18q21, and PPARG on 3p25. Some of the genes found to be associated with one or several bone density phenotypes have not previously been studied for osteoporosis. Analyses stratified by sex confirmed association of BMD phenotypes with CDH9 and DCC in women. The sex specificity of our findings is fully in line with the previous observations of both sex-specific linkage 1321 and sex-specific genetic association with BMD 2541.

Our phenotypic subgroup combining BMD and the ultrasound trait, BUA, demonstrated that the interleukin 1 receptor-like 1 (IL1RL1) gene on 2q12 was associated with all bone mass traits. In contrast there was a notable absence of pleiotropic associations between BMD and hip geometry (femoral neck length and width, neck-shaft angle, and shaft width). This is an important observation in light of the current debates about the role of femoral bone geometry in fracture prediction above and beyond BMD, as well as a criticism of DXA-derived bone geometric traits, which are perceived as partially captured in the BMD measurement. We have clearly shown that BMD traits (measured at the femoral neck, trochanter, and at the lumbar spine) share associations with SNPs in several genes, whereas there are no such SNPs that share associations with both hip BMD and hip geometric traits. Despite substantial shared genetic determinants between hip BMD and geometric traits (preliminary unpublished results from our group), these phenotypes are at least in part governed by distinct gene variants. Since genes contributing to variation in BMD do not always contribute to osteoporotic fractures 42, there is a necessity to genetically dissect geometric phenotypes along with BMD to derive a valid endophenotype of osteoporosis and to better encircle the complex heritability of this disease 31.

When we queried the FHS 100K association results, several SNPs in previously well-studied candidate genes for osteoporosis were identified (FBAT or GEE p < 0.05), such as rs1801133 in MTHFR; rs1884052 and rs3778099 in ESR1; rs4988300 in LRP5; rs2189480 in VDR; rs2075555 in COLIA1; rs10519297 and rs2008691 in CYP19. We also examined SNPs in other candidate genes of interest to our group such as SNPs in PPARG (rs10510418 and rs2938392) that were associated with BMD and ultrasound and SNPs in ANKH (rs2454873 and rs379016) that were associated with femoral neck section modulus.

We observed limited overlap between the results of the association and linkage analyses, which may, in part, be attributed to a low power of the linkage analysis. The resolution of linkage peaks depends not only upon the extracted information content of the genotyped markers, but also upon the useful recombination information contained in the sample (which is related to family structure and sample size). Despite low power for detecting modest genetic effects in linkage analyses, we identified several loci with a substantial indication of linkage, such as LOD scores ≥3.0 on chromosomes 15 (around 55,573,538 bp) and 22 (around 44,039,280 bp) for bone geometry phenotypes. Also, LOD scores for BMD phenotypes generally increased after restricting our sample to the members of the Offspring Cohort only, despite potential loss of power due to smaller sample size. This observation indirectly supports our previous finding of refined LOD scores in subgroups of age 21 and may be explained by the younger age of the Offspring Cohort members. The BMD phenotypes in that cohort might be less affected by aging-related artifacts in DXA scans, such as osteophytes, ossification of ligaments, and abdominal aortic calcification. We also observed little overlap between the results of the linkage analysis reported above and our previous studies of the same phenotypes 2331, which may, in part, be attributed to low power of these linkage analyses and limited overlap between the 100K sample and the relatively large sample of FHS pedigrees published before 2331.

To our knowledge, this is a first genome-wide association study of many of the most commonly measured osteoporosis-related phenotypes. It is distinct in many respects. First, this is a population-based sample rather than a study of affected families. Second, our study included adults from two generations of the same families. The longitudinal nature of the FHS and ongoing recruitment of the third generation (grandchildren of the Original Cohort) will enable us to validate the results of this study with an even larger sample of 3 generations. Third, results of this collaborative effort (FHS 100K SNP study) are publicly available on the web http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?id=phs000007, which makes it a unique resource contributing to ongoing efforts of other groups.

One of the strengths of this study lies in the unbiased approach that allowed simultaneous exploration of multiple genetic loci rather than focusing on a single-gene association. However osteoporosis is likely to be a polygenic condition with a strong environmental component. Therefore, an obvious limitation of this study is the absence of analyses testing gene-environment and gene-by-gene interactions 43, which is a demanding element of study design and statistical analyses. We did not formally test for associations of the SNPs with some known environmental factors, such as smoking or estrogen replacement therapy in women, for which FHS has well-documented longitudinal data, because of sample size considerations.

A second potential limitation of this study is that there may be unmeasured factors that would affect the association or impact on generalizability of our results. Sampling of homogenous individuals is needed to minimize these effects. We attempted to approach this problem in part by stratifying our sample by gender; however, a more involved stratification, such as younger men or postmenopausal women, was not performed due to sample size restrictions. A third potential shortcoming of this study is the limited coverage of the genome by the 100K SNP array. These SNPs are undersampled from coding regions (both non-synonymous and synonymous) and oversampled from regions outside genes, relative to SNPs in the overall HapMap database 44. Current plans of the FHS include the genotyping and analysis of over half a million genetic variants in approximately 9,000 men and women from three generations of the FHS. This higher density genotyping will provide wider coverage of genomic variation and thus will clearly have advantages for GWA and linkage studies in the future.

Fourth, at the present there is no clear-cut strategy for prioritization of the regions identified by our approach to whole genome analyses. We implemented one possible SNP prioritization strategy based on a somewhat arbitrary grouping of multiple phenotypes. Pleiotropy across different (but probably related) phenotypes may be an indication of genes that influence more than one biological pathway and thus may affect overall risk of osteoporosis.

Finally, we are hesitant to regard any result reported in this paper as genome-wide significant. Although observing evidence for association at a significance level of .0001 in both the GEE and FBAT test is very encouraging, it does not meet the standard of genome-wide significance. Not only do we have a large number of statistical tests for each phenotype, but we also have numerous phenotypes 30. None of our findings should be considered true positive results until they are replicated in independent samples and/or functional studies are performed that provide the mechanism of the SNP effect.

An unbiased genome-wide strategy to detect genetic associations with skeletal traits provides an opportunity to identify novel pathways for replication in other human populations. This study is thus hypothesis-generating and awaits replication in comparable samples. On one hand, ample data on genetic markers and associated traits may form the basis of new phenotype definitions, aimed at defining phenotypic groupings that are characterized by the sharing of genetic linkage and association. Such data may uncover biological mechanisms not obvious at the phenotypic level 42. On the other hand, a composite genetic risk score that combines the effects of multiple loci may prove to be of more practical utility for a common and complex chronic disease such as osteoporosis, 43 rather than results of single-gene association studies. Identification of such a molecular profile and creation of an array for genetic screening will be a very useful tool for prevention and management of osteoporosis.

BMD = bone mineral density; Chr = chromosome; DXA = dual-energy X-ray absorptiometry; FBAT = family based association testing; GEE = generalized estimating equations; GWA = Genome wide association; HWE = Hardy-Weinberg equilibrium; IBD = identity-by-descent; LD = linkage disequilibrium; LOD = logarithm of the odds; MAF = Minor allele frequency; QUS = bone quantitative ultrasound; SNP = single nucleotide polymorphism; VCA = variance component analysis.

The authors declare that they have no competing interests.

All authors have made substantial contributions to conception and design or acquisition of phenotypic data. JD participated in the design of the study and performed statistical analyses; DPK, DEK, SD, JMM, and KLL contributed to the analysis and interpretation of data. DK and DPK drafted the manuscript and revised it critically for important intellectual content. All authors read and approved the final manuscript.