Abstract

Quantitative ultrasound of the heel captures heel bone properties that independently predict fracture risk and, with bone mineral density (BMD) assessed by X-ray (DXA), may be convenient alternatives for evaluating osteoporosis and fracture risk. We performed a meta-analysis of genome-wide association (GWA) studies to assess the genetic determinants of heel broadband ultrasound attenuation (BUA; n = 14 260), velocity of sound (VOS; n = 15 514) and BMD (n = 4566) in 13 discovery cohorts. Independent replication involved seven cohorts with GWA data (in silico n = 11 452) and new genotyping in 15 cohorts (de novo n = 24 902). In combined random effects, meta-analysis of the discovery and replication cohorts, nine single nucleotide polymorphisms (SNPs) had genome-wide significant (P < 5 × 10−8) associations with heel bone properties. Alongside SNPs within or near previously identified osteoporosis susceptibility genes including ESR1 (6q25.1: rs4869739, rs3020331, rs2982552), SPTBN1 (2p16.2: rs11898505), RSPO3 (6q22.33: rs7741021), WNT16 (7q31.31: rs2908007), DKK1 (10q21.1: rs7902708) and GPATCH1 (19q13.11: rs10416265), we identified a new locus on chromosome 11q14.2 (rs597319 close to TMEM135, a gene recently linked to osteoblastogenesis and longevity) significantly associated with both BUA and VOS (P < 8.23 × 10−14). In meta-analyses involving 25 cohorts with up to 14 985 fracture cases, six of 10 SNPs associated with heel bone properties at P < 5 × 10−6 also had the expected direction of association with any fracture (P < 0.05), including three SNPs with P < 0.005: 6q22.33 (rs7741021), 7q31.31 (rs2908007) and 10q21.1 (rs7902708). In conclusion, this GWA study reveals the effect of several genes common to central DXA-derived BMD and heel ultrasound/DXA measures and points to a new genetic locus with potential implications for better understanding of osteoporosis pathophysiology.

INTRODUCTION

Bone structure in vivo has largely been evaluated using the attenuation of a photon beam by hydroxyapatite, the principal mineral in bone. This is positively related to the mass of hydroxyapatite in the path of the beam conventionally termed bone mineral content and normalized to bone area to produce an entity termed areal bone mineral density (BMD). To allow for the reduced attenuation of the beam by overlying non-bone tissues in central areas of the body, two photon beam energies are used, resulting in a clinical technique termed dual-energy X-ray absorptiometry (DXA), which at peripheral skeletal sites is termed pDXA.

Over the past 60 years, ultrasonic material analysis has been developed as a method of determining material properties of a variety of structures. In the last 30 years, this methodology has been applied to the in vivo assessment of bone structure and fragility termed quantitative ultrasound (QUS). This consists of the use of two separate ultrasound measurement techniques, velocity of sound (VOS) and broadband ultrasound attenuation (BUA). While much remains to be discovered about the exact physical determinants of QUS measures in the intact living calcaneum (1), cadaver studies have established a strong correlation of such indices with bone quantity and trabecular structure (2). Assessment of bone properties in the heel using QUS can predict the risk of prevalent osteoporotic fractures, such as those in the spinal vertebrae, comparably with DXA of the spine or hip, the so-called gold standard clinical techniques (3–5). Pearson correlation coefficients of heel QUS or pDXA with central DXA of the hip or spine in population-based studies are modest, typically in the range of 0.4–0.6 (6). Moreover, twin- and family-based studies have found genetic correlations of the order of 0.3–0.6 and environmental correlations of the order of 0.1–0.3 (7–9); yet relative risk estimates for fracture using QUS are of similar magnitude to those derived from central DXA (5,10,11). A recent meta-analysis showed that heel QUS predicts risk of various fractures (hip, vertebral and any clinical fractures) independently from hip BMD (12). Overall, these results suggest that QUS of the calcaneum might capture additional genetic determinants of bone structure beyond those associated with central DXA.

A genetic contribution to osteoporosis is well established with heritability estimates reaching 84% for central BMD (13), 74% for heel QUS (7,14), 47% for bone loss (15), and 48% for hip fracture (16). Previous genome-wide association (GWA) studies have identified several chromosomal regions associated with BMD in the hip and lumbar spine regions (17,18). The most recent meta-analysis of GWA studies, performed in the context of the Genetics Factors for Osteoporosis (GEFOS) consortium, identified 56 genome-wide significant loci (32 new) associated with hip/spine BMD (19). Fourteen out of these 56 BMD-associated loci were also associated with fracture risk in a case–control meta-analysis involving ∼31 000 fracture cases among 133 000 individuals (19). Using data from the GEFOS consortium, we aimed to extend the findings for central DXA-derived BMD phenotypes by searching for single nucleotide polymorphisms (SNPs) associated with heel QUS or heel DXA measures across the human genome.

RESULTS

Participant characteristics are summarized in Table 1 and key features of the discovery and replication phases are summarized in Figure 1. In aggregate, the initial discovery phase meta-analysis in 11 cohorts (Supplementary Material, Table S1) identified 42 loci of at least suggestive significance in relation to heel bone measures, of which 9 overlapped with loci previously found to be potentially associated with hip or spine BMD in the GEFOS-BMD meta-analysis (19). Regional conditional analyses results were available for QUS measures from 9 cohorts (comprising 7 of the initial discovery cohorts and a further 2 new cohorts that joined later). Based on the results of the conditional analyses (that identified two secondary signals for the QUS measures) and final combined meta-analysis of the unconditional results from all 13 discovery cohorts, a total of 25 independent SNPs (Table 2) were selected for replication in the next phase (i.e. in silico studies and de novo genotyping). Including the two secondary signals, the selected SNPs comprise 15 SNPs that were primarily associated with either BUA or VOS, and 12 SNPs that were associated with heel DXA BMD (Table 2).

Table 1.

Characteristics of studies that contributed to GWAS discovery and replication of SNP associations with heel QUS/DXA BMD measures

Stage\cohort Country Demographics
 
