Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture - PubMed (original) (raw)

. 2015 Oct 1;526(7571):112-7.

doi: 10.1038/nature14878. Epub 2015 Sep 14.

Vincenzo Forgetta 1 2, Yi-Hsiang Hsu 3 4 5, Karol Estrada 4 5 6 7, Alberto Rosello-Diez 8, Paul J Leo 9, Chitra L Dahia 10 11, Kyung Hyun Park-Min 12, Jonathan H Tobias 13 14, Charles Kooperberg 15, Aaron Kleinman 16, Unnur Styrkarsdottir 17, Ching-Ti Liu 18, Charlotta Uggla 19, Daniel S Evans 20, Carrie M Nielson 21 22, Klaudia Walter 23, Ulrika Pettersson-Kymmer 24 25, Shane McCarthy 23, Joel Eriksson 19 26, Tony Kwan 27, Mila Jhamai 6, Katerina Trajanoska 6 28, Yasin Memari 23, Josine Min 14, Jie Huang 23, Petr Danecek 23, Beth Wilmot 29 30, Rui Li 1 2, Wen-Chi Chou 3 4, Lauren E Mokry 2, Alireza Moayyeri 31 32, Melina Claussnitzer 3 4 5 33, Chia-Ho Cheng 3, Warren Cheung 27 34, Carolina Medina-Gómez 6 28 35, Bing Ge 27, Shu-Huang Chen 27, Kwangbom Choi 36, Ling Oei 6 28 35, James Fraser 37, Robert Kraaij 6 28 35, Matthew A Hibbs 36 38, Celia L Gregson 39, Denis Paquette 37, Albert Hofman 28 35, Carl Wibom 40, Gregory J Tranah 21 22, Mhairi Marshall 9, Brooke B Gardiner 9, Katie Cremin 9, Paul Auer 41, Li Hsu 15, Sue Ring 42, Joyce Y Tung 16, Gudmar Thorleifsson 43, Anke W Enneman 6, Natasja M van Schoor 44, Lisette C P G M de Groot 45, Nathalie van der Velde 6 46, Beatrice Melin 40, John P Kemp 9 14, Claus Christiansen 47, Adrian Sayers 39, Yanhua Zhou 18, Sophie Calderari 48 49, Jeroen van Rooij 6 35, Chris Carlson 15, Ulrike Peters 15, Soizik Berlivet 37, Josée Dostie 37, Andre G Uitterlinden 6 28 35, Stephen R Williams 50, Charles Farber 50, Daniel Grinberg 51 52 53, Andrea Z LaCroix 54, Jeff Haessler 15, Daniel I Chasman 4 55, Franco Giulianini 55, Lynda M Rose 55, Paul M Ridker 4 55, John A Eisman 56 57 58, Tuan V Nguyen 56 58, Jacqueline R Center 56 58, Xavier Nogues 59 60 61, Natalia Garcia-Giralt 59 60, Lenore L Launer 62, Vilmunder Gudnason 63 64, Dan Mellström 19, Liesbeth Vandenput 19, Najaf Amin 65, Cornelia M van Duijn 65, Magnus K Karlsson 66, Östen Ljunggren 67, Olle Svensson 68, Göran Hallmans 25, François Rousseau 69 70, Sylvie Giroux 70, Johanne Bussière 70, Pascal P Arp 6, Fjorda Koromani 6 28, Richard L Prince 71 72, Joshua R Lewis 71 72, Bente L Langdahl 73, A Pernille Hermann 74, Jens-Erik B Jensen 75, Stephen Kaptoge 31, Kay-Tee Khaw 76, Jonathan Reeve 77 78, Melissa M Formosa 79, Angela Xuereb-Anastasi 79, Kristina Åkesson 66 80, Fiona E McGuigan 80, Gaurav Garg 80, Jose M Olmos 81 82, Maria T Zarrabeitia 83, Jose A Riancho 81 82, Stuart H Ralston 84, Nerea Alonso 84, Xi Jiang 85, David Goltzman 86, Tomi Pastinen 27 34, Elin Grundberg 27 34, Dominique Gauguier 48 49, Eric S Orwoll 22 87, David Karasik 3 88, George Davey-Smith 14; AOGC Consortium; Albert V Smith 63 64, Kristin Siggeirsdottir 63, Tamara B Harris 89, M Carola Zillikens 6, Joyce B J van Meurs 6 28, Unnur Thorsteinsdottir 17 64, Matthew T Maurano 90, Nicholas J Timpson 14, Nicole Soranzo 23, Richard Durbin 23, Scott G Wilson 32 71 91, Evangelia E Ntzani 92 93, Matthew A Brown 9, Kari Stefansson 64 94, David A Hinds 16, Tim Spector 32, L Adrienne Cupples 18 95, Claes Ohlsson 19, Celia M T Greenwood 2 34 96 97; UK10K Consortium; Rebecca D Jackson 98, David W Rowe 85, Cynthia A Loomis 99, David M Evans 9 14, Cheryl L Ackert-Bicknell 36, Alexandra L Joyner 8, Emma L Duncan 9 100, Douglas P Kiel 3 4 5 33, Fernando Rivadeneira 6 28 35, J Brent Richards 1 2 32

