Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture

Abstract

The extent to which low-frequency (minor allele frequency (MAF) between 1-5%) and rare (MAF ≤ 1%) variants contribute to complex traits and disease in the general population is mainly unknown. Bone mineral density (BMD) is highly heritable, a major predictor of osteoporotic fractures, and has been previously associated with common genetic variants, as well as rare, population-specific, coding variants. Here we identify novel non-coding genetic variants with large effects on BMD (ntotal = 53,236) and fracture (ntotal = 508,253) in individuals of European ancestry from the general population. Associations for BMD were derived from whole-genome sequencing (n = 2,882 from UK10K (ref. 10); a population-based genome sequencing consortium), whole-exome sequencing (n = 3,549), deep imputation of genotyped samples using a combined UK10K/1000 Genomes reference panel (n = 26,534), and de novo replication genotyping (n = 20,271). We identified a low-frequency non-coding variant near a novel locus, EN1, with an effect size fourfold larger than the mean of previously reported common variants for lumbar spine BMD (rs11692564(T), MAF = 1.6%, replication effect size = +0.20 s.d., Pmeta = 2 × 10(-14)), which was also associated with a decreased risk of fracture (odds ratio = 0.85; P = 2 × 10(-11); ncases = 98,742 and ncontrols = 409,511). Using an En1(cre/flox) mouse model, we observed that conditional loss of En1 results in low bone mass, probably as a consequence of high bone turnover. We also identified a novel low-frequency non-coding variant with large effects on BMD near WNT16 (rs148771817(T), MAF = 1.2%, replication effect size = +0.41 s.d., Pmeta = 1 × 10(-11)). In general, there was an excess of association signals arising from deleterious coding and conserved non-coding variants. These findings provide evidence that low-frequency non-coding variants have large effects on BMD and fracture, thereby providing rationale for whole-genome sequencing and improved imputation reference panels to study the genetic architecture of complex traits and disease in the general population.

Extended Data Figure 1. Discovery single variant meta-analysis

a. Overall study design b. From top to bottom, quantile-quantile plots for the sex-combined single SNV meta-analysis, sex-stratified single SNV meta-analysis (forearm phenotype consists solely of female-only cohorts), and sex-combined single SNV conditional meta-analysis. Plots depicts p-values prior (blue) and after conditional analysis (red). c. From top to bottom, Manhattan plots for sex-combined meta-analysis for lumbar spine BMD, femoral neck BMD, and forearm BMD. Each plot depicts variants from the UK10K/1KG reference panel with MAF > 0.5% across the 22 autosomes (odd=grey, even=black) against the –log10 p-value from the meta-analysis of 7 cohorts (dots). Also depicted is the subset variants from the reference panel that are also present in Estrada et al. (2012) with p value < 5e-6 (diamonds). Variants with MAF < 5% and p < 1.2e-6 are also depicted (red). d. Quantile-Quantile plots for the sex-combined meta-analysis of lumbar spine, femoral neck, and forearm BMD for SNVs present across both exome-sequenced and genome-wide cohorts i.e. SNV absent from all exome-sequenced cohort are not shown. e. Manhattan plot for the Meta-Analysis of Sex-Combined results for Lumbar Spine BMD for SNVs present in exome-sequenced and genome-wide cohorts i.e. SNV absent from all exome-sequenced cohort are not shown (from top to bottom: lumbar spine, forearm and femoral neck BMD).

Extended Data Figure 2. Forest Plots by Cohort for Genome-wide Significant Loci from Discovery Meta-analysis

Forest plots for three BMD phenotypes are shown. Title of each plot includes gene overlapping the SNV and its genomic position on build hg19. P-values are from fixed-effect meta-analysis (see Supplemental Information).

Extended Data Figure 3. Gene Expression in Human and Mouse

a. Quantification of Dock8 expression and its temporal pattern through RNA-seq in cultured calvarial murine osteoblasts across day 2 through to day 18 of osteoblast development. Bglap is shown for comparison, which encodes osteocalcin a critical protein in osteoblasts. b. Quantification of expression of genome-wide significant genes and their temporal pattern through RNA-seq in cultured calvarial murine osteoblasts across day 2 through to day 18 of osteoblast development. c. Expression of EN1 mRNA in human cells presented as percent of GAPDH mRNA. d. Expression of En1 in control and sdEn1 mice in purified osteoblast culture. For osteoblast marker gene expression, total mRNAs were purified from osteoblast cultures at day 10 and measured using quantitative real-time PCR. mRNA levels were normalized relative to GAPDH mRNA. e. Real-time PCR expression of control and sdEn1 as compared to 18S mRNA in whole vertebral bone extract. All data are shown as mean± SEM. Significance computed by student unpaired t-test.

Extended Data Figure 4. Histological Assessment of En1^Cre–expressing cells in in skeletal cells of the vertebra