Heel QUS/DXA BMD outcomes 
  N Females Age (years) Weight (kg) Height (cm) BUA (dB/MHz) VOS (m/s) Heel BMD (g/cm2
(%) Mean (SD) Mean (SD) Mean (SD) N Mean (SD) N Mean (SD) N Mean (SD) 
GWAS discovery  
 EPIC UK 2630 56 62.1 (8.6) 80.5 (15.4) 167 (9) 2630 83 (19) 2630 1632 (40) – – 
 FHS USA 3229 58 64.6 (11.9) 76.8 (17.2) 166 (10) 3229 73 (21) 3225 1548 (38) – – 
 HKOS China 730 100 48.7 (15.4) 54.8 (10.4) 155 (7) 730 74 (22) 730 1551 (41) – – 
 NSPHS06 Sweden 495 55 51.4 (19.1) 71.9 (12.8) 164 (10) 495 96 (21) – – – – 
 RSI Netherlands 1615 54 66.5 (8.2) 74.3 (11.8) 169 (9) 1615 112 (13) 1615 1525 (37) – – 
 SHIP Germany 1198 54 58.0 (13.5) 80.2 (15.8) 168 (9) 1198 115 (15) 1198 1565 (35) – – 
 SHIP-TREND Germany 687 56 50.8 (13.6) 78.7 (15.1) 170 (9) 687 116 (14) 687 1571 (33) – – 
 TWINSUK1 UK 1701 100 46.2 (12.1) 65.8 (12.5) 163 (6) 1701 76 (18) 1701 1658 (49) – – 
 TWINSUK23 UK 1975 100 46.9 (12.5) 66.1 (12.2) 163 (6) 1975 76 (18) 1975 1653 (50) – – 
 H2SS Korea 1753 53 60.8 (6.6) 61.9 (10.0) 158 (8) – – 1753 1591 (45) – – 
 AGES Iceland 3179 58 76.4 (5.4) 75.8 (14.3) 167 (9) – – – – 3179 0.491 (0.152) 
 CroatiaKorcula Croatia 878 64 56.3 (14.2) 79.0 (14.2) 168 (9) – – – – 878 0.443 (0.098) 
 CroatiaSplit Croatia 499 57 49.3 (14.7) 80.6 (16.3) 172 (9) – – – – 499 0.459 (0.101) 
 Subtotal  20 569 66 60.3 (11) 73.3 (14.1) 165 (9) 14 260 86 (18) 15 514 1593 (42) 4556 0.478 (0.138) 
In silico replication 
 AOGCa Australia/UKb 1955 100 69.6 (8.6) 69.6 (17.3) 158 (16) – – – – – – 
 B-PROOF Netherlands 1092 59 74.0 (6.7) 76.0 (12.4) 168 (9) 1092 69 (17) 1091 1535 (32) – – 
 HABC USA 1493 48 74.8 (2.9) 73.8 (14.3) 167 (9) 1493 73 (18) 1493 1541 (30) – – 
 MICROS Italy 588 45 46.0 (16.6) 70.2 (14.9) 167 (9) 588 73 (16) 588 1544 (29) – – 
 MrOS-USA USA 3925 73.9 (5.9) 83.1 (12.7) 175 (7) 3925 79 (17) 3925 1551 (30) – – 
 SOF USA 2103 100 80.1 (4.2) 66.3 (12.5) 158 (6) 2103 59 (17) 2103 1527 (30) – – 
 YFS Finland 1265 58 37.9 (5.0) 75.8 (15.5) 172 (9) 1265 80 (16) 1265 1559 (29) 1250 0.560 (0.110) 
 HCS-AUS Australia 986 49 66.2 (7.6) 79.4 (15.5) 166 (9) – – – – 986 0.538 (0.166) 
 Subtotal  13 407 52 69.2 (6.9) 75.4 (14.2) 167 (9) 10 466 73 (17) 10 465 1544 (30) 2236 0.550 (0.138) 
De novo replication 
 AUSTRIOS-B Austria 448 85 83.6 (5.9) 62.0 (12.3) 156 (8) 448 90 (17) 448 1496 (36) – – 
 CABRIO-C Spain 1274 62 62.4 (9.2) 73.7 (13.1) 161 (8) 1274 70 (23) 1273 1545 (41) – – 
 CAIFOS Australia 1113 100 80.0 (2.6) 67.5 (12.1) 157 (6) 1113 101 (9) 1113 1516 (28) – – 
 CALEX-FAM Finland 983 79 37.0 (22.4) 64.3 (16.9) 164 (11) 983 83 (16) – – – – 
 EMAS Europeb 2870 59.9 (11.0) 83.1 (13.6) 173 (7) 2870 80 (19) 2870 1550 (34) – – 
 EPICNOR UK 5723 54 63.6 (9.2) 73.2 (12.4) 167 (9) 5723 79 (20) 5718 1638 (43) – – 
 EPOLOS Poland 684 56 53.4 (16.0) 73.2 (13.7) 166 (10) 684 112 (13) 684 1548 (35) – – 
 FLOS Italy 1000 84 59.8 (12.7) 64.8 (12.3) 163 (9) 1000 58 (7) 1000 1503 (83) – – 
 GEOS Canada 5495 100 55.8 (10.3) 65.4 (11.9) 158 (6) 5495 111 (10) 5495 1546 (32) – – 
 LASA Netherlands 894 51 75.6 (6.5) 74.2 (12.6) 166 (9) 894 71 (20) 894 1611 (44) – – 
 MrOS-SWE Sweden 1718 75.4 (3.2) 80.6 (12.0) 175 (7) 1718 81 (21) 1718 1555 (38) – – 
 OPRA Sweden 821 100 75.2 (0.1) 67.6 (11.3) 160 (6) 821 102 (10) 821 1523 (27) – – 
 OSTEOSII Greece 307 87 50.5 (12.6) 74.1 (15.7) 163 (7) 307 112 (16) 307 1556 (36) – – 
 PEAK25 Sweden 857 100 25.5 (0.2) 64.5 (11.2) 168 (6) 857 118 (11) 857 1575 (32) – – 
 SWS UK 715 100 29.7 (3.7) 72.4 (14.8) 163 (7) 714 72 (13) 715 1548 (27) – – 
 Subtotal  24 902 64 60.2 (10.0) 71.6 (12.7) 165 (8) 24 901 89 (16) 23 913 1568 (40) – – 
Total  58 878 62 62.3 (9.7) 73.0 (13.6) 165 (8) 49 627 85 (17) 49 892 1570 (39) 6792 0.502 (0.138) 
Stage\cohort Country Demographics
 
Heel QUS/DXA BMD outcomes 
  N Females Age (years) Weight (kg) Height (cm) BUA (dB/MHz) VOS (m/s) Heel BMD (g/cm2
(%) Mean (SD) Mean (SD) Mean (SD) N Mean (SD) N Mean (SD) N Mean (SD) 
GWAS discovery  
 EPIC UK 2630 56 62.1 (8.6) 80.5 (15.4) 167 (9) 2630 83 (19) 2630 1632 (40) – – 
 FHS USA 3229 58 64.6 (11.9) 76.8 (17.2) 166 (10) 3229 73 (21) 3225 1548 (38) – – 
 HKOS China 730 100 48.7 (15.4) 54.8 (10.4) 155 (7) 730 74 (22) 730 1551 (41) – – 
 NSPHS06 Sweden 495 55 51.4 (19.1) 71.9 (12.8) 164 (10) 495 96 (21) – – – – 
 RSI Netherlands 1615 54 66.5 (8.2) 74.3 (11.8) 169 (9) 1615 112 (13) 1615 1525 (37) – – 
 SHIP Germany 1198 54 58.0 (13.5) 80.2 (15.8) 168 (9) 1198 115 (15) 1198 1565 (35) – – 
 SHIP-TREND Germany 687 56 50.8 (13.6) 78.7 (15.1) 170 (9) 687 116 (14) 687 1571 (33) – – 
 TWINSUK1 UK 1701 100 46.2 (12.1) 65.8 (12.5) 163 (6) 1701 76 (18) 1701 1658 (49) – – 
 TWINSUK23 UK 1975 100 46.9 (12.5) 66.1 (12.2) 163 (6) 1975 76 (18) 1975 1653 (50) – – 
 H2SS Korea 1753 53 60.8 (6.6) 61.9 (10.0) 158 (8) – – 1753 1591 (45) – – 
 AGES Iceland 3179 58 76.4 (5.4) 75.8 (14.3) 167 (9) – – – – 3179 0.491 (0.152) 
 CroatiaKorcula Croatia 878 64 56.3 (14.2) 79.0 (14.2) 168 (9) – – – – 878 0.443 (0.098) 
 CroatiaSplit Croatia 499 57 49.3 (14.7) 80.6 (16.3) 172 (9) – – – – 499 0.459 (0.101) 
 Subtotal  20 569 66 60.3 (11) 73.3 (14.1) 165 (9) 14 260 86 (18) 15 514 1593 (42) 4556 0.478 (0.138) 
In silico replication 
 AOGCa Australia/UKb 1955 100 69.6 (8.6) 69.6 (17.3) 158 (16) – – – – – – 
 B-PROOF Netherlands 1092 59 74.0 (6.7) 76.0 (12.4) 168 (9) 1092 69 (17) 1091 1535 (32) – – 
 HABC USA 1493 48 74.8 (2.9) 73.8 (14.3) 167 (9) 1493 73 (18) 1493 1541 (30) – – 
 MICROS Italy 588 45 46.0 (16.6) 70.2 (14.9) 167 (9) 588 73 (16) 588 1544 (29) – – 
 MrOS-USA USA 3925 73.9 (5.9) 83.1 (12.7) 175 (7) 3925 79 (17) 3925 1551 (30) – – 
 SOF USA 2103 100 80.1 (4.2) 66.3 (12.5) 158 (6) 2103 59 (17) 2103 1527 (30) – – 
 YFS Finland 1265 58 37.9 (5.0) 75.8 (15.5) 172 (9) 1265 80 (16) 1265 1559 (29) 1250 0.560 (0.110) 
 HCS-AUS Australia 986 49 66.2 (7.6) 79.4 (15.5) 166 (9) – – – – 986 0.538 (0.166) 
 Subtotal  13 407 52 69.2 (6.9) 75.4 (14.2) 167 (9) 10 466 73 (17) 10 465 1544 (30) 2236 0.550 (0.138) 
De novo replication 
 AUSTRIOS-B Austria 448 85 83.6 (5.9) 62.0 (12.3) 156 (8) 448 90 (17) 448 1496 (36) – – 
 CABRIO-C Spain 1274 62 62.4 (9.2) 73.7 (13.1) 161 (8) 1274 70 (23) 1273 1545 (41) – – 
 CAIFOS Australia 1113 100 80.0 (2.6) 67.5 (12.1) 157 (6) 1113 101 (9) 1113 1516 (28) – – 
 CALEX-FAM Finland 983 79 37.0 (22.4) 64.3 (16.9) 164 (11) 983 83 (16) – – – – 
 EMAS Europeb 2870 59.9 (11.0) 83.1 (13.6) 173 (7) 2870 80 (19) 2870 1550 (34) – – 
 EPICNOR UK 5723 54 63.6 (9.2) 73.2 (12.4) 167 (9) 5723 79 (20) 5718 1638 (43) – – 
 EPOLOS Poland 684 56 53.4 (16.0) 73.2 (13.7) 166 (10) 684 112 (13) 684 1548 (35) – – 
 FLOS Italy 1000 84 59.8 (12.7) 64.8 (12.3) 163 (9) 1000 58 (7) 1000 1503 (83) – – 
 GEOS Canada 5495 100 55.8 (10.3) 65.4 (11.9) 158 (6) 5495 111 (10) 5495 1546 (32) – – 
 LASA Netherlands 894 51 75.6 (6.5) 74.2 (12.6) 166 (9) 894 71 (20) 894 1611 (44) – – 
 MrOS-SWE Sweden 1718 75.4 (3.2) 80.6 (12.0) 175 (7) 1718 81 (21) 1718 1555 (38) – – 
 OPRA Sweden 821 100 75.2 (0.1) 67.6 (11.3) 160 (6) 821 102 (10) 821 1523 (27) – – 
 OSTEOSII Greece 307 87 50.5 (12.6) 74.1 (15.7) 163 (7) 307 112 (16) 307 1556 (36) – – 
 PEAK25 Sweden 857 100 25.5 (0.2) 64.5 (11.2) 168 (6) 857 118 (11) 857 1575 (32) – – 
 SWS UK 715 100 29.7 (3.7) 72.4 (14.8) 163 (7) 714 72 (13) 715 1548 (27) – – 
 Subtotal  24 902 64 60.2 (10.0) 71.6 (12.7) 165 (8) 24 901 89 (16) 23 913 1568 (40) – – 
Total  58 878 62 62.3 (9.7) 73.0 (13.6) 165 (8) 49 627 85 (17) 49 892 1570 (39) 6792 0.502 (0.138) 

aThe AOGC cohort contributed to in silico lookups of SNP-fracture associations only.

bThe EMAS study comprises cohorts in Belgium, Estonia, Hungary, Italy, Poland, Spain, Sweden and UK.

Table 2.

Summary of P-values for association of SNPs in 25 loci with heel BUA, VOS or heel DXA BMD in GWAS discovery/replication meta-analysis

Locus SNP Closest gene Genetic function Discovery P-valuesb
 
Replication P-valuesb
 
Combined P-valuesb
 
BUA VOS DXA BUA VOS DXA BUA VOS DXA 
Combined P < 5 × 10−8  9 cohorts, 14 258 participants 21 cohorts, 35 082 participants 30 cohorts, 49 335 participants 
 2p16.2 rs11898505 SPTBN1 Intronic, regulatory region 7.78×10−8 2.92×10−8 7.68×10−1 6.66×10−12 1.10×10−4 9.63×10−2 4.24×10−13 6.25×10−6 2.65×10−1 
 6q22.33 rs7741021 RSPO3 Intronic, regulatory region 8.52×10−7 1.72×10−7 7.69×10−6 1.19×10−18 2.54×10−21 1.49×10−3 9.26×10−21 9.58×10−20 4.11×10−8 
 6q25.1 rs4869739 CCDC170 Intronic 5.25×10−10 4.75×10−11 7.73×10−10 1.02×10−3 3.92×10−8 3.82×10−1 1.93×10−9 2.64×10−18 1.21×10−2 
 6q25.1 rs3020331c ESR1 Intronic 1.27×10−2 7.94×10−6 2.01×10−4 3.04×10−10 3.79×10−17 1.95×10−1 2.91×10−9 6.64×10−15 1.26×10−3 
 6q25.1 rs2982552 ESR1 Intronic, regulatory region 2.87×10−2 3.31×10−6 3.83×10−4 6.16×10−17 1.14×10−18 1.00×10−1 1.70×10−10 7.32×10−16 1.21×10−4 
 7q31.31 rs2908007 WNT16 Upstream 8.59×10−21 5.02×10−23 4.31×10−11 1.31×10−22 2.06×10−39 3.47×10−2 4.32×10−35 1.62×10−59 1.34×10−9 
 10q21.1 rs7902708 MBL2/DKK1 Intronic 8.23×10−3 1.46×10−7 9.51×10−1 1.02×10−8 6.99×10−9 2.60×10−3 1.30×10−8 5.29×10−15 2.47×10−1 
11q14.2 rs597319 TMEM135 Intronic 2.62×10−4 1.18×10−8 5.05×10−3 2.01×10−12 2.70×10−17 2.20×10−2 8.23×10−14 4.86×10−26 3.05×10−4 
 19q13.11 rs10416265 GPATCH1 Non-synonymous coding 8.30×10−7 2.99×10−8 1.15×10−1 5.84×10−8 2.92×10−5 3.45×10−1 2.37×10−13 4.08×10−12 6.72×10−2 
Combined P ≥ 5 x10−8  9 cohorts, 14 258 participants 21 cohorts, 35 082 participants 30 cohorts, 49 335 participants 
 5p13.3 rs9292469 NPR3 Upstream 3.09×10−6 6.01×10−3 9.27×10−1 5.95×10−1 1.69×10−1 9.96×10−1 1.43×10−1 6.12×10−1 9.42×10−1 
 7p15.2 rs11520772 TAX1BP1 Intronic 9.71×10−7 4.84×10−4 6.24×10−1 8.43×10−2 1.32×10−1 5.48×10−1 2.86×10−4 7.07×10−3 8.79×10−1 
 7p14.1 rs6974574c EPDR1 Upstream 5.81×10−3 1.34×10−5 2.56×10−4 2.51×10−4 4.84×10−3 7.31×10−1 8.25×10−5 3.89×10−5 9.25×10−3 
 7q11.23 rs38664 UPK3B Intronic 9.10×10−4 1.52×10−6 6.60×10−1 4.39×10−2 1.58×10−2 5.35×10−1 3.25×10−4 1.02×10−7 8.79×10−1 
 13q12.3 rs3000634 USPL1 Upstream 2.10×10−5 1.27×10−7 2.18×10−1 6.80×10−3 1.91×10−1 5.38×10−1 8.12×10−1 8.00×10−2 1.70×10−1 
 13q14.11 rs9533090 AKAP11 Upstream 3.78×10−2 5.04×10−3 5.05×10−10 7.60×10−3 2.44×10−4 6.44×10−1 1.02×10−3 1.40×10−5 6.97×10−3 
 16q24.1 rs7188801 FOXL1 Upstream 3.32×10−4 3.09×10−6 2.16×10−2 3.91×10−1 1.66×10−2 5.48×10−1 9.70×10−3 7.62×10−6 2.90×10−2 
    9 cohorts, 14 258 participants 6 cohorts, 10 466 participants 15 cohorts, 24 723 participants 
 2p21 rs17032452 CAMKMT Intronic 8.73×10−1 5.30×10−1 1.74×10−6 5.49×10−1 4.24×10−1 3.59×10−1 6.26×10−1 9.67×10−1 1.56×10−3 
 3p14.2 rs6414591 C3orf67 Upstream 3.49×10−1 2.39×10−1 1.72×10−6 1.31×10−1 9.13×10−2 6.83×10−1 7.86×10−1 8.17×10−1 9.22×10−2 
 5q31.2 rs11959305 TGFBI Intronic 1.89×10−2 1.82×10−2 6.84×10−8 6.52×10−1 2.80×10−1 8.61×10−1 8.47×10−2 7.74×10−3 1.15×10−1 
 7p15.3 rs7787266 STEAP1B Intronic 4.08×10−1 4.93×10−1 2.53×10−6 2.93×10−1 3.14×10−1 6.21×10−1 1.97×10−1 2.70×10−1 9.71×10−3 
 9q21.33 rs10868487 GAS1 Downstream 6.10×10−1 3.81×10−1 2.37×10−6 2.61×10−1 2.50×10−1 7.03×10−1 8.57×10−1 7.46×10−1 8.92×10−2 
 13q31.1 rs9574655 SPRY2 Downstream 2.58×10−1 1.38×10−1 9.09×10−8 8.59×10−1 6.43×10−1 8.67×10−2 5.81×10−1 4.72×10−1 3.50×10−1 
 16q12.2 rs923220 IRX5 Upstream 1.24×10−3 7.98×10−3 6.05×10−7 7.85×10−1 7.28×10−1 9.34×10−1 2.58×10−2 3.95×10−2 1.56×10−2 
 20q11.22 rs3746429 EDEM2 Missense variant 4.42×10−1 8.27×10−1 3.80×10−7 2.07×10−1 9.14×10−2 3.35×10−1 7.23×10−1 3.93×10−1 4.35×10−4 
 21q22.2 rs2836789 FLJ45139 Upstream 1.56×10−1 1.36×10−2 1.51×10−6 2.09×10−3 4.09×10−2 5.07×10−1 3.77×10−3 2.27×10−3 1.57×10−3 
Locus SNP Closest gene Genetic function Discovery P-valuesb
 
Replication P-valuesb
 
Combined P-valuesb
 
BUA VOS DXA BUA VOS DXA BUA VOS DXA 
Combined P < 5 × 10−8  9 cohorts, 14 258 participants 21 cohorts, 35 082 participants 30 cohorts, 49 335 participants 
 2p16.2 rs11898505 SPTBN1 Intronic, regulatory region 7.78×10−8 2.92×10−8 7.68×10−1 6.66×10−12 1.10×10−4 9.63×10−2 4.24×10−13 6.25×10−6 2.65×10−1 
 6q22.33 rs7741021 RSPO3 Intronic, regulatory region 8.52×10−7 1.72×10−7 7.69×10−6 1.19×10−18 2.54×10−21 1.49×10−3 9.26×10−21 9.58×10−20 4.11×10−8 
 6q25.1 rs4869739 CCDC170 Intronic 5.25×10−10 4.75×10−11 7.73×10−10 1.02×10−3 3.92×10−8 3.82×10−1 1.93×10−9 2.64×10−18 1.21×10−2 
 6q25.1 rs3020331c ESR1 Intronic 1.27×10−2 7.94×10−6 2.01×10−4 3.04×10−10 3.79×10−17 1.95×10−1 2.91×10−9 6.64×10−15 1.26×10−3 
 6q25.1 rs2982552 ESR1 Intronic, regulatory region 2.87×10−2 3.31×10−6 3.83×10−4 6.16×10−17 1.14×10−18 1.00×10−1 1.70×10−10 7.32×10−16 1.21×10−4 
 7q31.31 rs2908007 WNT16 Upstream 8.59×10−21 5.02×10−23 4.31×10−11 1.31×10−22 2.06×10−39 3.47×10−2 4.32×10−35 1.62×10−59 1.34×10−9 
 10q21.1 rs7902708 MBL2/DKK1 Intronic 8.23×10−3 1.46×10−7 9.51×10−1 1.02×10−8 6.99×10−9 2.60×10−3 1.30×10−8 5.29×10−15 2.47×10−1 
11q14.2 rs597319 TMEM135 Intronic 2.62×10−4 1.18×10−8 5.05×10−3 2.01×10−12 2.70×10−17 2.20×10−2 8.23×10−14 4.86×10−26 3.05×10−4 
 19q13.11 rs10416265 GPATCH1 Non-synonymous coding 8.30×10−7 2.99×10−8 1.15×10−1 5.84×10−8 2.92×10−5 3.45×10−1 2.37×10−13 4.08×10−12 6.72×10−2 
Combined P ≥ 5 x10−8  9 cohorts, 14 258 participants 21 cohorts, 35 082 participants 30 cohorts, 49 335 participants 
 5p13.3 rs9292469 NPR3 Upstream 3.09×10−6 6.01×10−3 9.27×10−1 5.95×10−1 1.69×10−1 9.96×10−1 1.43×10−1 6.12×10−1 9.42×10−1 
 7p15.2 rs11520772 TAX1BP1 Intronic 9.71×10−7 4.84×10−4 6.24×10−1 8.43×10−2 1.32×10−1 5.48×10−1 2.86×10−4 7.07×10−3 8.79×10−1 
 7p14.1 rs6974574c EPDR1 Upstream 5.81×10−3 1.34×10−5 2.56×10−4 2.51×10−4 4.84×10−3 7.31×10−1 8.25×10−5 3.89×10−5 9.25×10−3 
 7q11.23 rs38664 UPK3B Intronic 9.10×10−4 1.52×10−6 6.60×10−1 4.39×10−2 1.58×10−2 5.35×10−1 3.25×10−4 1.02×10−7 8.79×10−1 
 13q12.3 rs3000634 USPL1 Upstream 2.10×10−5 1.27×10−7 2.18×10−1 6.80×10−3 1.91×10−1 5.38×10−1 8.12×10−1 8.00×10−2 1.70×10−1 
 13q14.11 rs9533090 AKAP11 Upstream 3.78×10−2 5.04×10−3 5.05×10−10 7.60×10−3 2.44×10−4 6.44×10−1 1.02×10−3 1.40×10−5 6.97×10−3 
 16q24.1 rs7188801 FOXL1 Upstream 3.32×10−4 3.09×10−6 2.16×10−2 3.91×10−1 1.66×10−2 5.48×10−1 9.70×10−3 7.62×10−6 2.90×10−2 
    9 cohorts, 14 258 participants 6 cohorts, 10 466 participants 15 cohorts, 24 723 participants 
 2p21 rs17032452 CAMKMT Intronic 8.73×10−1 5.30×10−1 1.74×10−6 5.49×10−1 4.24×10−1 3.59×10−1 6.26×10−1 9.67×10−1 1.56×10−3 
 3p14.2 rs6414591 C3orf67 Upstream 3.49×10−1 2.39×10−1 1.72×10−6 1.31×10−1 9.13×10−2 6.83×10−1 7.86×10−1 8.17×10−1 9.22×10−2 
 5q31.2 rs11959305 TGFBI Intronic 1.89×10−2 1.82×10−2 6.84×10−8 6.52×10−1 2.80×10−1 8.61×10−1 8.47×10−2 7.74×10−3 1.15×10−1 
 7p15.3 rs7787266 STEAP1B Intronic 4.08×10−1 4.93×10−1 2.53×10−6 2.93×10−1 3.14×10−1 6.21×10−1 1.97×10−1 2.70×10−1 9.71×10−3 
 9q21.33 rs10868487 GAS1 Downstream 6.10×10−1 3.81×10−1 2.37×10−6 2.61×10−1 2.50×10−1 7.03×10−1 8.57×10−1 7.46×10−1 8.92×10−2 
 13q31.1 rs9574655 SPRY2 Downstream 2.58×10−1 1.38×10−1 9.09×10−8 8.59×10−1 6.43×10−1 8.67×10−2 5.81×10−1 4.72×10−1 3.50×10−1 
 16q12.2 rs923220 IRX5 Upstream 1.24×10−3 7.98×10−3 6.05×10−7 7.85×10−1 7.28×10−1 9.34×10−1 2.58×10−2 3.95×10−2 1.56×10−2 
 20q11.22 rs3746429 EDEM2 Missense variant 4.42×10−1 8.27×10−1 3.80×10−7 2.07×10−1 9.14×10−2 3.35×10−1 7.23×10−1 3.93×10−1 4.35×10−4 
 21q22.2 rs2836789 FLJ45139 Upstream 1.56×10−1 1.36×10−2 1.51×10−6 2.09×10−3 4.09×10−2 5.07×10−1 3.77×10−3 2.27×10−3 1.57×10−3 

P-values smaller than the genome-wide significance threshold (P < 5 × 10−8) or suggestive significance threshold (P < 5 × 10−6) are indicated in bold typefacea.

The association statistics for a new locus at chr 11q14.2 are italicized.

aThe P-values in the GWAS discovery are based on a fixed-effect meta-analysis model, while those in the replication and combined analyses are based on a random-effects meta-analysis model.

bThe number of cohorts and participants contributing to the analysis of each SNP at each stage slightly varied depending on quality control filters as well as successful imputation or de novo genotyping of the particular SNP. Figure 1 and Supplementary Material, Figure S3 show the exact numbers that were available for each SNP at each stage for the confirmed loci.

cSecondary signals at the discovery phase following conditional analyses within the region (see Supplementary Material, Fig. S6 for the regional association plots).

Figure 1.

Flow chart summarizing key features of the discovery and replication phases.

Figure 1.

Flow chart summarizing key features of the discovery and replication phases.

Associations between the 15 SNPs that were considered for replication primarily on the basis of their association with heel BUA or VOS are shown in Figure 2. The SNP characteristics are summarized in Table 2. In the combined meta-analysis of the discovery and replication cohorts using a random-effects model, 9 SNPs showed genome-wide significant associations, of which 7 were previously reported to be associated with central DXA BMD (19). Two of the SNPs (rs7741021 and rs2908007) also showed genome-wide significant association with heel DXA BMD (Table 2). Three SNPs on chromosome 6q25.1 (rs4869739, rs3020331 and rs2982552) mapped to intronic or regulatory regions around the ESR1 (estrogen receptor 1) and CCDC170 (coiled-coil domain containing 170, previously known as C6orf97) genes (Fig. 3), and five other SNPs mapped to loci within or near previously identified osteoporosis susceptibility genes, including 2p16.2 (SPTBN1, rs11898505), 6q22.33 (RSPO3, rs7741021), 7q31.1 (WNT16, rs2908007), 10q21.1 (DKK1, rs7902708) and 19q13.11 (GPATCH1, rs10416265). We identified a new locus on chromosome 11q14.2 (TMEM135, rs597319) significantly associated with both BUA and VOS (P < 8.23 × 10−14).

Figure 2.

Summary of SNP associations with heel BUA or VOS in GWAS discovery meta-analysis and replication in independent samples of participants. The pooled estimates in the GWAS discovery are based on a fixed-effect meta-analysis model, while those in the replication and combined analyses are based on a random-effects meta-analysis model. Allele b indicates the effect allele, and the presence of two alleles in this column indicates that a proxy SNP with r2 > 0.8 (except for 16q24.1 locus for which r2 = 0.6) was used for the replication analyses.

Figure 2.

Summary of SNP associations with heel BUA or VOS in GWAS discovery meta-analysis and replication in independent samples of participants. The pooled estimates in the GWAS discovery are based on a fixed-effect meta-analysis model, while those in the replication and combined analyses are based on a random-effects meta-analysis model. Allele b indicates the effect allele, and the presence of two alleles in this column indicates that a proxy SNP with r2 > 0.8 (except for 16q24.1 locus for which r2 = 0.6) was used for the replication analyses.

Figure 3.

Association of SNPs at chromosome 6q25.1 region with heel BUA, VOS, and heel DXA BMD in meta-analysis of discovery cohorts before (left column) and after (right column) adjusting for the most significant SNP in the region (i.e. unconditional and conditional analyses respectively); as well as the unconditional results for a novel locus for heel bone properties at chromosome 11q14.2. (The conditional analyses led to the identification of the highlighted secondary signal for association of 6q25.1 with VOS. Conditional analyses results for heel DXA BMD were not available from the three relevant discovery cohorts. Color versions of the above figures have been made available in Supplementary Material, Fig. S6.).

Figure 3.

Association of SNPs at chromosome 6q25.1 region with heel BUA, VOS, and heel DXA BMD in meta-analysis of discovery cohorts before (left column) and after (right column) adjusting for the most significant SNP in the region (i.e. unconditional and conditional analyses respectively); as well as the unconditional results for a novel locus for heel bone properties at chromosome 11q14.2. (The conditional analyses led to the identification of the highlighted secondary signal for association of 6q25.1 with VOS. Conditional analyses results for heel DXA BMD were not available from the three relevant discovery cohorts. Color versions of the above figures have been made available in Supplementary Material, Fig. S6.).

Subsidiary comparisons with fixed-effect meta-analysis results (Supplementary Material, Table S2 and Figs S2 and S3) suggested two additional genome-wide significant loci; one at 7p14.1 upstream of EPDR1 (rs6974574, P < 4.92 × 10−8 for BUA and VOS) and the other at 13q14.11 upstream of AKAP11 (rs9533090, P = 5.33 × 10−8 for VOS), although there was statistically significant between-study heterogeneity in these two loci for the respective phenotypes (Supplementary Material, Table S3), necessitating some caution in generalizing the fixed-effect meta-analysis results. Figure 4 provides a comparison of the magnitudes of association of the 25 SNPs with heel bone measures and central DXA BMD, suggesting generally concordant associations in the overlapping genome-wide significant or suggestive loci.

Figure 4.

Comparison of magnitudes of associations of 25 SNPs with heel bone properties and central DXA BMD. The SNP associations with central DXA BMD are based on lookup of previously published results from GEFOS.

Figure 4.

Comparison of magnitudes of associations of 25 SNPs with heel bone properties and central DXA BMD. The SNP associations with central DXA BMD are based on lookup of previously published results from GEFOS.

We further tested if the genome-wide significant or suggestive genetic loci were associated with fracture risk based on data available from 25 cohorts with up to 54 245 participants, among whom there were 14 958 cases of any fracture (excluding fractures of the skull and extremities, i.e. fingers and toes), 10 663 non-vertebral fractures and 3220 clinical vertebral fractures (Supplementary Material, Table S4). Ten of 10 SNPs associated with heel bone properties at P < 5 × 10−6 showed the expected directions of association with any fracture outcome based on the point estimates (Fig. 5). Furthermore, 6 of these 10 SNPs showed nominally significant (P < 0.05) associations with fractures, including three SNPs with P < 0.005 (i.e. corrected for multiple comparisons using Bonferroni method) at 6q22.33 (rs7741021), 7q31.31 (rs2908007), and 10q21.1 (rs7902708). Fixed-effect meta-analysis gave similar results (Supplementary Material, Fig. S4).

Figure 5.

Per-allele odds ratios for association with fracture risk for 10 SNPs that were associated with heel BUA, VOS or heel DXA BMD at P < 5 × 10−6 in combined meta-analyses using a random-effects model. The pooled estimates are based on a random-effects meta-analysis model. FXANY = any fracture; FXNONVERT = non-vertebral fracture; FXVERT = vertebral fracture. Allele b indicates the effect of allele, and the presence of two alleles in this column indicates that a proxy SNP with r2 > 0.8 was used for the replication analyses.

Figure 5.

Per-allele odds ratios for association with fracture risk for 10 SNPs that were associated with heel BUA, VOS or heel DXA BMD at P < 5 × 10−6 in combined meta-analyses using a random-effects model. The pooled estimates are based on a random-effects meta-analysis model. FXANY = any fracture; FXNONVERT = non-vertebral fracture; FXVERT = vertebral fracture. Allele b indicates the effect of allele, and the presence of two alleles in this column indicates that a proxy SNP with r2 > 0.8 was used for the replication analyses.

Supplementary Material, Figure S5 presents forest plots of the study-specific results and summary estimates by random-effects meta-analysis for the 15 SNPs that were considered for replication primarily on the basis of their association with heel BUA or VOS in GWA discovery meta-analysis, suggesting generally consistent results across cohorts for a majority of the SNPs. Supplementary Material, Figure S6 shows the regional association plots within a one megabase window of the top SNP in each locus in the GWA discovery meta-analysis, demonstrating strong credible association signals for a number of SNPs underlying the loci selected for replication.

Subsidiary investigation of potential sex differences in the association of SNPs and heel BUA or VOS measures did not reveal convincing evidence of potentially important differences, considering the secondary nature of the hypothesis and multiple comparisons done (Supplementary Material, Fig. S7).

DISCUSSION

This is the first large-scale collaborative GWA study for heel bone properties assessed by quantitative ultrasound and DXA of the heel. Its conception was inspired by the observational evidence of association of heel QUS measures and fracture risk (12), independent of central DXA BMD measures (20), demonstration of a reasonably high genetic heritability of heel QUS measures (7), and suggestions of pleiotropic effects of genes in the determination of bone phenotypes (8). Indeed, consistent with the expected similarities and differences in the physical properties of bone determined by DXA and QUS and prior evidence of moderate genetic correlations between the measures (7–9), we found evidence for some genetic loci common to heel QUS measures and central DXA BMD as well as a novel locus for heel QUS at 11q14.2 (TMEM135, rs597319) that had not been previously identified as associated with BMD or other bone phenotype.

Seven of nine genome-wide significant loci found in the present study were previously reported to be associated with BMD of the hip and/or spine (Fig. 4). This complements our previous findings (17–19) and lends support to the hypothesis of partially shared genetic determinants between QUS and BMD measures (7–9). A comparison of the standardized effect sizes (Fig. 4) also revealed existence of some quantitative differences for some SNPs. For example, in the 7q31.31 locus (WNT16), the effect of rs2908007 on heel measures was about three times as great as its effect on hip or spine BMD, supporting Karasik et al.'s finding that there is significant pleiotropy in the effects of genes on bone phenotypes at different measurement sites (8). Similar quantitative differences were also observed for rs7741021 at the 6q22.33 locus (RSPO3). In the absence of bias and assuming minimal type II errors (i.e. adequate power), such quantitative differences in effect sizes of SNPs at different skeletal sites might indicate heterogeneity in genetically mediated responses of the skeleton to environmental stimuli, including for example, ground reaction forces that are particularly high at the heel but are dampened at more proximal sites such as the lumbar spine (21,22).

Perhaps the most intriguing finding was that we identified a new locus for bone phenotypes on chromosome 11q14.2 (rs597319) near the transmembrane protein 135 (TMEM135) gene, that was genome-wide significant for both BUA and VOS. The TMEM135 gene was first identified in a human lung adenocarcinoma cell line cDNA library (23). It has been suggested that it is critically involved in the process of osteoblastogenesis from human multipotent adipose tissue-derived stem cells (24). Marrow fat cells and osteoblasts share a common stromal precursor and there is currently great interest in the role of increased marrow fat in osteoporotic conditions and the metabolic inter-relationships between these neighboring cell types (25). In depth protein sequence analysis showed that TMEM135 is a multi-transmembrane protein with seven transmembrane helices of high confidence. Homologies exist between TMEM135 and the transmembrane region of frizzled-4 (24), a known component of the Wnt signaling pathway (26). ENCODE project (27) data show that two SNPs in the intronic region of TMEM135 and close to our lead signal (rs502580 and rs603140, both with high linkage disequilibrium with rs597319 [r2 >0.92], and both highly associated with QUS outcomes in our discovery cohorts [P ∼ 1.3 × 10–7 for both]) are associated with changes in MIF-1 and Cart1 motifs in osteoblastic cell lines. Interestingly, both of these transcription factors have been previously shown to be associated with skeletal development and bone density (28,29). Furthermore, TMEM135 was previously reported to be associated with longevity in Caenorhabditis elegans models (30) as well as with longevity and walking speed in humans (31). In summary, the associations observed in our study might be the results of direct effects of increased osteoblastogenesis on heel bone properties, or indirect effects mediated through increased mechanical loading of the calcaneum, associated with faster movements.

The other genetic loci with significant associations with heel bone measures have previously been reported to be associated with BMD or fractures. The ESR1 gene has been shown to be related to osteoporosis susceptibility in both candidate gene (32) and GWA studies (18,33). SNPs in SPTBN1 gene were significantly associated with central DXA BMD in a previous meta-analysis of GEFOS cohorts (18), as were SNPs in WNT16, DKK1, and GPATCH1 genes in the recent GEFOS-BMD meta-analysis (19). The RSPO3 gene has recently been suggested as a bone-related locus by a GWA study of extreme low and high BMD populations (34). The spectrin, beta, non-erythrocytic 1 (SPTBN1) gene located at chromosome 2p16.2 codes for the β-subunit of spectrin, which is a molecular scaffold protein essential in linking plasma membrane to the actin cytoskeleton. Spectrin plays an important role in determination of cell shape, positioning of transmembrane proteins, resilience of membranes to mechanical stress, and organization of organelles and molecular traffic in cells. β-Subunits coded by SPTBN1 are responsible for most of the spectrin-binding activity. Despite several GWA studies confirming the association between SPTBN1 and osteoporosis (18,19,33,35), its role in bone pathophysiology is unclear.

The estrogen receptor 1 (ESR1) gene located at chromosome 6q25.1 codes for the estrogen receptor type 1 (also known as ER-α). Two isoforms of estrogen receptors in humans (α and β) are encoded by two different genes (ESR1 and ESR2) and have distinct tissue and cell patterns of expression. Estrogen receptor is a DNA-binding transcription factor that regulates the activity of many different genes. Estrogen is well known to inhibit bone resorption through both direct and indirect actions on osteoclasts, and it is a major anabolic steroid in bone, particularly evident in the establishment of peak bone mass. Postmenopausal bone loss is complex, involving many genetically regulated processes. After menopause, bone is lost rapidly but variably for several years by most women as osteoclastic bone resorptive activity increases in association with osteocyte apoptosis (36). In an osteoporosis GWA study by deCODE Genetics in 2008 (33), several markers close to ESR1 were reported to show association with BMD, including intronic variants and upstream SNPs close to CCDC170 (previously known as C6orf97). This association was replicated in both GEFOS-BMD meta-analyses (18,19), and we found three-independent SNPs in this region associated with heel BUA and VOS. Most recently, this locus has been shown to be more associated with cortical volumetric BMD (as opposed to trabecular BMD), which implies a role of ESR1 products in osteoblastogenesis and cortical porosity (37).

The wingless-type MMTV integration site family, member 16 (WNT16) gene located at chromosome 7q31.31 is part of the Wnt/LRP pathway, which is a known major anabolic pathway in bone (38). The effects of activation of this pathway include differentiation of mesenchymal precursors into osteoblasts, osteoblast proliferation, bone mineralization, and avoidance of osteoblast apoptosis, and inhibition of osteoclastogenesis through effects on expression of OPG and RANKL. Other members of this pathway such as LRP5, LRP4, SOST, WLS, DKK1 and CTNNB1 have previously been associated with BMD at genome-wide significance level (18,19,33,35).

The variant rs7902708 on chromosome 10q21.1 locates between the MBL2 and DKK1 genes and is in close linkage disequilibrium with another SNP in this locus (rs1373004, R2 = 0.87 in HapMap CEU population) that was previously found to have a significant association with BMD and fracture risk in GWA meta-analyses (19). Since the MBL2 (mannose-binding lectin 2) gene product is active in the innate immune system, it is more likely that these variants have a cis regulatory effects on Dickkopf-1 (DKK1), which is a known Wnt signaling pathway inhibitor (39). Several functional studies have showed the role of DKK1 in osteolytic bone lesions in patients with advanced multiple myeloma (40) and its inverse relationship with bone mass has been shown in knockout mouse models (41). A similar relationship to the Wnt signaling pathway has also been proposed for the RSPO3 gene (21). Although GPATCH1 was also found to be associated with hip and spine BMD in a previous GEFOS meta-analysis (19), there is no functional information about it in genomic databases.

Caution must be exercised in interpreting the results of the heel DXA BMD analyses because there were less than 7000 participants contributing to the combined meta-analysis. The obtained results, however, were consistent with the work of Portero et al., suggesting that heel DXA BMD and BUA measure comparable properties of the calcaneum, which reflect the amount of bone mineral in the field of view of the detector (2).

While the current study had limited statistical power in the meta-analysis of SNP associations with fracture outcomes, it was nevertheless encouraging to observe nominally statistically significant and expected directions of associations with fractures for six SNPs associated with heel bone measures, including three SNPs at 6q22.33 (rs7741021), 7q31.31 (rs2908007) and 10q21.1 (rs7902708) whose P-values for association surpassed the multiple testing chance-corrected threshold of P < 0.005. The concordant findings may, albeit indirectly, suggest that some of the genetic susceptibility to fracture could partly be mediated through bone properties (e.g. structural or material) captured by QUS or DXA measures; but larger well-powered studies are needed to appropriately assess such relevance.

In conclusion, the present GWA study reveals the effect of several genes common to central DXA-derived BMD and heel ultrasound/pDXA measures and points to a new genetic locus with potential implications for better understanding of osteoporosis pathophysiology. Quantitative differences seen in the standardized effect sizes of some SNPs at different skeletal sites are potentially indicative of heterogeneity in genetically mediated responses of the skeleton to environmental stimuli, including ground reaction forces that are particularly high at the heel than at central sites.

MATERIALS AND METHODS

Study subjects and measurements

The GEFOS consortium is an international collaboration of investigators dedicated to identify the genetic determinants of osteoporosis (http://www.gefos.org/). In particular, the GEFOS consortium extended the breadth of its predecessor, the Genetic Markers for Osteoporosis (GENOMOS) consortium, into meta-analysis of GWA discovery studies. In the current GEFOS/GENOMOS project, we performed GWA discovery and replication of genetic loci associated with heel bone properties, including QUS (measures: BUA and VOS) and DXA (measure: heel BMD).

The discovery phase comprises 13 cohort studies with GWA data and relevant heel bone phenotypes (including BUA in 14 260 participants from 9 cohorts; VOS in 15 514 participants from 9 cohorts; and heel DXA BMD in 4556 participants from 3 cohorts) arising from populations across North America, Europe and East Asia. Independent replication was performed using summary results from seven cohorts with GWA data (in silico n = 11 452) and analysis of individual-level data from 15 other cohorts in the GENOMOS consortium that were centrally genotyped for candidate polymorphisms by the Kbioscience laboratory in the UK (de novo n = 24 902). Characteristics of the study cohorts/participants are summarized in Table 1. All studies were approved by institutional ethics review committees at the relevant organizations and all participants provided written informed consent. Further descriptive information about the participating cohorts is available from the GEFOS/GENOMOS websites (http://www.gefos.org/?q=studies and http://www.genomos.eu/index.php?page=cohorts).

Genotyping and imputation methods

All the discovery cohorts were genotyped using commercially available Affymetrix (Affymetrix Inc., Santa Clara, CA, USA) or Illumina (Illumina Inc., San Diego, CA, USA) genotyping arrays. Quality control was performed independently for each study according to standard manufacturer protocols and within study procedures. To facilitate meta-analysis, each group performed genotype imputation with IMPUTE or MACH software using genotypes from the HapMap Phase II release 22, NCBI build 36 (CEU or CHB/JPT as appropriate) as reference panels. Each imputation software estimates an overall imputation quality score for each SNP. These quality scores and minor allele frequencies for up to ∼2.5 million SNPs available from each cohort were considered in the meta-analysis.

Association analyses

In the discovery phase, each cohort conducted analyses according to a standard prespecified analysis plan under an additive (i.e. per allele) genetic model. Phenotypes for the association analyses were defined as the sex-specific standardized residuals from linear regression of each outcome variable (BUA, VOS or heel BMD) on age, age-squared, weight, height and machine type (if more than one machine was used). The assumption of normality of residuals in the linear regression model was checked within each cohort for each phenotype and no deviations were reported. The SNP–phenotype associations in each study were adjusted for potential confounding by population substructure using principal components as appropriate; pedigree and twin-based studies—additionally—corrected for family structure. The final results submitted to the Coordinating Center for meta-analysis were the per-allele regression coefficients with corresponding standard errors and P-values for the associations of up to 2.5 million SNPs and standardized residuals of each outcome variable. Analysis of imputed genotypes used either the dosage information from MACH or the genotype probabilities from IMPUTE. The replication analyses used the same analytical procedures as above where applicable (e.g. using study-specific standardized residuals as outcomes).

Meta-analysis

Meta-analysis of the GWA discovery summary results was conducted in two-independent collaborating centers (Cambridge, UK and Boston, USA). Because of potentially limited power to detect sex-specific associations, we prespecified the primary analyses to involve meta-analysis of the pooled data (i.e. males and females combined). Quality control filters applied for exclusions of SNPs from the meta-analysis were: imputation quality score of <0.3 for MACH and <0.4 for IMPUTE, average minor allele frequency of <1% across studies, and SNPs missing from >50% of the cohorts contributing to each outcome. Inverse-variance fixed-effects meta-analysis (using METAL software) was conducted in the discovery set with double genomic correction (42) to control for potential inflation of the test statistics in individual studies and in the meta-analysis. The genome-wide level of statistical significance was set at P < 5 × 10−8 and suggestive level of significance at 5 × 10−8P < 5 × 10−6. There were no extreme genomic inflation factors noted in the discovery phase studies or in the GWA meta-analysis (Supplementary Material, Table S1). QQ plots for the combined GWAS meta-analysis results are provided in Supplementary Material, Figure S1.

To help refine the choice of SNPs to be taken forward for replication, conditional analyses were conducted within a 1 megabase window of the best-associated SNP in each locus in the discovery cohorts, if there was more than one SNP with a suggestive level of significance. These secondary analyses took the SNP in the locus with the lowest P-value and conditioned the analysis of all of the other SNPs in the locus by including it in the regression models. In addition, for loci containing SNPs previously associated with hip or spine BMD in GEFOS (19), we performed additional conditioning on the nearby “BMD SNP”.

The DerSimonian and Laird random-effects model was used for meta-analysis of studies in the replication set and also in the final combined analysis of the discovery and replication studies (43). For each SNP included in the replication phase, we meta-analyzed its association with all three phenotypes, simply for completeness, but interpreted the findings while taking into account the primary outcome that the SNP was associated with in the discovery phase. Fixed-effect meta-analysis results were used for subsidiary comparison. We also conducted meta-analysis of the associations of SNPs with fracture outcomes, using only SNPs that were associated with BUA, VOS or heel DXA BMD at P < 5 × 10−6 in the combined analyses, to assess their potential relevance to this clinical outcome.

SUPPLEMENTARY MATERIAL

Supplementary Material is available at HMG online.

FUNDING

This research and the Genetic Factors for Osteoporosis (GEFOS) consortium have been funded by the European Commission (HEALTH-F2-2008-201865-GEFOS). Several other sources of funding and people have supported work in the contributing cohorts as acknowledged in Supplementary Material, Table S2.

ACKNOWLEDGEMENTS

We thank all study participants for making this work possible and also acknowledge the contributions of several staff in the participating cohorts as detailed in Supplementary Material, Table S2.

Conflict of Interest statement. None declared.

REFERENCES

1
Langton
C.M.
The 25th anniversary of BUA for the assessment of osteoporosis: time for a new paradigm?
Proc. Inst. Mech. Eng. H
 , 
2011
, vol. 
225
 (pg. 
113
-
125
)
2
Portero
N.R.
Arlot
M.E.
Roux
J.P.
Duboeuf
F.
Chavassieux
P.M.
Meunier
P.J.
Evaluation and development of automatic two-dimensional measurements of histomorphometric parameters reflecting trabecular bone connectivity: correlations with dual-energy x-ray absorptiometry and quantitative ultrasound in human calcaneum
Calcif. Tissue Int.
 , 
2005
, vol. 
77
 (pg. 
195
-
204
)
3
Gluer
C.C.
Eastell
R.
Reid
D.M.
Felsenberg
D.
Roux
C.
Barkmann
R.
Timm
W.
Blenk
T.
Armbrecht
G.
Stewart
A.
, et al.  . 
Association of five quantitative ultrasound devices and bone densitometry with osteoporotic vertebral fractures in a population-based sample: the OPUS Study
J. Bone Miner. Res.
 , 
2004
, vol. 
19
 (pg. 
782
-
793
)
4
Krieg
M.A.
Barkmann
R.
Gonnelli
S.
Stewart
A.
Bauer
D.C.
Del Rio
B.L.
Kaufman
J.J.
Lorenc
R.
Miller
P.D.
Olszynski
W.P.
, et al.  . 
Quantitative ultrasound in the management of osteoporosis: the 2007 ISCD Official Positions
J. Clin. Densitom.
 , 
2008
, vol. 
11
 (pg. 
163
-
187
)
5
Nelson
H.D.
Haney
E.M.
Dana
T.
Bougatsos
C.
Chou
R.
Screening for osteoporosis: an update for the U.S. Preventive Services Task Force
Ann. Intern. Med.
 , 
2010
, vol. 
153
 (pg. 
99
-
111
)
6
Gonnelli
S.
Cepollaro
C.
Gennari
L.
Montagnani
A.
Caffarelli
C.
Merlotti
D.
Rossi
S.
Cadirni
A.
Nuti
R.
Quantitative ultrasound and dual-energy X-ray absorptiometry in the prediction of fragility fracture in men
Osteoporos. Int.
 , 
2005
, vol. 
16
 (pg. 
963
-
968
)
7
Howard
G.M.
Nguyen
T.V.
Harris
M.
Kelly
P.J.
Eisman
J.A.
Genetic and environmental contributions to the association between quantitative ultrasound and bone mineral density measurements: a twin study
J. Bone Miner. Res.
 , 
1998
, vol. 
13
 (pg. 
1318
-
1327
)
8
Karasik
D.
Hsu
Y.H.
Zhou
Y.
Cupples
L.A.
Kiel
D.P.
Demissie
S.
Genome-wide pleiotropy of osteoporosis-related phenotypes: the Framingham Study
J. Bone Miner. Res.
 , 
2010
, vol. 
25
 (pg. 
1555
-
1563
)
9
Lee
M.
Czerwinski
S.A.
Choh
A.C.
Demerath
E.W.
Sun
S.S.
Chumlea
W.C.
Towne
B.
Siervogel
R.M.
Unique and common genetic effects between bone mineral density and calcaneal quantitative ultrasound measures: the Fels Longitudinal Study
Osteoporos. Int.
 , 
2006
, vol. 
17
 (pg. 
865
-
871
)
10
Bauer
D.C.
Gluer
C.C.
Cauley
J.A.
Vogt
T.M.
Ensrud
K.E.
Genant
H.K.
Black
D.M.
Broadband ultrasound attenuation predicts fractures strongly and independently of densitometry in older women. A prospective study. Study of Osteoporotic Fractures Research Group
Arch. Intern. Med.
 , 
1997
, vol. 
157
 (pg. 
629
-
634
)
11
Bauer
D.C.
Ewing
S.K.
Cauley
J.A.
Ensrud
K.E.
Cummings
S.R.
Orwoll
E.S.
Quantitative ultrasound predicts hip and non-spine fracture in men: the MrOS study
Osteoporos. Int.
 , 
2007
, vol. 
18
 (pg. 
771
-
777
)
12
Moayyeri
A.
Adams
J.E.
Adler
R.A.
Krieg
M.A.
Hans
D.
Compston
J.
Lewiecki
E.M.
Quantitative ultrasound of the heel and fracture risk assessment: an updated meta-analysis
Osteoporos. Int.
 , 
2012
, vol. 
23
 (pg. 
143
-
153
)
13
Arden
N.K.
Baker
J.
Hogg
C.
Baan
K.
Spector
T.D.
The heritability of bone mineral density, ultrasound of the calcaneus and hip axis length: a study of postmenopausal twins
J. Bone Miner. Res.
 , 
1996
, vol. 
11
 (pg. 
530
-
534
)
14
Karasik
D.
Myers
R.H.
Hannan
M.T.
Gagnon
D.
McLean
R.R.
Cupples
L.A.
Kiel
D.P.
Mapping of quantitative ultrasound of the calcaneus bone to chromosome 1 by genome-wide linkage analysis
Osteoporos. Int.
 , 
2002
, vol. 
13
 (pg. 
796
-
802
)
15
Zhai
G.
Andrew
T.
Kato
B.S.
Blake
G.M.
Spector
T.D.
Genetic and environmental determinants on bone loss in postmenopausal Caucasian women: a 14-year longitudinal twin study
Osteoporos. Int.
 , 
2009
, vol. 
20
 (pg. 
949
-
953
)
16
Michaelsson
K.
Melhus
H.
Ferm
H.
Ahlbom
A.
Pedersen
N.L.
Genetic liability to fractures in the elderly
Arch. Intern. Med.
 , 
2005
, vol. 
165
 (pg. 
1825
-
1830
)
17
Ralston
S.H.
Uitterlinden
A.G.
Genetics of osteoporosis
Endocr. Rev.
 , 
2010
, vol. 
31
 (pg. 
629
-
662
)
18
Rivadeneira
F.
Styrkarsdottir
U.
Estrada
K.
Halldorsson
B.V.
Hsu
Y.H.
Richards
J.B.
Zillikens
M.C.
Kavvoura
F.K.
Amin
N.
Aulchenko
Y.S.
, et al.  . 
Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies
Nat. Genet.
 , 
2009
, vol. 
41
 (pg. 
1199
-
1206
)
19
Estrada
K.
Styrkarsdottir
U.
Evangelou
E.
Hsu
Y.H.
Duncan
E.L.
Ntzani
E.E.
Oei
L.
Albagha
O.M.
Amin
N.
Kemp
J.P.
, et al.  . 
Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture
Nat. Genet.
 , 
2012
, vol. 
44
 (pg. 
491
-
501
)
20
Moayyeri
A.
Kaptoge
S.
Dalzell
N.
Luben
R.N.
Wareham
N.J.
Bingham
S.
Reeve
J.
Khaw
K.T.
The effect of including quantitative heel ultrasound in models for estimation of 10-year absolute risk of fracture
Bone
 , 
2009
, vol. 
45
 (pg. 
180
-
184
)
21
Han
X.H.
Jin
Y.R.
Seto
M.
Yoon
J.K.
A WNT/beta-catenin signaling activator, R-spondin, plays positive regulatory roles during skeletal myogenesis
J. Biol. Chem.
 , 
2011
, vol. 
286
 (pg. 
10649
-
10659
)
22
Zheng
H.F.
Tobias
J.H.
Duncan
E.
Evans
D.M.
Eriksson
J.
Paternoster
L.
Yerges-Armstrong
L.M.
Lehtimaki
T.
Bergstrom
U.
Kahonen
M.
, et al.  . 
WNT16 influences bone mineral density, cortical bone thickness, bone strength, and osteoporotic fracture risk
PLoS Genet.
 , 
2012
, vol. 
8
 pg. 
e1002745
 
23
Liu
F.
Li
Y.
Yu
Y.
Fu
S.
Li
P.
Cloning of novel tumor metastasis-related genes from the highly metastatic human lung adenocarcinoma cell line Anip973
J. Genet. Genomics
 , 
2007
, vol. 
34
 (pg. 
189
-
195
)
24
Scheideler
M.
Elabd
C.
Zaragosi
L.E.
Chiellini
C.
Hackl
H.
Sanchez-Cabo
F.
Yadav
S.
Duszka
K.
Friedl
G.
Papak
C.
, et al.  . 
Comparative transcriptomics of human multipotent stem cells during adipogenesis and osteoblastogenesis
BMC. Genomics.
 , 
2008
, vol. 
9
 pg. 
340
 
25
Fazeli
P.K.
Horowitz
M.C.
MacDougald
O.A.
Scheller
E.L.
Rodeheffer
M.S.
Rosen
C.J.
Klibanski
A.
Marrow fat and bone--new perspectives
J. Clin. Endocrinol. Metab.
 , 
2013
, vol. 
98
 (pg. 
935
-
945
)
26
Robitaille
J.
MacDonald
M.L.
Kaykas
A.
Sheldahl
L.C.
Zeisler
J.
Dube
M.P.
Zhang
L.H.
Singaraja
R.R.
Guernsey
D.L.
Zheng
B.
, et al.  . 
Mutant frizzled-4 disrupts retinal angiogenesis in familial exudative vitreoretinopathy
Nat. Genet.
 , 
2002
, vol. 
32
 (pg. 
326
-
330
)
27
Bernstein
B.E.
Birney
E.
Dunham
I.
Green
E.D.
Gunter
C.
Snyder
M.
An integrated encyclopedia of DNA elements in the human genome
Nature
 , 
2012
, vol. 
489
 (pg. 
57
-
74
)
28
Iioka
T.
Furukawa
K.
Yamaguchi
A.
Shindo
H.
Yamashita
S.
Tsukazaki
T.
P300/CBP acts as a coactivator to cartilage homeoprotein-1 (Cart1), paired-like homeoprotein, through acetylation of the conserved lysine residue adjacent to the homeodomain
J. Bone Miner. Res.
 , 
2003
, vol. 
18
 (pg. 
1419
-
1429
)
29
Onodera
S.
Sasaki
S.
Ohshima
S.
Amizuka
N.
Li
M.
Udagawa
N.
Irie
K.
Nishihira
J.
Koyama
Y.
Shiraishi
A.
, et al.  . 
Transgenic mice overexpressing macrophage migration inhibitory factor (MIF) exhibit high-turnover osteoporosis
J. Bone Miner. Res.
 , 
2006
, vol. 
21
 (pg. 
876
-
885
)
30
Exil
V.J.
Silva
A.D.
Benedetto
A.
Exil
E.A.
Adams
M.R.
Au
C.
Aschner
M.
Stressed-induced TMEM135 protein is part of a conserved genetic network involved in fat storage and longevity regulation in Caenorhabditis elegans
PLoS ONE
 , 
2010
, vol. 
5
 pg. 
e14228
 
31
Lunetta
K.L.
D'Agostino
R.B.
Sr.
Karasik
D.
Benjamin
E.J.
Guo
C.Y.
Govindaraju
R.
Kiel
D.P.
Kelly-Hayes
M.
Massaro
J.M.
Pencina
M.J.
, et al.  . 
Genetic correlates of longevity and selected age-related phenotypes: a genome-wide association study in the Framingham Study
BMC. Med. Genet.
 , 
2007
, vol. 
8
 
Suppl. 1
pg. 
S13
 
32
Li
W.F.
Hou
S.X.
Yu
B.
Li
M.M.
Ferec
C.
Chen
J.M.
Genetics of osteoporosis: accelerating pace in gene identification and validation
Hum. Genet.
 , 
2010
, vol. 
127
 (pg. 
249
-
285
)
33
Styrkarsdottir
U.
Halldorsson
B.V.
Gretarsdottir
S.
Gudbjartsson
D.F.
Walters
G.B.
Ingvarsson
T.
Jonsdottir
T.
Saemundsdottir
J.
Center
J.R.
Nguyen
T.V.
, et al.  . 
Multiple genetic loci for bone mineral density and fractures
N. Engl. J. Med.
 , 
2008
, vol. 
358
 (pg. 
2355
-
2365
)
34
Duncan
E.L.
Danoy
P.
Kemp
J.P.
Leo
P.J.
McCloskey
E.
Nicholson
G.C.
Eastell
R.
Prince
R.L.
Eisman
J.A.
Jones
G.
, et al.  . 
Genome-wide association study using extreme truncate selection identifies novel genes affecting bone mineral density and fracture risk
PLoS Genet.
 , 
2011
, vol. 
7
 pg. 
e1001372
 
35
Styrkarsdottir
U.
Halldorsson
B.V.
Gretarsdottir
S.
Gudbjartsson
D.F.
Walters
G.B.
Ingvarsson
T.
Jonsdottir
T.
Saemundsdottir
J.
Snorradottir
S.
Center
J.R.
, et al.  . 
New sequence variants associated with bone mineral density
Nat. Genet.
 , 
2009
, vol. 
41
 (pg. 
15
-
17
)
36
Tomkinson
A.
Reeve
J.
Shaw
R.W.
Noble
B.S.
The death of osteocytes via apoptosis accompanies estrogen withdrawal in human bone
J. Clin. Endocrinol. Metab.
 , 
1997
, vol. 
82
 (pg. 
3128
-
3135
)
37
Paternoster
L.
Lorentzon
M.
Lehtimaki
T.
Eriksson
J.
Kahonen
M.
Raitakari
O.
Laaksonen
M.
Sievanen
H.
Viikari
J.
Lyytikainen
L.P.
, et al.  . 
Genetic determinants of trabecular and cortical volumetric bone mineral densities and bone microstructure
PLoS Genet.
 , 
2013
, vol. 
9
 pg. 
e1003247
 
38
Kubota
T.
Michigami
T.
Ozono
K.
Wnt signaling in bone metabolism
J. Bone Miner. Metab.
 , 
2009
, vol. 
27
 (pg. 
265
-
271
)
39
Fedi
P.
Bafico
A.
Nieto
S.A.
Burgess
W.H.
Miki
T.
Bottaro
D.P.
Kraus
M.H.
Aaronson
S.A.
Isolation and biochemical characterization of the human Dkk-1 homologue, a novel inhibitor of mammalian Wnt signaling
J. Biol. Chem.
 , 
1999
, vol. 
274
 (pg. 
19465
-
19472
)
40
Kocemba
K.A.
Groen
R.W.
van Andel
H.
Kersten
M.J.
Mahtouk
K.
Spaargaren
M.
Pals
S.T.
Transcriptional silencing of the Wnt-antagonist DKK1 by promoter methylation is associated with enhanced Wnt signaling in advanced multiple myeloma
PLoS ONE
 , 
2012
, vol. 
7
 pg. 
e30359
 
41
MacDonald
B.T.
Joiner
D.M.
Oyserman
S.M.
Sharma
P.
Goldstein
S.A.
He
X.
Hauschka
P.V.
Bone mass is inversely proportional to Dkk1 levels in mice
Bone
 , 
2007
, vol. 
41
 (pg. 
331
-
339
)
42
Willer
C.J.
Li
Y.
Abecasis
G.R.
METAL: fast and efficient meta-analysis of genomewide association scans
Bioinformatics
 , 
2010
, vol. 
26
 (pg. 
2190
-
2191
)
43
DerSimonian
R.
Laird
N.
Meta-analysis in clinical trials
Control Clin. Trials
 , 
1986
, vol. 
7
 (pg. 
177
-
188
)

Author notes

Coordination and writing group.
Deceased.
Main principal investigator from each cohort and coordinators for the current project.
§
GEFOS/GENOMOS Coordinating Center and principal investigator (A.G.U.).

Supplementary data