An atlas of genetic influences on osteoporosis in humans and mice

Abstract

Osteoporosis is a common aging-related disease diagnosed primarily using bone mineral density (BMD). We assessed genetic determinants of BMD as estimated by heel quantitative ultrasound in 426,824 individuals, identifying 518 genome-wide significant loci (301 novel), explaining 20% of its variance. We identified 13 bone fracture loci, all associated with estimated BMD (eBMD), in ~1.2 million individuals. We then identified target genes enriched for genes known to influence bone density and strength (maximum odds ratio (OR) = 58, P = 1 × 10^-75) from cell-specific features, including chromatin conformation and accessible chromatin sites. We next performed rapid-throughput skeletal phenotyping of 126 knockout mice with disruptions in predicted target genes and found an increased abnormal skeletal phenotype frequency compared to 526 unselected lines (P < 0.0001). In-depth analysis of one gene, DAAM2, showed a disproportionate decrease in bone strength relative to mineralization. This genetic atlas provides evidence linking associated SNPs to causal genes, offers new insight into osteoporosis pathophysiology, and highlights opportunities for drug development.

Figure 1.. Manhattan plot of genome-wide association results for eBMD in the UK Biobank.

The dashed red line denotes the threshold for declaring genome-wide significance (6.6×10⁻⁹). 1,103 conditionally independent SNPs at 515 loci passed the criteria for genome-wide significance in n=426,824 UK Biobank participants. 301 novel loci (defined as > 1 Mbp from previously reported genome-wide significant BMD variants) reaching genome-wide significance are displayed in blue. Previously reported loci that reached genome-wide significance are displayed in red, and previously reported loci failing to reach genome-wide significance in our study are shown in black.

Figure 2.. Fine-mapping SNPs and target gene selection diagram.

a) For each 500 Mbp region around a conditionally independent lead SNP (p<6.6×10⁻⁹ after conditional independence testing; n=426,824 UK Biobank participants) we applied statistical fine-mapping to calculate log₁₀ Bayes factors for each SNP as a measure of their posterior probability for causality. Conditional independence testing was implemented using GCTA-COJO^, and log₁₀ Bayes factors were estimated using FINEMAP. SNPs that were conditionally independent lead SNPs or that had log₁₀ Bayes factors > 3 were considered our fine-mapped SNPs that we then used for target gene identification. b) Target Genes were identified if: 1) It was the gene closest to a fine-mapped SNP. 2) A fine-mapped SNP was in its gene body. 3) A fine-mapped SNP was coding. 4) The gene mapped closest to a fine-mapped SNP which resided in an SaOS-2 ATAC-seq peak. 5) A fine-mapped SNP was present in a Hi-C osteoblast or osteocyte promoter interaction peak, therefore being closer to a target gene in three-dimensions than linearly on the genome.

Figure 3.. SNPs at genome-wide significant loci are enriched for bone-relevant open chromatin sites.

Comparison of eBMD-associated SNPs in terms of enrichment for DHSs from primary osteoblasts, and ATAC-seq peaks from SaOS-2 osteosarcoma cells. Odds ratios were computed relative to all SNPs at genome-wide significant regions. Enrichments for missense protein coding SNPs are shown as baselines. a) Enrichments for conditionally independent (COJO) or log₁₀ Bayes factor >3 (FINEMAP); note the latter set contains nearly twice the number of SNPs. b) Ranking SNPs by log₁₀ Bayes factor (x-axis) showed increasing enrichment. 95% confidence interval (shaded region) was calculated by a two-sided Fisher’s Exact Test.

Figure 4.

Target Gene Identification Workflow.

Figure 5.. Reduction of DAAM2 protein resulted in reduced mineralization in SaOS-2 cells.

Mineralization quantification in control cells and DAAM2 exon 2 double-stranded break (DSB) induced cells in either the presence of osteogenic factors (treated) or absence (untreated). a) Dot plot of n=6 independent experiments ± standard error of the mean (SEM) from Alizarin red staining in (b) to quantify mineralization; Bar=5mm. ***p=1.3×10⁻¹⁵ compared to untreated control cells and ^&&&p=9.3×10⁻¹⁵ (left) and 8.2×10⁻¹³ (right) compared to treated control cells determined by one-way ANOVA (F=49.7, df=5) and Bonferroni post-hoc tests.

Figure 6.. Biomechanical analyses of mice with Daam2 knockdown.

a) Femur biomechanical analysis. Destructive 3-point bend testing (Instron 5543 load frame) of femurs from wild-type (WT, n_female=3, n_male=4), Daam2^+/tm1a (n_female=6, n_male=4) and Daam2^tm1a/tm1a (n_female=5, n_male=9) mice. Graphs show yield load, maximum load, fracture load, stiffness (gradient of the linear elastic phase) and toughness (energy dissipated prior to fracture). Female data are shown on the left and male data on the right. Data are shown as mean ± standard error of the mean (SEM). Female maximum load analyses for WT versus Daam2^tm1a/tm1a (**) and Daam2^+/tm1a versus Daam2^tm1a/tm1a (^#) had statistically significant differences (one-way ANOVA p=3.0×10⁻³, F=10.29, df=13, Tukey’s post-hoc test **p<0.01 and ^#p<0.05). Male maximum load analyses for WT versus Daam2^tm1a/tm1a (***) and Daam2^+/tm1a versus Daam2^tm1a/tm1a had statistically significant differences [one-way ANOVA p<1.0×10⁻⁴ (GraphPad Prism does not report smaller p-values), F=50.11, df=16, Tukey’s post-hoc test ***p<1.0×10⁻³ and ^###p<1.0×10⁻³]. Male fracture load analyses for WT vs Daam2^tm1a/tm1a (***) and Daam2^+/tm1a vs Daam2^tm1a/tm1 (^##) had statistically significant differences (one-way ANOVA p=3.0×10⁻⁴, F=15.49, df=16, Tukey’s post-hoc test ***p<1.0×10⁻³ and ^##p<0.01). b) Vertebra biomechanical analyses. Destructive compression testing (Instron 5543 load frame) of caudal vertebrae from WT (n_female=3, n_male=4), Daam2^+/tm1a (n_female=6, n_male=4) and Daam2^tm1a/tm1a (n_female=5, n_male=9) mice. Graphs show yield load, maximum load and stiffness. Data are shown as mean ± SEM. Female yield load analysis for WT versus Daam2^tm1a/tm1a (**) had a statistically significant difference (one-way ANOVA p=6.5×10⁻³, F=8.26, df=13, Tukey’s post-hoc test **p<0.01). Female maximum load analyses for WT versus Daam2^tm1a/tm1a (**) and WT versus Daam2^+/tm1a (*) had statistically significant differences (one-way ANOVA p=2.9×10⁻³, F=10.45, df=13, Tukey’s post-hoc test **p<0.01 and *p<0.05). Male maximum load analysis for WT vs Daam2^tm1a/tm1a (*) had a statistically significant difference (one-way ANOVA p=0.04, F=4.10, df=16, Tukey’s post-hoc test *p<0.05). c) Bone quality analysis from rapid throughput screening mouse knockouts. The graph demonstrates the physiological relationship between bone mineral content and stiffness in caudal vertebrae from P112 female WT mice (n=320). The blue line shows the linear regression (Pearson’s r=0.21, p=1.2×10⁻⁴) and the grey box indicates ± 2 standard deviations (SD). The mean value for female Daam2^tm1a/tm1a [n=2 from initial OBCD screen (Supplementary Note)] mice is shown in orange (−2.14 SD).