a. Lineage history of En1^Cre–expressing cells in skeletal cells of the vertebra. The En1^Cre allele was combined with the R26^LSL-YFP reporter allele and examined using frozen fluorescent immunohistochemistry and alkaline phosphatase (AP) staining. Cell nuclei were detected with DAPI. YFP-expressing cells have expressed CRE (En1) at some time in their history. A. Control animals lacking the R26^LSL-YFP reporter show low background YFP signal (green). B. In En1^Cre/+; R26^LSL-YFP/+ mice YFP-expressing cells are detected in the growth plate chondrocytes of the vertebra (*), trabecular bone lining cells (arrow) and osteocytes (arrow head). Note, high fluorescent background staining in the marrow space. C. The same section is shown stained for AP activity using the fast red substrate. Strong activity is present in the hypertrophic chondrocytes of the growth plate and trabecular bone lining cells (arrow). D. Alignment of the AP and YFP images shows that the trabecular lining cells co-express AP and YFP. b. Colocalization of En1 and Alkaline Phosphatase expression. Images of lumbar vertebrae sections (growth plate and trabecular bone regions, 40x) from two-month old En1^lacZ/+ mice. (see Figure 3b), stained for LacZ and Alkaline phosphatase (AP), false-coloured as indicated. Double-positive cells are indicated by arrows, while single-positive cells are indicated by arrowheads (LacZ+) or asterisks (AP+). Except for some chondrocytes, most AP+ cells are also LacZ+, i.e. express En1. The bone marrow was digitally removed, as it contains no AP+ cells.

Extended Data Figure 5. MicroCT Results for control (En1^flox/+) and self-deleting En1 knockout (sdEn1, En1^Cre/flox) animals

a. Trabecular Bone MicroCT images from Lumbar Vertebra 5. b. Morphological characteristics at lumbar vertebra 4,5, and 6 (from bottom to top). c. Morphological characteristics of left femur trabecular bone and d. left femur cortical bone. e. MicroCT parameter results for the comparison of control type and sdEn1 animals at lumbar vertebra 5, femur trabecula, and femur cortical bone. Horizontal lines denote mean of observations. Significance between control and sdEn1 is calculated using an unpaired t-test.

Extended Data Figure 6. Novel association from 7q31.3

a. Chromatin interaction data from Hi-C performed in H1 ES cells of a 2 Mb region encompassing rs148771817 (red and identified by arrow) and WNT16. b. The left axis denotes the association P-value (red and green lines at P = 1.2 × 10⁻⁵ and 1.2 × 10⁻⁸, respectively). The novel genome-wide significant SNV, rs148771817, within an intron of CPED1, and the lead genome wide-significant SNV rs7776725 upstream to WNT16 (within FAM3C) are in low LD with each other. c. Allele frequency versus absolute effect size (in standard deviations) for forearm BMD of all previously identified genome-wide significant variants (blue) and the novel variant within CPED1 (red), rs148771817 from replication meta-analysis. The blue line denotes the mean of effect sizes for previously reported forearm BMD variants. d. Meta-analysis summary statistics of rs148771817 conditioned on rs7776725.

Extended Data Figure 7. Regional Plots of Genome-Wide Significant Loci from Single-SNV Association Tests for forearm and femoral neck BMD

Each regional plot depicts SNVs within 1 Mb of a locus’ lead SNV (x-axis) and their associated meta-analysis p value (-log10). SNVs are color coded according to r² with the lead SNV (labelled, r² calculated from UK10K whole genome sequencing dataset). Recombination rate (blue line), and the position of genes, their exons and the direction of transcription are also displayed (below plot).

Extended Data Figure 8. Regional Plots of Genome-Wide Significant Loci from Single-SNV Association Tests from Lumbar Spine BMD

Each regional plot depicts SNVs within 1 Mb of a locus’ lead SNV (x-axis) and their associated meta-analysis p value (-log10). SNVs are color coded according to r² with the lead SNV (labelled, r² calculated from UK10K whole genome sequencing dataset). Recombination rate (blue line), and the position of genes, their exons and the direction of transcription are also displayed (below plot).

Extended Data Figure 9. Region-based association tests using skatMeta for windows of 30 SNVs and window step of 20 SNVs

a. From top to bottom, quantile-quantile plots for forearm BMD (FA), femoral neck BMD (FN), and lumbar spine (LS) BMD. For each MAF range considered (<5% or <1%), analysis was conducted across all variants, variant overlapping coding exons, and variants with GERP++ score > 1. b. From top to bottom, Manhattan plots forearm BMD, femoral neck BMD, and lumbar spine BMD. For each MAF range considered (<5% or <1%), analysis was conducted across all variants, variant overlapping coding exons, and variants with GERP++ score > 1. Blue and red lines at genome-wide suggestive [P = 1.2 × 10⁻⁶] and genome-wide significant [P = 1.2 × 10⁻⁸] thresholds, respectively.