Affiliations

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

Hou-Feng Zheng et al. Nature. 2015.

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.

PubMed Disclaimer

Figures

Extended Data Figure 1

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

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

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

Extended Data Figure 4. Histological Assessment of _En1Cre_–expressing cells in in skeletal cells of the vertebra

a. Lineage history of _En1Cre_–expressing cells in skeletal cells of the vertebra. The En1Cre allele was combined with the R26LSL-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 R26LSL-YFP reporter show low background YFP signal (green). B. In En1Cre/+; R26LSL-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 En1lacZ/+ 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

Extended Data Figure 5. MicroCT Results for control (En1flox/+) and self-deleting En1 knockout (sdEn1, En1Cre/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

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−5 and 1.2 × 10−8, 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

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 r2 with the lead SNV (labelled, r2 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

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 r2 with the lead SNV (labelled, r2 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

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−6] and genome-wide significant [P = 1.2 × 10−8] thresholds, respectively.

Extended Data Figure 10

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

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−8) 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

Figure 2. Genome-wide features of association signals

a, Box plots of the effect sizes of genome-wide significant SNVs (P < 1.2×10−8), pruned for LD (r2 < 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−8 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

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 En1lacZ/+ 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 En1Cre/flox (self-deleted En1, sdEn1) mutants and En1flox/+ 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 (En1lox/+, 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.

References

    1. Richards JB, Rivadeneira F, Inouye M, et al. Bone mineral density, osteoporosis, and osteoporotic fractures: a genome-wide association study. Lancet. 2008;371(9623):1505–1512. doi: 10.1016/S0140-6736(08)60599-1. - DOI - PMC - PubMed
    1. Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, et al. Multiple genetic loci for bone mineral density and fractures. N Engl J Med. 2008;358(22):2355–2365. doi: 10.1056/NEJMoa0801197. NEJMoa0801197 [pii]\r. - DOI - PubMed
    1. Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, et al. New sequence variants associated with bone mineral density. Nat Genet. 2008;41(1):15–17. doi: 10.1038/ng.284. - DOI - PubMed
    1. Rivadeneira F, Styrkársdottir U, Estrada K, et al. Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat Genet. 2009;41(11):1199–1206. doi: 10.1038/ng.446. - DOI - PMC - PubMed
    1. Duncan EL, Danoy P, Kemp JP, et al. Genome-wide association study using extreme truncate selection identifies novel genes affecting bone mineral density and fracture risk. PLoS Genet. 2011;7(4):e1001372. doi: 10.1371/journal.pgen.1001372. - DOI - PMC - PubMed

MeSH terms

Substances

Grants and funding

LinkOut - more resources