|Home | About | Journals | Submit | Contact Us | Français|
Previous genome-wide association studies (GWAS) have identified common variants in genes associated with variation in bone mineral density (BMD), although most have been carried out in combined samples of older women and men. Meta-analyses of these results have identified numerous SNPs of modest effect at genome-wide significance levels in genes involved in both bone formation and resorption, as well as other pathways. We performed a meta-analysis restricted to premenopausal white women from four cohorts (n= 4,061 women, ages 20 to 45) to identify genes influencing peak bone mass at the lumbar spine and femoral neck. Following imputation, age- and weight-adjusted BMD values were tested for association with each SNP. Association of a SNP in the WNT16 gene (rs3801387; p=1.7 × 10−9) and multiple SNPs in the ESR1/C6orf97 (rs4870044; p=1.3 × 10−8) achieved genome-wide significance levels for lumbar spine BMD. These SNPs, along with others demonstrating suggestive evidence of association, were then tested for association in seven Replication cohorts that included premenopausal women of European, Hispanic-American, and African-American descent (combined n=5,597 for femoral neck; 4,744 for lumbar spine). When the data from the Discovery and Replication cohorts were analyzed jointly, the evidence was more significant (WNT16 joint p=1.3 × 10−11; ESR1/C6orf97 joint p= 1.4 × 10−10). Multiple independent association signals were observed with spine BMD at the ESR1 region after conditioning on the primary signal. Analyses of femoral neck BMD also supported association with SNPs in WNT16 and ESR1/C6orf97 (p< 1 × 10−5). Our results confirm that several of the genes contributing to BMD variation across a broad age range in both sexes have effects of similar magnitude on BMD of the spine in premenopausal women. These data support the hypothesis that variants in these genes of known skeletal function also affect BMD during the premenopausal period.
Osteoporosis, a major public health problem due to the associated increased risk of fractures and consequent morbidity and excess mortality, is a skeletal disorder characterized by compromised bone strength predisposing to an increased risk of fracture. Low bone mass is a major contributor to reduced bone strength. In the United States osteoporosis accounts for 1.5 million new fracture cases each year, and more than 40 million people (primarily women) either have osteoporosis or are at high risk due to low bone mass. (1)
Previous studies have demonstrated that bone mineral density (BMD) in older adults is influenced by genes underlying both bone accrual (i.e., peak bone mass) and bone loss.(2) The contribution of genetic factors to variation in BMD is substantial;(3,4) however, only recently have specific genes and variants associated with BMD been identified and replicated. (5–12) Because of the small effect sizes of identified variants, the most robust results have been obtained by combining results across numerous studies, most including older women and men, to obtain very large sample sizes. These highly powered meta-analyses with very large discovery and replication samples have both replicated known genes or variants and detected novel genes or variants contributing to BMD variability.(13,14) However, these studies have not focused on the genetic determinants of premenopausal bone density.
Peak BMD, which generally occurs around the third to fourth decade of life at the lumbar spine, and much earlier at the hip,(15) is one of the strongest predictors of BMD in later life. Peak BMD is highly heritable, with estimates in Caucasian samples that 50–85% of its variation is genetically determined.(2) De-coupling the effects of peak BMD and premenopausal bone loss(16,17) on BMD in later life is important for understanding osteoporosis risk. Bone loss with age is influenced by substantial environmental as well as genetic factors. These factors may or may not also contribute to peak BMD. Thus, studies of older populations in whom bone loss is likely to have begun may not identify genes that uniquely contribute to acquisition and maintenance of peak BMD. Rather, in women changes in BMD related to the menopausal transition may be influenced by genetic variants that do not reflect peak BMD and may mask the ability to identify specific variants associated with this phenotype.
This study was designed to identify genes contributing to peak BMD in premenopausal women. While several drugs exist to inhibit bone loss, including bisphosphonates and RANKL inhibitors, there are currently few medicines aimed at enhancing bone formation. Thus, the identification of genes responsible for bone accrual may serve to identify drug targets whose manipulation could aid in the improvement of BMD, which has been a successful strategy for therapeutic development.(18)
Genotypic and phenotypic data were obtained from four cohorts of well-characterized European-American premenopausal women, with DXA BMD data available at both the lumbar spine (L1-L4 or L2-L4) and femoral neck (Table 1 and Supplemental Methods). All four Discovery cohorts also contributed GWAS results to the Estrada et al analysis(13); the Discovery cohorts for our study (except the Indiana Sister sample) represent subsets of larger samples also containing older individuals of both sexes. To obtain a sample of female subjects yet to reach menopause and as near to their peak BMD level as possible, the age range was limited to subjects between 20 and 45 years of age at the time BMD was obtained. When serial BMD measures within this range were available, the earliest measure was used. The precise age of menopause could not be uniformly determined in all participating cohorts; therefore, 45 years was chosen because it is the inflection point in the trend of BMD with age(19) and because population-based prospective data show this is where major bone loss begins.(20) Exclusion criteria for all studies included inability to have BMD measured because of obesity, irregular menses or a history of pregnancy or lactation within three months prior to BMD measurement. Women taking oral contraceptives were not excluded. Additional cohort-specific criteria for exclusion of peri- and post-menopausal women are detailed in the Supplemental Material.
Genotyping was performed independently in each sample and results of the Discovery GWAS studies were previously reported.(7,8,21,22) Each study independently imputed their autosomal data to the CEU HapMap II release 22 SNP map; 2,543,887 genotyped and imputed SNPs were then available for analysis in each cohort. Association analyses of each SNP with the two BMD phenotypes were conducted in each cohort separately. Each study included appropriate covariates such as age, weight, and familial relationships, as needed (see Supplemental Material for details). All cohorts contributing to the meta-analysis had study approval from their appropriate institutional research ethics boards, and written informed consent was provided by all subjects.
Results from the different Discovery Sample datasets were combined using inverse variance fixed-effects meta-analysis as implemented in METAL.(23) The contribution of each study to the meta-analysis was weighted by the standard error of the SNP association parameter (beta of the SNP effect). Double genomic control correction was applied to account for possible inflation of association test statistics due to cryptic stratification among or between the Discovery cohorts.
SNPs with p-values < 5 × 10−5 were prioritized for further review and potential inclusion in the replication studies. All SNPs meeting this threshold were reviewed in detail to determine whether: 1) at least 3 of the 4 studies in the Discovery Sample supported the association (β in same direction); 2) there was at least 1 other SNP within the surrounding 100 kb region that supported the association (p<10−4); and 3) the MAF was >10%, so as to have sufficient power to replicate the association in the Replication Sample. We then pruned the set of SNPs in a chromosomal region to remove all SNPs in linkage disequilibrium (r2>0.3) with the SNP having the most significant p-value in the region. Using this approach, a total of 23 SNPs were selected for genotyping in an independent Replication Sample.
For replication, we identified 7 additional studies that included premenopausal women for whom BMD measurements had been made at the femoral neck and lumbar spine sites. The majority of these samples were independent of those employed in the larger, more heterogeneous BMD study by Estrada(13) (Table 1 and Supplemental Methods). The primary inclusion and exclusion criteria for the Replication Sample were identical to those applied to the Discovery Sample (details of each Replication cohort are provided in the Supplemental Material). Each dataset included in the Replication Sample was analyzed separately, and summary results were then meta-analyzed in the same manner as for the Discovery Sample.
A total of 23 independent SNPs were evaluated for association in the Replication Sample. A joint analysis of the data from the Discovery and Replication Samples was conducted for the 23 SNPs on which data were available for all the contributing studies (termed the Joint Sample). The significance of these joint meta-analysis results was evaluated using a genome-wide threshold of 5 × 10−8. Due to the racial diversity of the Replication Sample, the Joint Analysis was also performed in only those cohorts that were Caucasian and non-Hispanic.
The GWAS meta-analyses revealed two regions (WNT16 and ESR1/C6orf97) for which multiple SNPs provided evidence of association with lumbar spine BMD. We performed conditional analyses to test whether there was evidence to support other independent association signals for BMD in these regions (see Supplemental Materials for details).
The Discovery sample consisted of 4,061 women with BMD measurements. Results from the GWAS meta-analysis of lumbar spine BMD are presented in Figure 1A. Following double genomic control, there was no evidence to suggest any inflation of the association results (genomic inflation factor = 1.001; Supplemental Figure 1A). Details of SNPs with evidence of association with lumbar spine BMD are provided in Table 2.
Two SNPs in the lumbar spine BMD GWAS exceeded conventional criteria for genome-wide significance (p<5 × 10−8). The most significant result was on chromosome 7 with a SNP in the last intron of WNT16 (rs3801387; p-value=1.7×10−9). The association extends to SNPs with similar minor allele frequency (MAF=0.26–0.29, p<10−7) across a linkage disequilibrium block containing WNT16 and the neighboring gene FAM3C (Figure 2A).
The second genome-wide significant result was obtained with a SNP on chromosome 6 in C6orf97 (rs4870044; p-value = 1.3 × 10−8). For the chromosome 6 region including ESR1 and C6orf97, the evidence for association extended to nearby SNPs with a range of MAF; these SNPs included rs6930633 (p=1.0 × 10−6; MAF=0.33) and rs2982575 (p=8.3 × 10−7; MAF=0.47, Figure 2B).
Several other SNPs associated with suggestive p-values of 5 × 10−5 or less were genotyped in the Replication Sample (Table 2). These included 2 additional SNPs in the ESR1/C6orf97 region, and 9 SNPs at genomic locations distinct from WNT16 and ESR1/C6orf97.
Results from the GWAS meta-analysis of femoral neck BMD are presented in Figure 1B. Following double genomic control, there was no evidence of genomic inflation (genomic inflation factor = 1.003; Supplemental Figure 1B). Further details on SNPs with evidence of association are provided in Table 3. No SNPs in the femoral neck BMD GWAS exceeded conventional criteria for genome-wide significance.
Several SNPs had p-values making them eligible for consideration in the panel to be genotyped in the Replication Sample (p<5×10−5; Table 3). This included rs3801387 in WNT16 (p=2.68 × 10−6), which had already met our criteria for inclusion in the replication genotyping based on the results with the lumbar spine BMD phenotype (p=1.7×10−9). Several associations unique to the femoral neck BMD meta-analysis were also identified. These included rs1566045 on chromosome 16, 120 kb 5’ of LOC100652974 (p=2.57 × 10−7) and rs7386 on chromosome 11, 28 kb 5’ of METTL12 (p=6.12 × 10−6).
There were 5,597 women in the Replication Sample, of whom 4,158 were of European descent. Only femoral neck BMD was available for the ALSPAC cohort (Replication n=4,744 for lumbar spine BMD). A total of 23 SNPs in 20 distinct chromosomal regions were selected for genotyping in the Replication Sample. Results for the 23 SNPs are presented for lumbar spine BMD (Table 2) and femoral neck BMD (Table 3).
In the joint meta-analysis of both the Discovery and Replication Samples, the evidence for the lumbar BMD association with SNPs in WNT16 and ESR1/C6orf97 became more significant. In the WNT16 region, the p-value in the Joint Analysis for rs3801387 was 1.3 × 10−11 (Table 2). In the ESR1/C6orf97 region, the p-value for rs4870044 improved to 1.4 × 10−10 and another SNP, rs6930633, now met genome-wide significance (p = 1.16 × 10−8; Table 2). For most other SNPs, the Joint Analysis did not substantially improve the evidence of association observed in the Discovery Sample. Plots of the effect size estimates for each of the cohorts and for the Joint Analysis as a whole are provided for WNT16 SNP rs3801387 (Figure 3A) and for the SNPs in the ESR1/C6orf97 region (Figure 3B).
The Joint Analysis of the femoral neck BMD phenotype (Table 3) did not result in any SNPs achieving genome-wide significance. However, for several SNPs (rs3801387 in WNT16 and rs4870044 in the ESR1/C6orf97 region) the evidence for association improved following the Joint Analysis (p = 4.3 × 10−7 and 1.0 × 10−5, respectively; Table 3).
The Replication Sample included an African American and a Hispanic cohort. Given the potential heterogeneity that could be introduced with their inclusion, we performed the Joint Analysis without these cohorts included. The sample now including 8,219 Caucasian and non-Hispanic subjects (85% of the original) yielded similar association results, although for a few SNPs, the evidence of association was substantially improved (ESR1, rs6930633: p = 8.9 × 10−11 and rs2982575: p=4.6 × 10−8).
Conditional analyses were performed to determine whether the associations detected at the WNT16 or C6orf97 loci might be the result of multiple independent genetic factors in these regions. For WNT16, including the SNP (rs3801387) with the most significant lumbar spine BMD GWAS p-value as a covariate in association analysis resulted in nearly complete reduction of the association signal in this region (compare Figures 2A and and4A),4A), suggesting that there is no evidence for a second, independent genetic factor in this region. In contrast, for the ESR1/C6orf97 region, the evidence of association with the other SNPs in the region was largely unaffected when conditioning on any of the 3 individual SNPs from this region included in the replication set (compare Figure 2B and Figures 4B–D). This supports the hypothesis that there is more than one genetic factor influencing peak BMD in the ESR1/C6orf97 region. Similar results were observed in the Replication Sample and Joint Analysis when the same 3 SNPs were analyzed (Table 4 and Supplemental Figure 2).
This is the first GWAS analysis of an exclusively premenopausal female sample to identify genes associated with peak BMD. The two-stage design included an initial Discovery stage consisting of genome-wide analysis of a sample of premenopausal women of European descent for lumbar spine and femoral neck BMD phenotypes. This was followed by a second-stage analysis of 23 SNPs in a Replication Sample using the same BMD phenotypes. The Replication Sample was more diverse and included women of European, Hispanic-American, and African-American descent. We identified SNP associations in the Discovery Sample that reached genome-wide significance. The joint analysis of the Discovery and Replication Samples resulted in greater evidence of association for three SNPs in two distinct regions (WNT16 and ESR1/C6orf97). For two SNPs in ESR1, the evidence of association improved substantially when the analysis was limited to only the Caucasian, non-Hispanic cohorts.
We found the greatest evidence for association in this study between lumbar spine BMD and the SNP rs3801387, in the last intron of the WNT16 gene. The evidence of association extended to SNPs in linkage disequilibrium with rs3801387 in the neighboring gene, FAM3C (Figure 2A), which encodes a less well-characterized interleukin-like molecule. Suggestive evidence of association to this region was also observed in analyses of femoral neck BMD, suggesting an effect on BMD at multiple skeletal sites.
WNT16 is a member of the Wnt-frizzled signaling pathway, and is believed to be a ligand for transmembrane receptors of the frizzled gene family.(13). Functionally, Wnt16 is involved in specification of the sclerotomal somite compartment, which houses vertebral and vascular smooth muscle cell precursors and has been proposed to signal via the non-canonical pathway(24), which is conserved in vertebrates. Wnt16 knockout mice have reduced bone strength and increased susceptibility to fracture(25). SNPs in the human WNT16/FAM3C region, and rs3801387 in particular, have been associated with BMD and bone phenotypes at other skeletal sites, including total body and skull BMD(26), as well as radial BMD, cortical thickness of the tibia as measured by pQCT and forearm fracture.(25) SNPs in WNT16 have also been associated with bone phenotypes across the human lifespan. Previous studies have shown an association with BMD in early childhood(26) and in heterogeneous samples primarily composed of elderly individuals of European descent.(13) The effect size observed for rs3801387 association in our Joint Analysis with lumbar spine BMD is larger than that observed for lumbar spine BMD in the large sample of primarily older women from the GEFOS consortium reported by Estrada et al(13). A similar trend is observed for estimated effect sizes of top SNPs from our study with those observed in the larger, more heterogeneous sample (Table 5). This suggests that earlier periods of life might be important with respect to which genetic variants in WNT16 and other BMD genes identified by the GWAS approach have a lasting effect.
Two recent studies(25,26) have detected multiple independent genetic factors at WNT16 contributing to BMD at skeletal sites not examined in our study, total body and radial BMD. In our conditional analyses, we did not detect evidence to support a second independent genetic factor in that region influencing peak spine BMD. However, our study was not sufficiently powered to exclude the presence of a second signal in this region.
We also detected significant evidence of association with lumbar spine BMD and SNPs in the chromosomal region encompassing the ESR1, estrogen receptor alpha gene. ESR1 has multiple isoforms and a distinct, non-overlapping open reading frame near its 5’ terminus, named C6orf97, whose biological role is poorly understood. The SNPs associated with lumbar spine BMD are located within the C6orf97 sequence or the non-coding region 5’ of the longer ESR1 isoforms. Other investigators have reported complex patterns of association with genetic variants at or near the ESR1 locus with bone-related phenotypes such as radial BMD and heel ultrasound parameters.(27) ESR1 encodes a nuclear transcription factor, with evidence suggesting a possible skeletal role in periosteal mechanosensitivity and expansion.(28,29)
Because the chromosome 6 SNPs associated with lumbar spine BMD had differing MAF (0.29, 0.33 and 0.49) and low levels of linkage disequilibrium between them, we performed conditional association analysis to test whether these SNPs represented independent genetic factors influencing peak BMD. Association results were unchanged when either of the two nearby ESR1/C6orf97 SNPs was included in the association model along with the SNP under consideration. This confirms the finding of the prior study on radius and heel BMD phenotypes(27) that there are multiple independent factors in this chromosomal region. These results are also consistent with the results from Estrada et al,(13) who detected a second association achieving genome-wide significance in the ESR1 region in a more heterogeneous sample of men and women with a wider age range. Investigators of common-variant associations with complex disease should continue to pursue detailed analyses of key findings using conditional analyses and other approaches to enable the most complete dissection possible of the role of common variants in complex diseases.
The association results were less significant for femoral neck than for lumbar spine BMD. This may be due to several factors. First, our analyses tended to identify smaller effect sizes for femoral neck BMD as compared with lumbar spine BMD for the same SNPs (Table 5), suggesting we have lower power to detect associations with femoral neck BMD. Second, there are higher measurement errors at the femoral neck skeletal site. Third, the age range of subjects in our analyses were further beyond the typical age of peak bone mass at the femoral neck than at the lumbar spine; peak lumbar spine bone mass occurs at 33- 40 years in women, while femoral neck peak bone mass occurs earlier (16–19 years).(15) This suggests that even though we included only premenopausal women, the femoral neck BMD may be influenced by both bone acquisition and bone loss. Fourth, the spine and femoral neck sites could be differentially affected by lifestyle factors such as weight-bearing or high-impact exercise. However, suggestive p-values for several SNPs were observed at both sites, notably for SNPs in the WNT16 and ESR1 regions, suggesting effects of variants at these loci across both BMD phenotypes.
We also compared our results with those obtained in the largest BMD GWAS meta-analysis published to date from the GEFOS consortium(13), which analyzed men and women across a wide age range. The two genes in which we identified SNPs with genome-wide evidence of association, WNT16 and C6orf97/ESR1, also contained SNPs associated with BMD in the GEFOS study. As shown in Table 5, we generally did not find strong evidence that the same SNPs in these genes were associated with BMD in the current study as in the GEFOS paper(13) ; and only 4 of the 56 loci detected in the joint female GEFOS analysis(13) were observed to have p-values below 5 × 10−5 in the present study for one or both BMD phenotypes (Supplemental Table 1). This may be due to several factors. One might be that some genes contribute to BMD variation primarily in premenopausal women. Another factor could be that because we limited our sample to premenopausal women, we had limited power to detect genetic factors of small effect. The Discovery Sample used in our study had 67% power at our inclusion threshold to detect associations of similar effect (standardized regression coefficient = 0.10) to the larger effects observed in the much larger meta-analysis of Estrada et al(13) for common SNPs (MAF=0.30–0.40). Thus, we were reasonably powered to detect SNP associations that were at least as large as those identified in the phenotypically more heterogeneous GEFOS study, but may not have been powered to detect these effects if they were smaller in our sample.
In conclusion, we identified SNPs in two genes that were significantly associated with BMD in premenopausal women, who are likely to be at or around peak BMD. These genes have been previously reported to affect BMD at the same skeletal sites in heterogeneous samples of older and younger subjects. Importantly, there is evidence in premenopausal women to support allelic heterogeneity in ESR1. Future studies will focus on the identification of the genetic variants that underlie these replicated associations.
Study design: DLK, DK, MJE, BDM, JBR, DPK, TF
Study conduct: SG, SAC, ACC, FR, JHT, KA, MJE, BDM, JBR, DPK
Data collection: FM, SG, MP, SAC, NJT, DAL, BT, TDS, FR, JCP, MAC, JHT, KA, MJE, BDM, JBR, DPK
Data analysis: DLK, HFZ, LYA, CTL, FM, JPK, SG, DL, SAC, JB
Data interpretation: All
Drafting manuscript: DLK, TF
Revising manuscript content: All
Approving final version of manuscript: All
DLK and TF take responsibility for the integrity of the data analysis.
These studies were supported by R01AG041517, M01 RR-00750, R01AR/AG41398, R01R050066, N01-HC-25195, N02-HL-6-4278, R01AR043351, P01HL045522, R01MH078111, R01MH083824, R01HD012252, R01AR052147, R01DK064391, R01AR43351, R01AG18728, R01HL088119, R01AR046838, U01 HL084756, P30DK072488, and F32AR059469. Genotyping was provided by CIDR (HHSN268200782096C). Support was also provided by the Canadian Institutes of Health Research (CIHR); Merck Frosst Canada Ltd.; Eli Lilly Canada Inc.; Novartis Pharmaceuticals Inc.; The Alliance: Sanofi-Aventis & Procter and Gamble Pharmaceuticals Canada Inc.; Servier Canada Inc.; Amgen Canada Inc.; The Dairy Farmers of Canada; and The Arthritis Society. Drs. Richards and Zheng are supported by the CIHR, the Lady Davis Institute, Fonds de Rechercheen Santé du Québec, Ministere Développement Économique, Innovation et Exportation du Québec, the UK Medical Research Council (092731) and Wellcome Trust (092731; WT088806), British Heart Foundation (SP/07/008/24066), European Commission Framework (FP7/2007–2013), ENGAGE project HEALTH-F4-2007-201413, FP5 GenomEUtwin Project (QLG2-CT-2002-01254), Arthritis Research Campaign, Chronic Disease Research Foundation, and a gift from the AT&T Foundation, C06 RR017515. The National Institute for Health Research (NIHR) comprehensive Biomedical Research Centre award to Guy’s and St Thomas’ NHS Foundation Trust in partnership with King’s College London, and a Biotechnology and Biological Sciences Research Council project grant (G20234).
All authors state that they have no conflicts of interest.