Extended Data Figure 10. Single Variant Analysis of Signals from Region-based Tests

a. Drop-one SNV and drop-one cohort for genome-wide significant 30 SNV windows for forearm BMD from skatMeta analysis. (A, B) For given 30 SNV window, the –log10(p) of skatMeta test for 29 SNVs, excluding (i.e. dropping) the SNV at position labeled on x-axis. (C, D) For given 30 SNV window, the –log10(p) of skatMeta test for 4 cohorts, excluding (i.e. dropping) cohort labeled on x-axis. b. Drop-one SNV and drop-one cohort for genome-wide significant 30 SNV windows for femoral neck BMD for skatMeta analysis. (A) For given 30 SNV window, the –log10(p) of skatMeta test for 29 SNVs, excluding (i.e. dropping) the SNV at position labeled on x-axis. (B) For given 30 SNV window, the –log10(p) of skatMeta test for 4 cohorts, excluding (i.e. dropping) cohort labeled on x-axis. c. Regional view of CPED1/WNT16 locus for forearm BMD. In top panel, significant SNVs from single variant meta-analysis (rs148771817 and rs79162867, in blue) overlap significant regions found using region-based test (red bars).

Figure 1. Association signals near Engrailed-1 for lumbar spine BMD

a, A topological domain includes associated variants and EN1, and chromatin interaction analysis with paired-end tag sequencing (ChIA-PET for CTCF in MCF-7 cell line) suggests a smaller interacting region containing EN1, and three genome-wide significant variants for lumbar spine BMD (in red). b, Association signals at the EN1 locus (green line at P = 1.2×10⁻⁸) for lumbar spine BMD. Red circles and triangles represent results from discovery and combined discovery and replication using fixed-effects meta-analysis (see Supplementary Information), respectively. c, Allele frequency versus absolute effect size for lumbar spine BMD for previously identified variants (blue) and the three EN1 novel variants (red). The blue line denotes the mean of previously reported effect sizes.

Figure 2. Genome-wide features of association signals

a, Box plots of the effect sizes of genome-wide significant SNVs (P < 1.2×10⁻⁸), pruned for LD (r² < 0.2) by MAF bin for discovery cohorts. Grey bars represent the values of beta not observed and for which we lack statistical power to observe (at α ≤ 1.2×10⁻⁸ and power ≥ 0.8). P-values per phenotype are from the non-parametric trend test across MAF bins (see Supplementary Information). b, Proportion of SNVs passing an FDR q-value 0.05 across different annotation features in discovery cohorts (green) vs. matched control variants (red). The rightmost three panels show enrichment across a range of evolutionary constraint scores, where green denotes SNVs above the threshold and red denotes variants below the threshold. Bars represent standard error (for methods refer to Supplementary Information).

Figure 3. Mouse En1 Functional Experiments

a, Left: Quantitative expression of En1 and its temporal pattern (RNA-seq) in cultured calvarial murine osteoblasts (n=3 per time point). Right: Confirmation of the expression of En1 in a separate RT-PCR experiment of cultured calvarial murine osteoblasts and lack of expression in osteoclasts matured from bone marrow derived precursor cells (Positive controls for osteoblasts (osteocalcin) and osteoclast (RANK) are also shown). b, Representative sections from lumbar vertebra 2 show the growth plate and bone marrow (GP and BM, left), cortical bone (CB, middle), and trabecular bone (TB, right) at 40x magnification from En1^lacZ/+ adult mice (n = 2) stained for β-gal activity (LacZ blue, En1+ cells) and alkaline phosphatase (AP, red late chondrocytes and actively calcifying tissues). In the periosteum (PO), all the LacZ+ cells were AP+; some AP- BM cells expressed LacZ. Some AP- proliferative chondrocytes in the GP expressed lacZ+, whereas most AP+ hypertrophic chondrocytes expressed LacZ. Some AP- osteocytes (Ocy) in CB and TB were LacZ+. c, Left: Histomorphometry images of lumbar vertebrae 5 show decreased trabecular bone volume and increased bone surface area occupied by osteoclast cells when comparing En1^Cre/flox (self-deleted En1, sdEn1) mutants and En1^flox/+ control mice. Right: Reconstructed μCT images show the mineral density in a control and an sdEn1 animal (right panels). d, Micro-CT (μCT) and histomorphometry measures within sdEn1 (n = 5) and controls (En1^lox/+, n = 6). By μCT, sdEn1 mutants exhibit decreased L5 trabecular number (Tb.N) and thickness (Tb.Th), as well as deceased bone volume fraction (BV/TV). Using histomorphometry, sdEn1 mutants exhibit increased osteoclastic area (TRAP/BS). Average for each measure denoted by solid horizontal line. P-value computed using paired t-test.