1. Genetics and Genomics
Download icon

Epigenetic profiling of growth plate chondrocytes sheds insight into regulatory genetic variation influencing height

  1. Michael Guo
  2. Zun Liu
  3. Jessie Willen
  4. Cameron P Shaw
  5. Daniel Richard
  6. Evelyn Jagoda
  7. Andrew C Doxey
  8. Joel Hirschhorn
  9. Terence D Capellini  Is a corresponding author
  1. Broad Institute of MIT and Harvard, United States
  2. Harvard Medical School, United States
  3. Harvard University, United States
  4. University of Waterloo, Canada
Research Article
  • Cited 16
  • Views 3,512
  • Annotations
Cite this article as: eLife 2017;6:e29329 doi: 10.7554/eLife.29329

Abstract

GWAS have identified hundreds of height-associated loci. However, determining causal mechanisms is challenging, especially since height-relevant tissues (e.g. growth plates) are difficult to study. To uncover mechanisms by which height GWAS variants function, we performed epigenetic profiling of murine femoral growth plates. The profiled open chromatin regions recapitulate known chondrocyte and skeletal biology, are enriched at height GWAS loci, particularly near differentially expressed growth plate genes, and enriched for binding motifs of transcription factors with roles in chondrocyte biology. At specific loci, our analyses identified compelling mechanisms for GWAS variants. For example, at CHSY1, we identified a candidate causal variant (rs9920291) overlapping an open chromatin region. Reporter assays demonstrated that rs9920291 shows allelic regulatory activity, and CRISPR/Cas9 targeting of human chondrocytes demonstrates that the region regulates CHSY1 expression. Thus, integrating biologically relevant epigenetic information (here, from growth plates) with genetic association results can identify biological mechanisms important for human growth.

https://doi.org/10.7554/eLife.29329.001

eLife digest

Humans vary considerably in height, a trait that is partly inherited from each individual's parents. Studies have identified hundreds of small changes in DNA that contribute to differences in human height. These small changes often swap out just one of the four letters that make up the DNA code. Some changes occur within genes, yet most occur in stretches of DNA that do not contain genes. These stretches of DNA likely control whether genes important for the growth of the skeleton are switched on or off.

Cartilage cells, also known as chondrocytes, are important for height. These cells are found near the end of bones in the growth plates, where the bones grow during childhood and adolescence. The on and off switches for growth genes in chondrocytes are unknown. Identifying these genetic switches could help scientists understand how hundreds of small changes in DNA help determine how tall a person will be.

Now, Guo, Liu, Willen et al. identify thousands of switches that turn on and off genes in chondrocytes. In their experiments, chondrocytes were removed from mouse bones and a type of genetic sequencing called ATAC-seq was used to identify the stretches of DNA that act as on/off switches for genes in these cells. Further analysis revealed that many of these same on/off switches occur in human chondrocytes too. Importantly, the experiments showed that many of the small changes in DNA that contribute to differences in human height are also found within these DNA switches.

Guo, Liu, Willen et al. next tested how these switches and height-linked genes interact, and found, for example, that one switch acts to shut off the gene for a protein called chondroitin sulfate synthase 1. People with mutations in this gene can have unusually short stature. However, people with DNA changes in its nearby switch can be taller or shorter than average based on how the DNA influences the switch. More studies might help scientists understand the evolutionary significance of these DNA changes. They will also help determine if the genetic switches in chondrocytes contribute to diseases that affect the bones and joints, like bone cancer or arthritis.

https://doi.org/10.7554/eLife.29329.002

Introduction

Human height, a main outcome of skeletal growth, is the product of many biological and environmental interactions spanning numerous cell and tissue types acting during pre- and post-natal development. The growth plate, located at the distal ends of long bones such as the femur, is among the most important tissues influencing height. In the growth plate, condensations of mesenchymal cells differentiate into chondrocytes in the process of endochondral ossification. Chondrocytes in the growth plate reside in longitudinally oriented columns divided in zones (resting, proliferative, pre-hypertrophic, and hypertrophic) each with important functions in allowing the bone to elongate and mature (see Liu et al., 2017 and Samsa et al., 2017 for reviews) (Liu et al., 2017; Samsa et al., 2017). This process of growth plate chondrocyte differentiation and maturation promotes elongation of long bones, which ultimately determines much of human height (Baron et al., 2015). These processes are under the influence of many extrinsic and intrinsic signals (Samsa et al., 2017; Lui et al., 2012; Kronenberg, 2003).

Height has been studied for centuries as a model genetic trait, as it is easily measured and highly heritable (typically, 70–90% of variation in height within a population is attributable to genetic variation) (Silventoinen et al., 2003; Perola et al., 2007; Visscher et al., 2006; Visscher et al., 2007; Galton, 1886). Height and body proportions are likely selected and thus are interesting from an evolutionary standpoint (Turchin et al., 2012). Studies of the genetics of height can shed insight into the mechanisms of childhood growth and into developmental biology of the skeleton (Baron et al., 2015; Dauber et al., 2014). Additionally, abnormalities of growth—short stature or overgrowth—are among the most common childhood disorders and their pathophysiology is poorly understood. Thus, genetic studies of height could further our understanding of the genetic basis of skeletal disease, evolution, biology, and normal human growth and development.

Recently, genome-wide association studies (GWAS) have been used to study the genetics of height, as well as many other traits and disorders (Wood et al., 2014). This approach systematically tests each common variant in the genome for association with a phenotype of interest, such as adult height, with the goal of identifying relevant biology, either from the identities of the genes near the associated variants, or by deducing mechanism from the associated variants themselves. GWAS for height are among the largest conducted thus far (sample sizes > 250,000) and have identified more associations than any other phenotype, with nearly 700 regions of the genome (or loci) robustly associated with height (Wood et al., 2014). The finding of so many loci influencing height, with the number expected to increase with larger sample sizes, supports the assertion that height is highly polygenic. Although each of the variants at these loci have only a small effect on height (typically much less than 0.5 cm per allele) (Wood et al., 2014), the genes within the associated loci as a group highlight important mechanisms influencing growth. For example, height variants are enriched for genes implicated in numerous biological processes relevant to growth, including embryonic stem cell function, long bone development, and spine length, among a number of other processes (Wood et al., 2014). Height variants are also preferentially located near genes that are differentially expressed in the growth plate (Lui et al., 2012).

However, as with all GWAS, moving from genetic associations to biological mechanisms at any individual locus can be challenging. There are two primary reasons for this difficulty. First, there is extensive linkage disequilibrium (LD) in the human genome, meaning that the genotypes at nearby genetic variants are often tightly correlated (Gabriel et al., 2002). The consequence of LD is that at any particular locus, there will be many variants that have indistinguishable associations, and determining which variant(s) are in fact causal (rather than associated due to LD with a causal variant) remains difficult. Second, most GWAS variants (>80%) reside in non-coding regions (Gusev et al., 2014), with no nearby coding variants that can account for the associations. Non-coding variants are challenging to interpret because there is no universal regulatory code to infer function/mechanism, and these variants can act to influence gene expression at long distances and in highly tissue-specific manners (Bernstein and ENCODE Project Consortium, 2012). Various approaches have been pursued for identifying the causal variants and/or mechanisms, including overlapping associated variants with epigenetic datasets in appropriate tissues/cell-types, statistical genetic methods leveraging subtle variations in LD patterns (called ‘fine mapping’), and in vitro and in vivo functional assays (Tak and Farnham, 2015). It is also becoming increasingly clear that GWAS variants tend to fall in regions that have particular epigenetically defined states (such as enhancers or repressors) (Finucane et al., 2015; Farh et al., 2015). Thus, the overlap between GWAS and the appropriate epigenetic annotations can help refine the identities of causal variants, helping to decipher mechanism. However, epigenetic data may not exist for the biologically relevant cell type, and there is no universally successful approach to translating non-coding associations from GWAS into biological mechanisms.

Although the large number of associated variants offers advantages in making biological inferences for the set of associated variants as a whole (Lango Allen et al., 2010), defining functional mechanisms of individual height GWAS variants faces specific challenges. Height variants likely act in many tissues and stages of development, including not only the growth plate but also embryonic cells, endocrine tissues such as the pituitary, and bone, among others (Wood et al., 2014). Human growth is also a whole-organism phenotype and is less amenable to study in cell-intrinsic models, making functional approaches particularly challenging. Finally, skeletal tissues can be highly heterogeneous and challenging to isolate because in their mature form they are surrounded by a hard extracellular matrix and strongly adhering soft tissues such as tendons, ligaments, and muscles. These factors make it difficult to isolate the large numbers of biologically relevant cells that are needed for functional experiments or for traditional epigenomic assays such as ChIP-seq (Furey, 2012).

In this study, we profile the epigenetic landscape of mouse growth plate chondrocytes using a new approach called ATAC-seq (Buenrostro et al., 2015, 2013), which allows for accurate profiling of open chromatin regions using relatively few cells (~50,000). We then use the epigenetic landscape in growth plate chondrocytes to gain new insights into the mechanisms of non-coding elements in the genome influencing expression of growth plate genes. Finally, we use these epigenetic marks, in the context of statistical fine-mapping results and gene expression profiles, to better understand height GWAS variants and the mechanisms by which they may function to influence height in humans. These studies illustrate a path that can help decipher mechanism from GWAS, even when the relevant cell type is relatively inaccessible in humans.

Results

Generation of mouse femoral growth plate epigenetic profiles

In order to identify open chromatin regions involved in long bone development, we performed ATAC-seq (Buenrostro et al., 2015, 2013) on E15.5 (Theiler Stage 23) mouse proximal and distal femoral growth zones, regions comprised predominately of the growth plate proper, the overlying perichondrium, and the cartilaginous super-structures that dictate adult bony morphology (Materials and methods). This stage of mouse embryogenesis corresponds to approximately E48-51 (Carnegie Stage 19) of human development and reflects a developmental window when the growth plate proper possesses all chondrocyte zones (resting, proliferative, pre-hypertrophic, hypertrophic zone) and is undergoing rapid cell proliferation and differentiation during endochondral ossification. Given that it remains unclear which sub-region(s) (i.e. chondrocyte zones) of the growth plate are under genetic influence to drive height variation, we performed our experiments on tissues that include all chondrocyte zones (including perichondrium) and cartilaginous tissues, excluding regions indicative of the osteogenic invasion in the developing bone collar (Figure 1a). We used ATAC-seq to isolate open chromatin regions from the DNA, which were subsequently subjected to next-generation sequencing (Materials and methods; Supplementary file 1). This approach yielded 28,257 proximal femur peaks and 27,958 distal femur peaks, with 23,764 shared peaks, 4493 peaks specific to the proximal femur, and 4194 specific to the distal femur (Figure 1b). The proximal and distal peaks span approximately 0.79% and 0.76% of the genome, respectively. We also assessed the genomic distribution of proximal and distal femur ATAC-seq peaks and found a characteristic enrichment near gene transcriptional start sites and more distal intronic and intergenic non-coding sequences (Figure 1—figure supplement 1).

Figure 1 with 2 supplements see all
ATAC-seq on in vivo collected mouse femoral growth plate zones.

(A) Pipeline of ATAC-seq on the embryonic mouse femur. Proximal (red) and distal (blue) samples were extracted from E15.5 femoral growth plates, as labeled by Col2a1 expression (insets). Tissues were harvested just proximal (red dotted lines in inset) or distal (blue dotted lines in inset) to the forming bone collar (outlined by dashed black lines). Cells were then treated to generate a single-cell suspension. Approximately 50,000 cells per sample were utilized for subsequent ATAC-seq. Open chromatin regions were then identified from the ATAC-seq data. ATAC-seq protocol adapted from Buenrostro et al. (2013). (B) Venn diagram shows the number of ATAC-seq peaks in the proximal and distal samples. (C) ATAC-seq peaks (after subtracting peaks shared in brain ATAC-seq) were analyzed for enrichment in gene sets using GREAT. Horizontal bars show enrichment p values (shown as –log10(p value)) for Human Phenotypes. The union of proximal and distal femur peaks was used.

https://doi.org/10.7554/eLife.29329.003

Since a main goal was to understand the gene regulatory functions of variants influencing human height, we next mapped the locations of each mouse femur open chromatin region to the human genome (Materials and methods). Our liftover resulted in 24,805 proximal femur peaks and 24,788 distal femur peaks, representing approximately 87.7% and 88.6% of sites identified in mouse, respectively. These sites that were capable of liftover to human likely reflect relatively conserved genes involved in vertebrate musculoskeletal development (Shubin et al., 1997; Petit et al., 2017). The union of the two sets contained 25,829 peaks (Figure 1—figure supplement 1).

To determine whether ATAC-seq peaks are enriched near gene sets/pathways with specific functions, we performed an enrichment analysis using Genomic Regions Enrichment of Annotations Tool (GREAT) (McLean et al., 2010), which identifies curated gene sets enriched near epigenetic peaks. Analyses of each proximal, distal, and union femur datasets yielded general GO Biological Process and GO Molecular Processes terms related to basic cell biological processes, as well as phenotypes related to a number of biological tissues including skeletal tissues (Supplementary file 24). They also yielded phenotype terms related to general skeletal development, chondrogenesis, and long bone skeletal phenotypes, such as ‘Metaphyseal widening’, ‘Abnormality of the femoral epiphyses’, and ‘Abnormality of the femoral head’ (Figure 1c). Thus, ATAC-seq datasets generated on each specific femoral growth plate in mice reflect in part highly relevant anatomy, as well as general aspects of skeletal development in humans (see Discussion).

The ATAC-seq peaks we identified also highlight biological processes that are important in the growth plate. For example, GREAT analyses uncovered enrichment for several relevant GO Biological Processes such as ‘Cartilage condensation’, ‘Positive regulation of chondrocyte differentiation’, and ‘Chondrocyte proliferation’ (Supplementary file 24). Additionally, we illustrate ATAC-seq peaks at several well-known growth plate genes including Col2a1, Fgfr1, and Pth1r (Figure 1—figure supplement 2) (Long and Ornitz, 2013).

ATAC-seq peaks overlap other skeletal and chondrocyte epigenetic data sets

We next wanted to validate that our ATAC-seq datasets reflect open chromatin regions important to mouse chondrocyte biology. To do so, we analyzed a previous dataset that used ChIP-seq to profile Sox9 binding in mouse rib chondrocytes (Ohba et al., 2015). We found extensive overlap between Sox9-occupied sites and our open chromatin ATAC-seq peaks. For example, 86.4% of Class I Sox9 elements (i.e. elements clustering near transcriptional start sites of highly expressed non-chondrocyte-specific genes), 44.6% of Class II Sox9 elements (i.e. elements directing chondrocyte-related gene activity through direct Sox9 dimer binding), and 93.8% of Sox9 super-enhancer regions (i.e. enhancer complexes involved in chondrocyte cell identity) overlap with our femur ATAC-seq peaks (p<0.001 for all three sets, Materials and methods, Figure 2a).

Mouse proximal and distal ATAC-seq datasets capture chondrocyte biology in mouse and humans.

(A) Proportion of mouse Sox9 Class I, Class II, and super enhancer sites (Ohba et al., 2015) overlapping with ATAC-seq peaks. (B) Proportion of mouse rib chondrocyte H3K27ac, H3K36me3, H3K4me2, H3K4me3, p300, and RNA Pol II sites overlapping with ATAC-seq peaks. (C) Proportion of human embryonic limb E47 H3K27ac (Cotney et al., 2013) and human BMDC H3K27ac sites (Kundaje et al., 2015) overlapping with ATAC-seq peaks. (D) ATAC-seq peaks at the Matn3 locus. Published rib chondrocyte peaks for H3K4me2, H3K4me3, and Sox9 are shown in shades of orange. Proximal ATAC-seq peaks are shown in red, with called peaks shown as bars below. Distal ATAC-seq peaks are shown in blue, with called peaks shown as bars below. The locations of previously published chondrocyte enhancers are shown as black bars.

https://doi.org/10.7554/eLife.29329.006

We also examined overlap between our dataset and various histone modifications in the mouse rib chondrocytes (Materials and methods) (Figure 2b) (Ohba et al., 2015). We found that 85.1% of RNA Pol II transcriptionally active promoters, 56.9% of H3K4me2 active promoter sites, 86.7% of H3K4me3 active promoter sites, 87.4% of H3K27ac active enhancer sites, and 59.4% of p300 histone acetyltransferase active enhancer sites overlapped with femur ATAC-seq peaks. Conversely, our ATAC-seq datasets displayed relatively lower overlap with markers of gene repression; that is, only 33.7% of H3K27me3-repressed enhancer sites overlap with femur ATAC-seq signals (Figure 2b). Based on permutation analyses, the enrichments of all these histone modification sets were significant at p<0.001. Thus, we found that femoral ATAC-seq peaks show marked evidence of active chromatin state and/or gene transcription.

We also wanted to determine if our ATAC-seq peaks reflect regulatory elements involved in human limb development and chondrocyte biology. First, we used a ChIP-seq dataset comprised of H3K27ac peaks found in E47 human limb buds (equivalent to E13.5 mouse limb buds) that are shared with mouse limb buds (Cotney et al., 2013). We found that 45.2% of these E47 human limb bud H3K27ac peaks overlapped with ATAC-seq peaks, indicating that many of our mouse femur ATAC-seq peaks reflect active human limb enhancers involved in chondrogenesis (p<0.001, Figure 2c). Second, we used an H3K27ac ChIP-seq dataset generated on human adult cultured bone marrow mesenchymal stem-cell-derived chondrocytes (BMDCs) (Herlofsen et al., 2013; Kundaje et al., 2015) and found that 13.3% of human adult BMDC H3K27ac peaks overlapped our mouse distal femur dataset (p<0.001) (Figure 2c), which indicates that our femur ATAC-seq dataset also captures cell-intrinsic properties of chondrogenesis. However, the greater overlap seen in the limb bud dataset as compared to the BMDCs suggests that our dataset may be capturing epigenetic properties present in developing in vivo tissues that are poorly modeled by cell culture systems.

Our ATAC-seq data were also able to identify 14 of 26 previously-discovered active regulatory elements residing near key cartilage genes detected at different time-points and using different assays (e.g. see Figure 2dSupplementary file 5). This suggests that mouse femoral ATAC-seq regions have in vivo enhancer activity as determined in independent assays. Overall, these data indicate that ATAC-seq datasets generated on the mouse femoral growth plates match previously expected patterns of regulatory usage important to chondrocyte biology.

Human height variants are enriched in femoral open chromatin regions

We next sought to use our ATAC-seq dataset to gain insight into height-associated variants identified by GWAS. To do so, we intersected our orthologous human ATAC-seq peaks with 688 height GWAS loci from Wood et al. (2014) (nine loci were excluded; see Materials and methods) (Wood et al., 2014). Since we do not know which variants are causal and which are associated with height due to LD, for each locus, we identified potentially causal ‘proxy’ SNPs: those with a correlation of r2 >0.5 to the lead SNP based on the European subset of 1000 Genomes Phase 3 (Auton et al., 2015). The number of proxy SNPs ranged from 1 to 3548 for each of these 688 loci (Figure 3—figure supplement 1).

Across the 688 height loci, we found that 317 loci (or ~46%) contained at least one variant that overlapped with an ATAC-seq peak. To test the significance of the observed overlap, we applied GoShifter (see Materials and methods), which performs a shuffling-based approach to evaluate the overlap between an epigenetic data set and GWAS data. We found strong enrichments for the ATAC-seq peak set and height GWAS data (p=0.0059) (Figure 3—figure supplement 1a). We repeated this analysis using the Sox9 binding and histone modifications in mouse rib chondrocytes, human BMDCs, and mouse embryonic limb datasets from above, as well as DNaseI hypersensitivity sites from mouse brain, liver, and lung (Bernstein and ENCODE Project Consortium, 2012; Ohba et al., 2015; Cotney et al., 2013; Yue et al., 2014). Only the mouse rib chondrocyte H3K4me2 peaks reached statistical significance (p=0.0237), and was far less significant than our ATAC-seq data. This suggests that our ATAC-seq data captures functional genetic variation influencing height that is not captured by the other datasets. Also, as no enrichment was observed for other mouse tissues from ENCODE, this suggests that the observed enrichment is unlikely due to artifacts stemming from our use of mouse tissues (e.g. from conserved regions being lifted over).

Figure 3 with 1 supplement see all
Height GWAS loci are enriched in ATAC-seq peaks.

(A) Enrichment of height GWAS loci in ATAC-seq peaks using GoShifter. Histogram (light blue bars) shows the number of overlaps observed in the 100,000 GoShifter permutations. Dotted line shows the number of overlaps observed in the ATAC-seq data. (B) GoShifter analysis p values for the ATAC-seq data, mouse rib chondrocyte Sox9 binding and histone modifications, human embryonic limb E47 H3K27ac, H3K27Ac of BMDCs, and brain, liver, and lung DNaseI hypersensitivity site samples from mouse ENCODE (C) Enrichment of height GWAS loci in ATAC-seq peaks using random matched GWAS loci. Histogram (light blue bars) shows the number of overlaps observed in 1000 matched random loci. Dotted line shows the number of overlaps observed in the ATAC-seq data.

https://doi.org/10.7554/eLife.29329.007

To further demonstrate the specificity of our ATAC-seq peaks for capturing height GWAS signals, we repeated the GoShifter analysis by examining the enrichment of GWAS results for other complex traits/diseases at our ATAC-seq peaks. These other traits/diseases were chosen as they have no obvious connection to height and included age of menarche (Day et al., 2017), body mass index (Locke et al., 2015), coronary artery disease risk (Schunkert et al., 2011), LDL cholesterol (Teslovich et al., 2010), mean corpuscular volume (van der Harst et al., 2012), schizophrenia risk (Ripke and Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014), and type 2 diabetes risk (Morris et al., 2012). These additional traits/diseases all demonstrated little to no enrichment at ATAC-seq peaks (Figure 3—figure supplement 1), suggesting that the growth plate ATAC-seq peaks are specifically capturing height/growth biology.

As a second approach to evaluate overlap between our ATAC-seq peak set and height GWAS data, we generated 1000 sets of random loci which were matched to the actual height GWAS loci based on factors including SNP density, minor allele frequency, number of proxy SNPs in the locus, and gene density (Materials and methods) (Pers et al., 2015). For each of the 1000 matched sets of random loci, we repeated the overlap with the ATAC-seq peaks. As none of the matched sets had as much overlap as our height GWAS loci, the enrichment we observed for the height GWAS data is highly significant (p<0.001; z-score: 11.88) (Figure 3c).

Notably, we observed that the actual number of SNPs that overlap an ATAC-seq peak for each set was higher than the number of represented loci. In total, there were 928 variants that overlapped with an ATAC-seq peak; these 928 variants were distributed across 317 loci (Supplementary file 6; Figure 3—figure supplement 1). This finding of more than one overlapping variant per locus may indicate multiple causal variants at each locus, as has been shown for other traits (Roman et al., 2015).

Fine-mapped height variants overlap with femoral growth plate open chromatin regions

An alternative method to narrow down potentially causal signals within GWAS loci is to use statistical fine-mapping, which leverages the association statistics and patterns of LD at each locus to identify which variants at a given locus are most likely to be causal. We applied PICS (Farh et al., 2015) to perform statistical fine-mapping of the (Wood et al., 2014) height GWAS results (Materials and methods). For each of 688 height-associated loci, we calculated credible set probabilities for all SNPs with an r2 >0.5 to the lead variant based on the European subset of 1000 Genomes Phase 3 (Auton et al., 2015). For each locus, we then determined the 95% credible set, which represents the minimum set of SNPs at that locus such that the truly causal SNP has at least a 95% probability of being in the set. Across 688 loci, this yielded a 95% credible set of 32,846 SNPs, which represents approximately a 48.2% reduction from the total number of variants in the GWAS loci. Intersecting these credible set SNPs with femur open chromatin regions revealed 468 overlapping variants spanning 192 of the 688 loci (Figure 4a).

Fine-mapping of GWAS height variants.

(left): Larger light green circle represents all height GWAS SNPs and their proxies (r2 >0.5). Smaller green circle represents the proportion of these SNPs in the PICS 95% credible set. (right) Larger light blue circle represents the 688 height GWAS loci. Medium blue circle represents 317 height GWAS loci overlapping an ATAC-seq peak. Smaller dark blue circle represents 192 loci where a PICS 95% credible set variant overlaps an ATAC-seq peak.

https://doi.org/10.7554/eLife.29329.009

A subset of human height variants in femoral open chromatin regions reside near genes differentially expressed in the growth plate

A previous study profiled expression of mouse growth plate genes and identified a set of 427 genes that were differentially expressed in the growth plate (Lui et al., 2012). That study also revealed that height GWAS loci (using a previous smaller height GWAS) are enriched near differentially expressed growth plate genes (Lui et al., 2012; Lango Allen et al., 2010). We expanded on those results using the most recent height GWAS from Wood et al. (2014) and indeed found that height GWAS loci are enriched near 427 differentially expressed growth plate genes (within 100 kb) when compared to the 1000 matched random GWAS loci (p<0.001, z-score: 8.54, Figure 5a).

Height GWAS loci are enriched in differentially expressed growth plate genes.

(A) Enrichment of height GWAS loci in differentially expressed growth plate genes. Dotted line shows number of GWAS loci with a nearby differentially expressed growth plate gene (±100 kb). Histogram (light blue bars) shows the number of overlaps observed in 1000 matched random loci. (B) Enrichment of height GWAS loci in differentially expressed growth plate genes and ATAC-seq peaks. Same as (A), except showing overlap for GWAS loci overlapping differentially expressed growth plate gene (±100 kb) and ATAC-seq peak.

https://doi.org/10.7554/eLife.29329.010

Since we detected enrichment of height GWAS variants near femoral open chromatin regions (see Figure 3) and reasoned that many of the variants likely function by influencing expression of growth plate genes, we next identified the number of GWAS variants that fit the criteria of residing in an open chromatin region near a differentially expressed gene (within 100 kb). We again used the 688 GWAS loci from Wood et al. (2014) and their proxy SNPs (r2 >0.5) and found that 46/688 loci (or ~6.7%) contained at least one GWAS variant that overlapped with an ATAC-seq peak and resided near a differentially expressed gene (Supplementary file 6). Using the same 1000 sets of matched random loci as described above, we found that these overlaps were highly significant (p<0.001, z-score: 6.36, Figure 5b). These loci likely represent combinations of functional variants that overlap an epigenomic element in growth plates that then modulates expression of nearby genes influencing height.

Wood et al., had previously observed that height GWAS SNPs that represent cis eQTLs in whole blood are enriched for expression in cartilage (Wood et al., 2014). Here, we sought to evaluate whether differentially expressed growth plate genes that are also eQTLs genes in whole blood are enriched at ATAC-seq peaks. Among the 427 differentially expressed growth plate genes, 157 had previously been identified as eQTL genes in whole blood (Westra et al., 2013). Using the GoShifter permutation analysis as described above, we demonstrated that the corresponding eQTL variants for these 157 genes are significantly enriched at ATAC-seq peaks (p=0.0004)

Of the 46 loci in the ATAC-seq femur comparison that fall near genes differentially expressed in growth plates, 44 different genes have been identified as possibly being modulated by these putative regulatory variants (note, some genomic loci have multiple genes and/or GWAS loci) (Supplementary file 6). These genes represent major signaling molecules (e.g. IHH, HHIP), transcription factors (e.g. RUNX2), extracellular matrix factors (e.g. CHSY1, EXTL1), among other factors involved in the growth plate in both humans and mice.

Known and novel targets of human height variants in putative chondrocyte regulatory regions

We next examined some of the most compelling loci from a genetic and mechanistic perspective. We filtered for loci where a height GWAS variant in the PICS 95% credible set overlapped an ATAC-seq peak, and was near a differentially expressed growth plate gene (±100 kb). In total, 59 variants distributed across 26 loci met these criteria. These 59 variants thus have strong evidence for being potentially causal variants and might be prioritized for downstream study (Supplementary file 6).

For example, at the CHSY1 locus, four variants (rs8042551, rs11639408, rs9920291, and rs3911964) met these criteria. Disruption of Chsy1 in mice and humans leads to severe skeletal phenotypes including long bone length reductions (Wilson et al., 2012; Temtamy et al., 1998; Tian et al., 2010; Li et al., 2010). Of the four variants, rs9920291 is predicted to alter HOXD13 binding, as suggested by data generated using universal protein binding microarray technology that demonstrates that the eight base-pair sequence matching the immediate rs9920291 locus disrupts in vitro binding by HOXD13 (Figure 6a, Supplementary file 7), a homeodomain transcription factor that is expressed in the growth plate (Kuss et al., 2014; Reno et al., 2016) and has known activating and repressive roles in skeletal development (Johnson and Tabin, 1997). We performed a series of functional assays to test the effects of allelic variation at rs9920921. First, using a luciferase reporter assay in T/C-28a2 human chondrocytes, we found that the open chromatin region surrounding rs9920291 acts as a repressor of expression (p=7.83×10−6 and 3.6 × 10−4 for variants carrying the reference and alternate alleles respectively) (Figure 6b). In addition, the alternate T allele, (i.e. the height-increasing allele of rs9920291) makes the regulatory element act as a weaker repressor (p=0.0014). This is consistent with data from GTEx (Lonsdale J, GTEx Consortium, 2013), where rs9920291 acts as an eQTL for CHSY1 in transformed fibroblasts (p=1.0×10-7; Figure 6—figure supplement 1), with the T allele correlating with increased CHSY1 expression levels. Additionally, we leveraged the fact that in HEK-293FT cells, rs9920291 is heterozygous, along with several additional SNPs in the coding region of CHSY1. Assaying a heterozygous exonic coding SNP (rs28364839) and a 3’UTR SNP (rs11433), we demonstrate that rs9920291 displays allelic skew in expression (p<0.001 for both assayed variants; Figure 6d).

Figure 6 with 2 supplements see all
Intersections of GWAS height variants with ATAC-seq data identify novel putatively causal variants.

(A) ATAC-seq height GWAS variant intersections at the CHSY1 locus refine height GWAS loci to fewer putatively causal variants. Height GWAS SNPs (r2 >0.5 to lead SNP) are shown as black bars, and SNPs in the PICS 95% credible set are shown as red bars. Published H3K27ac ChIP-seq in E47 human limb buds and BMDCs are shown in orange and purple respectively. ATAC-seq peaks lifted over to human are shown for proximal femur (red) and distal femur (blue). Published Sox9 class II sites from mouse rib chondrocytes lifted over to human are shown in light blue. Four height SNPs overlap an ATAC-seq peaks, representing putatively causal variants. Of these four variants, rs9920291, a C to T substitution, significantly alters the predicted binding of HOXD13 (UniProbe: HOXD13 C allele z-score = 0.7825; HOXD13 T allele z-score = 5.1269 (also Supplementary file 7) (B) Results of a luciferase reporter assay. A 341 bp region surrounding rs9920291 and containing either the reference or alternate allele at rs9920291 was cloned into a luciferase reporter vector and transfected into the human chondrocyte cell line T/C-28a2. Expression was normalized (to 1.0) by luciferase activity of the empty vector luciferase construct. (C) Allelic skew of two CHSY1 genic variants (rs28364839 which is exonic and rs11433 which is in the 3’UTR) in HEK-293FT cells. Notably, HEK-293FT is heterozygous for rs9920291, along with both assayed variants. (D) Expression of CHSY1 following CRISPR-Cas9-mediated deletion of regions surrounding rs9920291 (see Figure 6—figure supplement 1) in T/C-28a2 chondrocyte cell lines (N = 3 biological replicates). CHSY1 expression following deletion of an approximate 135 bp region using guides 1 and 4 is shown in purple and expression following deletion of an approximate 17 bp region using guides 6 and 8 is shown in orange. Expression values are normalized to a control CRISPR-Cas9 vector (green). (E) CHSY1 expression following HOXD13 overexpression in T/C-28a2 cells (N = 3 biological replicates) (see also Figure 6—figure supplement 1). Expression is normalized by an empty overexpression vector. ***p<0.001 by unpaired t-test.

https://doi.org/10.7554/eLife.29329.011

To further demonstrate that the region surrounding rs9920291 displays regulatory activity for CHSY1, we performed CRISPR-Cas9 targeting of a 135 bp sequence spanning a conserved portion of the regulatory element containing rs9920291 (guides sg1 and sg4), or an 17 bp sequence immediately surrounding rs9920291 (guides sg6 and sg8). Either deletion led to a significant upregulation of CHSY1 expression in vitro in T/C-28a2 cells (p=7.83×10−6 for the 135 bp deletion and p=3.64×10−4 for the 17 bp deletion; Figure 6d), confirming the regulatory element’s repressor activity. Consistent with our computational predictions that rs9920291 perturbs HOXD13 binding, we found that overexpression of HOXD13 in T/C-28a2 cells significantly upregulated CHSY1 expression (p=1.63×10−5; Figure 6e,Figure 6—figure supplement 1).

We also highlight a previously identified height-associated variant (rs4911178) at the GDF5 locus, which resides in an ATAC-seq peak active in femoral growth plates (Figure 6—figure supplement 2). We have shown in separate work that the underlying enhancer drives GDF5 expression in growth plates and that rs4911178 has allele-specific activity(Capellini et al., 2017). Thus, our approach of integrating genetic and growth plate epigenetic data has the ability to identify previously discovered cis-regulatory targets that mediate human height variation.

Motif analysis of height variants that reside within distal femur ATAC-seq peaks

To expand upon our understanding of transcription factor binding at ATAC-seq open chromatin regions, we performed a series of motif enrichment analyses using HOMER (see Materials and methods). First, we performed a de novo motif analysis on the union of distal and proximal femoral ATAC-seq peaks in mouse to identify enriched transcription-factor-binding motifs. We identified a number of highly enriched motifs relevant to chondrogenesis and other skeletal developmental processes, including Atf3 (James et al., 2006), Msx2 (Amano et al., 2008), Pbx1 (Selleri et al., 2001), and Ctcf (Jerković et al., 2017) (Figure 7a, Supplementary file 8, sheet 1). We performed a similar analysis on the liftover human ATAC-seq regions and found enrichments for a similar set of motifs (Figure 7a, Supplementary file 8, sheet 2).

Motif disruption by height GWAS variants intersecting ATAC-seq peaks.

(A) Selected enriched motifs within mouse femoral ATAC-seq open chromatin regions. These motifs were also found following liftover of the peaks to the human genome (see Supplementary file 8, sheets 1 and 2). (B) Histogram showing distribution of log-odds motif scores of CTCF-disruption for height GWAS variants within ATAC-seq regions. Negative values indicate a disruptive effect of the variant while positive values indicate a 'conservation' of the motif hit by the variant. (C) Histogram of number of variants predicted to disrupt CTCF motifs within random subsamples of variants. Dotted line indicates number of CTCF-disrupting variants observed in the height GWAS variants overlapping ATAC-seq regions (n = 78, p<0.01). (D) An example of a variant (rs12529733) that overlaps ATAC-seq peaks and is predicted to alter CTCF binding. Features in the figure are the same as described in Figure 6A.

https://doi.org/10.7554/eLife.29329.014

We next characterized the effect of height GWAS variants that intersect with these motifs within ATAC-seq peaks. Further characterization of the GWAS variants within CTCF motifs indicated a tendency toward disruptive effects (Figure 7b). Comparing the height GWAS variants intersecting ATAC-seq peaks to 500 sets of random GWAS variants revealed an enrichment for variants that are predicted to impact CTCF-binding (p<0.01, z-score: 6.01, Figure 7c) (see Materials and methods). An example of a variant within a CTCF motif is shown in Figure 7d. Additionally, we detected a strong enrichment of height GWAS variants overlapping predicted PBX1 motifs within ATAC-seq regions, when compared to height GWAS variants across the genome (adjusted p value=1.91×10−108, Supplementary file 9; see Materials and methods).

Discussion

Like many anthropometric traits, height is highly heritable and influenced by a large number of distinct genetic inputs that influence the development and/or growth of a number of tissues (Visscher et al., 2010). Of these tissues, the growth plate plays a key role in determining overall height. Despite decades of research on chondrocyte growth plate biology, we have been limited in our understanding of how height GWAS variants—most of which are non-coding—act in the growth plate to influence height. This is due in part to the paucity of epigenetic data related to the regulatory landscape of the growth plate, as these datasets have been technically challenging to generate. Here, we leveraged a new method (ATAC-seq) (Buenrostro et al., 2015, 2013) to overcome previous technical challenges and reveal the open chromatin regulatory landscape of growth plate chondrocytes. We then analyzed height GWAS data in the context of the femur growth plate epigenetic landscape to gain insight into how height GWAS variants act at the growth plate.

Our analyses revealed that ATAC-seq datasets generated from mouse growth plate chondrocytes overlap previously identified signals of chondrocyte biology in the mouse, as well as in humans. The ATAC-seq peaks reside near genes enriched for human phenotypes related to skeletal biology, such as ‘metaphyseal widening.’ These peaks also reside near genes with known roles in chondrocyte biology and include biological processes such as ‘cartilage condensation’ and ‘chondrocyte proliferation.’ We also show that our open chromatin profiles overlap quite faithfully with locations of active regulatory enhancers and Sox9-binding sites as assessed using ChIP-seq on mouse rib chondrocytes (Ohba et al., 2015). Additionally, our datasets overlap with evidence of regulatory control during chondrogenesis and human limb embryogenesis (Cotney et al., 2013; Herlofsen et al., 2013; Kundaje et al., 2015). These findings indicate that general aspects of chondrocyte regulation, in part shared between humans and mice, are captured by our ATAC-seq growth plate datasets. While we note that appropriate human datasets generated from growth plate chondrocytes are still unavailable and that there is evidence of substantial functional divergence in the genic and non-coding usage in cell types between species (Bernstein and ENCODE Project Consortium, 2012), the strong overlap between our mouse ATAC-seq data and human embryonic limb and human BMDCs suggests considerable conservation in the regulatory profiles at the growth plate.

While our datasets capture aspects of general chondrocyte biology, we also found that ATAC-seq peaks acquired on femoral tissues show enrichment in human phenotypes related to the anatomy of femur. This finding points to a fine-grained, modularized control of skeletal element-specific biology. This is also supported by our independent identification of a number of previously published long bone growth plate chondrocyte enhancers in mice and humans, as well as the growing evidence that a number of important skeletal development genes, such as Bmp5 (Guenther et al., 2015), Gdf6 (Mortlock et al., 2003; Indjeian et al., 2016), and Gdf5 (Capellini et al., 2017; Chen et al., 2016) have conserved modularized cis-regulatory architectures. As height is encoded by a number of genetic inputs, undoubtedly variants that influence general, as well as element-specific regulatory enhancers (e.g. tibia versus femur), will be important factors shaping height variation in humans.

Given that over 80% of GWAS height loci are non-coding, our study demonstrates that height GWAS variants are significantly enriched in growth plate chondrocyte regulatory regions and that the majority reside in intergenic portions of the genome. Importantly, this enrichment of height GWAS loci was not observed nearly as strongly in other chondrocyte or limb bud tissues, which suggests that our dataset is capturing functional genetic variation influencing height that is not captured by other existing datasets. No enrichment was observed for DNaseI hypersensitivity sites from mouse brain, liver, or lung, which suggests that our enrichment with height GWAS variants is tissue-specific and also not due to an artifact from the liftover process from mouse to human. These ATAC-seq peaks overlapping height GWAS variants are enriched for nearby differentially expressed growth plate genes. GWAS results for other traits/diseases that are not clearly related to growth/height did not show enrichment at the ATAC-seq peaks, demonstrating that the ATAC-seq peaks are specifically capturing growth/height biology. Finally, our motif analysis additionally revealed that a number of these human height variants may alter the binding of transcription factors with important roles in chondrogenesis and gene regulatory processes.

Combining our statistical fine-mapping and functional fine-mapping has allowed us to narrow down to a much smaller subset of putatively causal variants for downstream functional testing. As an example, one putatively causal variant (rs9920291) resides in the Chondroitin Sulfate Synthase 1 (CHSY1) locus. Human mutations in CHSY1 result in Temtamy preaxial brachydactyly syndrome, characterized by short stature, preaxial brachydactyly, and hyperphalangism of digits 1–3 (45, 46), as well as a number of additional soft and hard tissue disorders. rs9920291 is in the fine-mapped credible set, resides in regulatory elements active in human adult chondrocytes and embryonic limb tissues, and displays evidence of binding by Sox9 in mouse chondrocytes. rs9920291 additionally modifies a critical nucleotide (C to T) and disrupts the experimentally determined in vitro binding motif of HOXD13, a homeodomain transcription factor that is expressed in chondrocytes (Kuss et al., 2014; Reno et al., 2016) and has known roles in skeletal development (Johnson and Tabin, 1997). We performed extensive experimental characterization of rs9920291 at the CHSY1 locus. We demonstrated that rs9920291 displays allelic regulatory activity in human chondrocyte cell lines. We also deleted the regulatory region surrounding rs9920291 using CRISPR-Cas9 and demonstrated that these deletions alter CHSY1 expression. Thus, our results elucidate how common regulatory variation at this locus influences height by modulating CHSY1 expression. Height, given that it is generally not a cell-intrinsic phenotype, has traditionally not been conducive to functional experimentation. Our results highlight how we can pinpoint likely causal variants and mechanisms to make experimentation more tractable for a challenging phenotype such as height. To our knowledge, along with GDF5, this is the only other height GWAS locus with functional validation (Capellini et al., 2017).

A number of other variants, listed in the browsable Supplementary file 6, are available for functional testing. When more extensive height GWAS data (ideally imputed to a dense reference panel) are generated and more sophisticated fine-mapping methods are developed, additional evidence of enrichment for fine-mapped variants may be revealed using our approach. When these experiments are combined with high-throughput functional assays to test variant functionality, a more refined picture of height biology will emerge.

Materials and methods

ATAC-seq methods

Request a detailed protocol

All experiments were performed following protocols approved by the Harvard University IACUC committee. Timed matings were established between FVB mice, and pregnant females were euthanized at E15.5 in order to acquire embryonic long bone growth plates. Embryos were dissected in PBS1X on ice under a dissection scope and the proximal and distal growth zones of the right and left femur were stripped clear of soft tissues. Each proximal or distal cartilaginous end was then micro-dissected from the bony diaphysis (Figure 1) and separately pooled from a single litter, consisting on average of eight animals. Two biological replicates were collected in line with previous ATAC-seq studies (Gehrke et al., 2015). The samples were collected in micro-centrifuge tubes containing 200 ul 5% FBS/DMEM. To generate a single-cell chondrocyte suspension, each pooled sample was then subjected to 1% Collagenase II (VWR 80056–222, Radnor, Pennsylavania) digestion for 2 hr at 37°C rocking, mixing every 30 min. After placing on ice, samples were filtered using a micro-centrifuge filter set-up by gently mashing the residual tissues through the filter followed by rinsing with 5% FBS/DMEM. Samples were then spun down at 500 g at 4°C for 5 min. We next performed cell counting methods using trypan blue and a hemocytometer and performed subsequent ATAC-seq steps on those samples that had cell death rates well below 25%. On average we acquired 500,000–1,000,000 living cells per harvest. Next, cells were re-suspended in concentrations of 50,000 cells in 1x PBS. Cell samples then subjected to the ATAC-seq protocol as described in Buenrostro et al., 2015Buenrostro et al., 2015, modifying the protocol by using 2 µl of transposase per reaction. The transposase reaction product was then purified using the Omega MicroElute DNA Clean Up Kit following manufacturers protocols, eluted in 10 µl of warmed ddH20, and stored at −20°C.

Next, samples were subjected to PCR amplification and barcoding following Buenrostro et al., 2015Buenrostro et al., 2015. All primers used in this step are listed in Supplementary file 1. Ten microliters of transposed DNA were then placed in a reaction containing NEBNext High-Fidelity PCR Master Mix, ddH20, and primers. After amplification, samples were transferred to micro-centrifuge tubes and subjected to the OMEGA Bead Purification Protocol following manufacturers protocol. The samples were eluted in 30 µl of TE, run on a nanodrop, diluted to 5 ng/µl and run on a bioanalyzer. Prior to sequencing sample concentrations were determined using the KAPA Library Quantification Complete Kit (KK4824). Samples were then sent out to the Harvard University Bauer Core Facility for sequencing on one lane of the Illumina NextSeq 500 (see Supplementary file 1).

Sequencing yielded approximately 400 million reads per lane and an average of 100 million per sample. Quality control statistics are presented in Supplementary file 1. Sequenced reads were next aligned to the mouse reference mm10 genome assembly using Bowtie2 (Langmead and Salzberg, 2012). Duplicated reads were removed and significant peaks determined by using MACS2 software (version 2.1.0) (Johnson and Tabin, 1997). Peaks were assessed for reproducibility between biological replicates using the IDR statistical test (Li et al., 2011) at various statistical cut-offs, with an IDR threshold of <0.05 selected to define reproducible peaks for all subsequent analyses. Datasets for other IDR score thresholds are available upon request and when used in all analyses show similar statistical enrichments. Intersections between proximal and distal femur peaks were done with IDR-filtered narrow peak files, with the proximal/distal-common set generated by merging overlapping peaks. Annotations of the narrow peak files generated on all datasets were made using HOMER annotatePeaks.pl (http://homer.salk.edu/homer/motif/).

Raw sequencing fastq files and processed peak bed files have been deposited on NCBI GEO (GSE100585).

GREAT analyses

Request a detailed protocol

GREAT (Genomic Regions Enrichment of Annotations Tool) (McLean et al., 2010) identifies curated gene sets that display enrichment near regulatory elements. The software first identifies likely target genes of each regulatory element based on distance. GREAT then tests the enrichment of these genes localized near regulatory elements in curated gene sets that reflect pathways and molecular processes. To increase the specificity of our results, we removed peaks that were also seen in ATAC-seq from the mouse brain. ATAC-seq results were then lifted over to hg19 for analysis. We used the ‘Significant By Region-based Binomial view’; otherwise, all default input parameters and output display settings available at the GREAT website (http://great.stanford.edu/) were used. GREAT v3.0.0 was used.

Additional epigenetic datasets

Request a detailed protocol

Mouse chondrocyte data: ChIP-seq profiles Sox9 binding in mouse rib chondrocytes was obtained from Ohba et al. (2015). We also obtained ChIP-seq of RNA Pol II, H3K4me2, H3K4me3, H3K27ac, p300 and H3K27me3 in mouse chondrocytes from Ohba et al. (2015).

Human data: H3K27ac ChIP-seq of E47 human limb buds that are shared with mouse was obtained from Cotney et al., 2013 (Cotney et al., 2013). H3K27ac ChIP-seq of human adult cultured BMDCs was obtained from http://egg2.wustl.edu/roadmap/web_portal/processed_data.html (Herlofsen et al., 2013; Kundaje et al., 2015). The consolidated epigenome BroadPeak file was used.

Mouse ENCODE data: DNA-seq data were downloaded from the mouse ENCODE project (Yue et al., 2014). All samples were from Mus musculus C57BL/6 and included brain male embryo (14.5 days) (ENCSR000COE), liver male adult (8 weeks) (ENCSR000CNI), and lung male adult (8 weeks) (ENCSR000CNM). Processed NarrowPeak bed files were downloaded for each tissue and then lifted over to hg19 as described above.

Epigenetic intersection analyses

Request a detailed protocol

All intersection analyses of ATAC-seq peaks with GWAS loci, genes, or ChIP-seq peaks were performed using BEDTools v2.18 (Quinlan and Hall, 2010) and custom R scripts (R version 3.1.1) (see Source code 1). For overlaps of genes with GWAS loci, an overlap was called for a given height GWAS locus if the gene’s transcriptional boundaries were within 100 kb of any SNP in the locus (defined as r2 >0.5 to lead SNP).

For overlap of ChIP-seq datasets, any overlap (>=1 bp) was considered to represent the same feature. To evaluate the significance of overlap, a permutation-based procedure was implemented in the R package regioneR (Gel et al., 2016). We used the ‘circularRandomizeRegions’ function to generate random permutations. In this function, each chromosome is ‘circularized’, and the peaks bed file is shifted at random along the circularized chromosome and evaluated for overlap. We used 1000 permutations to calculate the empiric p value.

GoShifter analysis

Request a detailed protocol

GoShifter v0.2 was used to evaluate the statistical enrichment of the overlap between height GWAS loci and ATAC-seq peaks (or other epigenetic peaks) (Trynka et al., 2015). The European (EUR) subset of 1000 Genomes Phase 3 (Auton et al., 2015) was used to identify proxy SNPs and calculate LD. An LD cutoff of r2 >0.5 to the index SNP was used. A window size of 100 kb was used to find LD SNPs (‘window’ option). All other parameters were set to the software default. 100,000 permutations were run. GoShifter analyses for additional traits/diseases were performed in the same fashion.

Generation of matched GWAS loci

Request a detailed protocol

SNPSNAP (Pers et al., 2015) was used to generate matched loci for the height GWAS loci. The European (EUR) subset of 1000 Genomes Phase 3 (Auton et al., 2015) was used to identify proxy SNPs and calculate LD. Loci were defined as variants with r2 > 0.5 to the index SNP. Loci were matched using default criteria: MAF ± 5%, gene density ± 50%, distance to nearest gene ±50% and LD buddies ± 50%. SNPSNAP was able to generate matched gene sets for 676 out of the 697 loci from Wood et al. (2014). 1000 matched sets of random GWAS loci were generated.

Differentially expressed growth plate genes

Request a detailed protocol

A list of 427 differentially expressed mouse genes was downloaded from Lui et al. (2012). These 427 differentially expressed genes were determined in the Lui et al. (2012) paper based on meeting two of three criteria: (1) spatial regulation in the growth plate, (2) temporal regulation during growth plate formation, and/or (3) specific expression in the growth plate (as compared to other tissues such as the lung, kidney, heart). Orthologous human genes and their respective positions were identified using Ensembl Biomart (Kinsella et al., 2011). The ‘Ensembl Genes 86’ database was used.

eQTL overlap analysis

Request a detailed protocol

Using the set of 427 differentially expressed growth plate genes (see above), we identified 157 genes that were also genes for cis-eQTLs in whole blood (FDR < 0.05) as identified by Westra et al. (2013). For each of these 157 eQTL genes, the corresponding eQTL variant with the strongest p value was used. Using this set of eQTL variants, GoShifter was run as described above to test for enrichment at ATAC-seq peaks.

GTEx data

Request a detailed protocol

The Genotype-Tissue Expression (GTEx) Project (Lonsdale J, GTEx Consortium, 2013) was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from the GTEx Portal on 10/27/2017.

Motif analysis

Request a detailed protocol

ATAC-seq peak centers from proximal/distal tissues were aggregated and padded outwards to a fixed length of 200 bp (based on the distribution of ATAC-seq peak sizes). HOMER (Jerković et al., 2017) de novo motif analysis was performed for the sequence set utilizing a 10x random shuffling as a background set. De novo motifs were compared to a vertebrate motif library provided by HOMER, which includes the JASPAR database (Visscher et al., 2010), with matches scored using Pearson’s correlation coefficient of vectorized motif matrices, with neutral frequencies (0.25) substituted for non-overlapping positions.

Lifted over human ATAC-seq peaks were then intersected with height GWAS variants from Wood et al. (2014). Variants were characterized for selected motifs which appeared in both the mm10 and hg19 de novo motifs using motifbreakR (Coetzee et al., 2015). Motif matches were scored using log-odds scoring (with 0.25 A/C/G/T background frequencies) for both reference and mutated sequences. Effects on position-weighted matrix (PWM) scoring were calculated as differences between the log-odds score of reference and mutant motif hits, with negative differences representing a disrupted motif. To assess enrichment for CTCF-disrupting variants in the ATAC-intersected variant set, 500 randomly sampled subsets (n = 928) of the height GWAS variants (excluding those intersecting ATAC-seq peaks) were generated and counted for CTCF-disrupting variants. Counts were standardized and statistical significance was assessed using a continuous distribution function of the standard normal distribution.

To assess enrichment of predicted transcription factor binding near variants intersecting ATAC-seq peaks, a set of hg19 DNA sequences was first constructed using 30 bp windows centered on each variant, generating an ATAC-intersected subset and whole GWAS superset. A utility from HOMER (findmotifs.pl) was then used to scan this sequence set using its provided vertebrate motif library to search for motif hits. Motifs from JASPAR 2016 were filtered (386 motifs) and counted within each of the two sets. Hypergeometric tests were performed for all motifs using Benjamini-Hochberg FDR correction (Benjamini and Hochberg, 1995).

Expression plasmids

Request a detailed protocol

CRISPR-Cas9 PX458 vector was obtained from Addgene (Cambridge, Massachusetts). pGL4.23 luciferase and pGL4.74 renilla vectors were acquired from Promega (Madison, Wisconsin). Plasmids encoding HOXD13-FLAG and HOXD13-HA were generously gifted from Dr. Vincenzo Zappavigna (University of Modena and Reggio Emilia) (Caronia et al., 2003).

Cell lines and culture conditions

Request a detailed protocol

Human embryonic kidney (HEK-293FT) cells were acquired from Dr. Pardis Sabeti (Harvard University and Broad Institute). T/C-28a2 human chondrocyte cells were acquired from Dr. Li Zeng (Tufts University) courtesy of Dr. Mary Goldring (The Hospital for Special Surgery). Both cell lines were cultured at 5% CO2 at 37°C in Dulbecco's Modified Eagle's Medium (DMEM), 10% fetal bovine serum (FBS), and 1% penicillin-streptomycin (P/S). Media was replaced every 2–3 days and the cells were sub-cultured every 5 days.

CRISPR-Cas9 targeting at the CHSY1 regulatory element

Request a detailed protocol

All guide RNAs surrounding the human CHSY1 regulatory element or rs9920291 were designed using MIT CRISPR Tools (http://crispr.mit.edu), synthesized by Integrated DNA Technologies, Inc (Coralville, Iowa), and cloned into the PX458 vector following published protocols (Ran et al., 2013). The sequence of all sgRNAs are listed in Supplementary file 10 and their locations with respect to the targeted CHSY1 regulatory element and rs9920291 are found in Figure 6—figure supplement 1. Guide RNAs were initially tested for ability to induce efficient deletions of the human element in cultured HEK-293FT cells. After two days of culture at 37°C, transfected HEK-293FT cells were examined under a GFP-microscope to verify successful transfection and GFP expression. DNA was then extracted using E.Z.N.A Tissue DNA Kit, and the CHSY1 regulatory element region was amplified using PCR and primers surrounding the guide RNA design sites (Primers: ‘rs9920291 Forward’ and ‘rs9920291 Reverse’, Supplementary file 10). Amplification products were isolated from 1% agarose gels (E.Z.N.A Gel Extraction Kit) and Sanger sequenced to verify deletions of various sizes (see below).

After confirmation that sgRNAs worked in HEK-293FT cells, we performed in vitro deletion in the T/C-28a2 chondrocyte cell line (Kokenyesi et al., 2000; Finger et al., 2003). The same sgRNAs were used as above. T/C-28a2 cells were maintained in DMEM (Gibco, Gaithersburg, Maryland) supplied with 10% FBS (Gibco) and 1% Pen/Strep (0.025%), and seeded in a six-well plate 1 day prior to transfection. After culturing at 37°C, we scanned cells under a GFP-microscope to verify successful transfection efficiency (i.e.,>70% of cells) and GFP expression (N = 3). DNA was then extracted using E.Z.N.A Tissue DNA Kit, and the CHSY1 regulatory element region was amplified using PCR primers surrounding the guide RNA design sites (MiSeq Primers Forward and MiSeq Primers Reverse, Supplementary file 10). Amplification products were isolated from 1% agarose gels (E.Z.N.A Gel Extraction Kit). Mi-Seq sequencing at the MGH CCIB DNA Core was used to verify successful targeting of the larger CHSY1 regulatory region (hg19: chr15:101,749,730-101,749,865) with a modification efficiency of approximately 10% and smaller CHSY1 regulatory region around rs9920291 (hg19: chr15:101,749,797–101,749,813) with a modification efficiency of approximately 8.5%. RNA was also extracted from control and CRISPR-Cas9 targeted T/C-28a2 cells. RNA was DNase-treated and converted to cDNA for qPCR (see below).

cDNA synthesis and qRT-PCR analysis

Request a detailed protocol

Total RNA was extracted from HEK-293FT and T/C-28a2 cells and prepared using the Trizol Reagent (Thermo Fisher Scientific, Springfield Township, New Jersey) and Direct-zolTM RNA Miniprep kit (ZYMO). Three micrograms of total RNA were used to synthesize first-strand cDNA using SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific). qRT-PCR analysis was then performed with specific primers and Applied Biosystems Power SYBR master mix (Thermo Fisher Scientific) with GAPDH as a reference gene. Primers used for qRT-PCR are listed in Supplementary file 10.

Cloning, construct preparation, in vitro luciferase reporter and overexpression assays

Request a detailed protocol

Prior to in vitro transfection of CHSY1 regulatory sequences into HEK-293FT and T/C-28a2 cells (see below), primers containing KpnI/HindIII linker sequences were used to amplify the putative regulatory element surrounding rs9920291. (see above and Supplementary file 10). The PCR protocol was initiated at 98°C for 30 s, followed by 34 cycles at 98°C for 20 s, 60°C for 20 s, and 72°C for 30 s, and then followed by a 72°C for 5 min final extension. Coriell (Camden, New Jersey) sample NA12286 was used to amplify the region containing the reference allele of rs9920291 and NA19119 was used to amplify the region containing the variant allele of rs9920291. Each amplicon was then ligated into a KpnI/HindIII containing pGL4.23 firefly luciferase vector using NEBuffer 2.1 (NEB B7202S) for digestion and T4 DNA Ligase (NEB M0202S) for ligation according to manufacturer's protocols (New England Biolabs, Beverly, Massachusetts). Ligates for constructs containing inserted or non-inserted (‘empty’) sequence were next transformed into DH10B cells, and streaked on ampicillin plates. Single colonies were then picked, PCR screened, sequenced, and purified using E.Z.N.A. Endo-Free Plasmid DNA Maxi Kits (Omega D6926-03). The PCR protocol used for screening is listed above. Sanger sequencing was used to confirm the orientation as well as sequence identity of each insert as well as the entire pGL4.23 firefly vector. pGL4.74 renilla vector was similarly transformed, amplified, purified and sequence verified. Prior to HOXD13 over-expression experiments, HOXD13 constructs were transformed into DH10B cells, streaked on ampicillin plates, individual colonies were grown and purified via Endo-Free Plasmid DNA Maxi Kits (Omega D6926-03). Purified plasmids underwent diagnostic digest and sequence verification.

HEK-293FT cells and T/C-28a2 chondrocytes were passaged using standard conditions (see above) prior to reporter gene transfection and HOXD13 over-expression experiments. Prior to transfection, cells were first cultured in DMEM with 10% FBS for 24 hr in 96-well dishes at a seeding density of 3 × 104 cells/well. The volume of media per well was then brought to 100 µl. For luciferase reporter experiments, cells were then transiently transfected with 100 ng firefly luciferase reporter vector or empty pGL4.23 luciferase vector, along with 5ng pGL4.74 renilla vector. Transfections were performed using Lipofectamine 2000 (Invitrogen 11668–019) at a Lipofectamine:DNA ratio of 4:1 according to manufacturer instructions. Luciferase activity was measured 48 hr after transfection using the 96-welll Dual-Glo Luciferase Assay System (Promega E2940) following manufacturer protocols on a SpectraMax L Microplate Reader (Molecular Devices; Cat# SpectraMax L Config) with a 1 min dark adapt, 5 s integration, and max range settings. For cell culture experiments, to compare expression between the reference and alternate allele CHSY1 constructs (Figure 6b), we performed at least four independent transfection experiments at each concentration containing eight technical replicates (i.e. individual wells of a 96-well plate) per construct. We first normalized each firefly luciferase value per well by its corresponding renilla luciferase value per well, and then compared the mean expression of the reference allele construct (normalized by empty vector) to that of the mean expression of the alternate allele construct (normalized by empty vector) using a two-sample directional Student's t-test. For HOXD13 overexpression experiments, HEK-293FT and T/C-28a2 cells were transfected as above but with 0, 2 and 4 µg of HOXD13 expression or control constructs (N = 3 biological replicates per condition), and after 24, 48, and 72 hr RNA was collected for qRT-PCR experiments (see above).

CHSY1 allelic skew analyses

Request a detailed protocol

Both HEK-293FT cells and T/C-28a2 chondrocytes were screened using Sanger sequencing to identify genotype rs9920291. Sanger sequencing results confirmed preliminary findings using the HEK-293 browser (http://hek293genome.org/v2/) that HEK-293FT cells were heterozygous at rs9920291, unlike T/C-28a2 cells which were homozygous for the reference C allele. Thus, only HEK-293FT cells could be used for allele-specific expression. Next, using the HEK-293 browser with follow-up verification using Sanger sequencing, we identified two heterozygous transcript variants (exon 3 coding variant [rs28364839] and 3’UTR variant [rs11433]).

For allelic skew analyses on untreated HEK-293FT samples, five independent biological replicates were grown and harvested. RNA was isolated separately from each sample using the Trizol Reagent (15596–026, Ambion by Life Technologies, Norwalk, Connecticut) and the RNA Clean and ConcentratorTM-5 Kit (supplied with DNase I, Zymo). Samples were then run on a bioanalyzer to achieve RNA Integrity Numbers greater than 8. These RNA samples were then reverse transcribed using a SuperScript IV First Strand cDNA Synthesis Reaction kit (18090010, Life Technologies) following manufacturer’s recommendations.

cDNA samples were then sent to EpigenDx for allelic skew analysis assay design and execution. SNPs in the coding regions of CHSY1 (rs11433 and rs28364839), were validated by EpigenDx (Hopkinton, Massachusetts). Pyrosequencing for SNP genotyping (PSQ H96A, Qiagen Pyrosequencing) is a real-time sequencing-based DNA analysis that quantitatively determines the genotypes of single or multiple mutations in a single reaction. Briefly, 1 ng of sample cDNA was used for PCR amplification. PCR was performed using 10X PCR buffer (Qiagen Inc., Maryland) at 3.0 mM MgCl2, 200M of each dNTP, 0.2 µM each of the forward and reverse primers (available through EpigenDx), and 0.75 U of HotStar DNA polymerase (Qiagen Inc.), per 30 µl reaction. PCR cycling conditions were: 94°C 15 min; 45 cycles of 94°C 30 s; 60°C 30 s; 72°C 30 s; 72°C 5 min. One of the PCR primer pairs was biotinylated to convert the PCR product to single-stranded DNA sequencing templates using streptavidin beads and the PyroMark Q96 Vacuum Workstation. Ten microliters of the PCR products were bound to streptavidin beads and the single strand containing the biotinylated primer was isolated and combined with a specific sequencing primer (available through EpigenDx). The primed single stranded DNA was sequenced using the Pyrosequencing PSQ96 HS System (Qiagen Pyrosequencing) following the manufacturer’s instructions (Qiagen Pyrosequencing). The genotypes of each sample were analyzed using Q96 software AQ module (Qiagen Pyrosequencing). Pyrosequencing results for each CHSY1 exonic or 3’UTR SNP were used to calculate the allelic ratios in the heterozygous state.

Statistical analysis

Request a detailed protocol

For regulatory element transfection experiments, CRISPR-Cas9 deletion experiments, and HOXD13 over-expression experiments, the means ± SEM of multiple independent measurements were calculated. The unpaired two-tailed Student’s t test was used to determine the significance of differences between means. P values smaller than 0.05 (*p<0.05, **p<0.01, ***p<0.001) were considered to be statistical significant.

UniProbe analysis

Request a detailed protocol

To find upstream transcription factors (TF) predicted to bind at rs9920291 near CHYS1 we used UniProbe (Hume et al., 2015). UniProbe (http://the_brain.bwh.harvard.edu/uniprobe/about.php) is built on experimental measurements of binding affinities between large numbers of expressed TFs and all possible 8-mer target oligonucleotides (Hume et al., 2015). For two 15 base-pair sequences, each centered at rs9920291 but different only at the C/T variant (Reference [REF] TGCTAAACAAAAATA and Alternate [ALT] TGCTAAATAAAAATA), we imported each sequence separately into the ‘Search for TF Binding Sites’ window on the UniProbe browser (http://the_brain.bwh.harvard.edu/uniprobe/) with the following conditions: Score Threshold: 0.4; Species: Homo sapiens. We acquired 18 TF predictions for REF and 61 TF predictions for ALT. We next imported the results into Supplementary file 7 (sheet 1: REF TGCTAAACAAAAATA and sheet 2: ALT TGCTAAATAAAAATA) and performed an intersection (sheet 3: REF (red) vs. ALT (blue)) to identify those TF-binding preferences that were reduced/gained between the two sequences. For a TF-binding site to be considered reduced or gained, the site must exhibit an enrichment score change of greater than at least 0.10, which corresponded to significant changes in z-score (e.g., Figure 6a). To identify the change in HOXD13 binding affinity, we re-ran UniProbe on the REF sequence but at a lower enrichment score of E >= 0.2. This yielded 341 TF predictions (sheet 4), which were then compared to the 61 TF predictions generated at E >= 0.4 (sheet 5). Finally, to generate the z-score, we downloaded HOXD13 experiment information from UniProbe (http://the_brain.bwh.harvard.edu/uniprobe/browse.php) and identified the corresponding z-score for each 8-mer.

Lists of potential upstream regulators exhibiting marked reduction/gain in binding affinity were then screened using expression and phenotypic data to identify those TFs also expressed or functionally required in growth plates. To carry out this analysis, we used data in VisiGene (http://genome.ucsc.edu/cgi-bin/hgVisigene), Eurexpress (http://www.eurexpress.org/ee/), Genepaint (http://www.genepaint.org/Frameset.html), and the Mouse Genome Informatics expression and phenotypic databases (http://www.informatics.jax.org). HOXD13 met the above criteria. HOXD13 had enrichment scores below the recommended UniProbe threshold of 0.4 for the REF sequence but much higher (above 0.4) for the ALT sequence, and this corresponded to a significant change in z-score (Figure 6a).

Data availability

The following data sets were generated

References

    1. Benjamini Y
    2. Hochberg Y
    (1995)
    Controlling the false discovery rate: a practical and powerful approach to multiple testing
    Journal of the Royal Statistical Society 57:289–300.
    1. Buenrostro JD
    2. Wu B
    3. Chang HY
    4. Greenleaf WJ
    (2015)
    ATAC-seq: A method for assaying chromatin accessibility genome-wide
    Current Protocols in Molecular Biology 109:21.29.1–21.2921.
    1. Day FR
    2. Thompson DJ
    3. Helgason H
    4. Chasman DI
    5. Finucane H
    6. Sulem P
    7. Ruth KS
    8. Whalen S
    9. Sarkar AK
    10. Albrecht E
    11. Altmaier E
    12. Amini M
    13. Barbieri CM
    14. Boutin T
    15. Campbell A
    16. Demerath E
    17. Giri A
    18. He C
    19. Hottenga JJ
    20. Karlsson R
    21. Kolcic I
    22. Loh PR
    23. Lunetta KL
    24. Mangino M
    25. Marco B
    26. McMahon G
    27. Medland SE
    28. Nolte IM
    29. Noordam R
    30. Nutile T
    31. Paternoster L
    32. Perjakova N
    33. Porcu E
    34. Rose LM
    35. Schraut KE
    36. Segrè AV
    37. Smith AV
    38. Stolk L
    39. Teumer A
    40. Andrulis IL
    41. Bandinelli S
    42. Beckmann MW
    43. Benitez J
    44. Bergmann S
    45. Bochud M
    46. Boerwinkle E
    47. Bojesen SE
    48. Bolla MK
    49. Brand JS
    50. Brauch H
    51. Brenner H
    52. Broer L
    53. Brüning T
    54. Buring JE
    55. Campbell H
    56. Catamo E
    57. Chanock S
    58. Chenevix-Trench G
    59. Corre T
    60. Couch FJ
    61. Cousminer DL
    62. Cox A
    63. Crisponi L
    64. Czene K
    65. Davey Smith G
    66. de Geus E
    67. de Mutsert R
    68. De Vivo I
    69. Dennis J
    70. Devilee P
    71. Dos-Santos-Silva I
    72. Dunning AM
    73. Eriksson JG
    74. Fasching PA
    75. Fernández-Rhodes L
    76. Ferrucci L
    77. Flesch-Janys D
    78. Franke L
    79. Gabrielson M
    80. Gandin I
    81. Giles GG
    82. Grallert H
    83. Gudbjartsson DF
    84. Guénel P
    85. Hall P
    86. Hallberg E
    87. Hamann U
    88. Harris TB
    89. Hartman CA
    90. Heiss G
    91. Hooning MJ
    92. Hopper JL
    93. Hu F
    94. Hunter DJ
    95. Ikram MA
    96. Im HK
    97. Järvelin MR
    98. Joshi PK
    99. Karasik D
    100. Kellis M
    101. Kutalik Z
    102. LaChance G
    103. Lambrechts D
    104. Langenberg C
    105. Launer LJ
    106. Laven JSE
    107. Lenarduzzi S
    108. Li J
    109. Lind PA
    110. Lindstrom S
    111. Liu Y
    112. Luan J
    113. Mägi R
    114. Mannermaa A
    115. Mbarek H
    116. McCarthy MI
    117. Meisinger C
    118. Meitinger T
    119. Menni C
    120. Metspalu A
    121. Michailidou K
    122. Milani L
    123. Milne RL
    124. Montgomery GW
    125. Mulligan AM
    126. Nalls MA
    127. Navarro P
    128. Nevanlinna H
    129. Nyholt DR
    130. Oldehinkel AJ
    131. O'Mara TA
    132. Padmanabhan S
    133. Palotie A
    134. Pedersen N
    135. Peters A
    136. Peto J
    137. Pharoah PDP
    138. Pouta A
    139. Radice P
    140. Rahman I
    141. Ring SM
    142. Robino A
    143. Rosendaal FR
    144. Rudan I
    145. Rueedi R
    146. Ruggiero D
    147. Sala CF
    148. Schmidt MK
    149. Scott RA
    150. Shah M
    151. Sorice R
    152. Southey MC
    153. Sovio U
    154. Stampfer M
    155. Steri M
    156. Strauch K
    157. Tanaka T
    158. Tikkanen E
    159. Timpson NJ
    160. Traglia M
    161. Truong T
    162. Tyrer JP
    163. Uitterlinden AG
    164. Edwards DRV
    165. Vitart V
    166. Völker U
    167. Vollenweider P
    168. Wang Q
    169. Widen E
    170. van Dijk KW
    171. Willemsen G
    172. Winqvist R
    173. Wolffenbuttel BHR
    174. Zhao JH
    175. Zoledziewska M
    176. Zygmunt M
    177. Alizadeh BZ
    178. Boomsma DI
    179. Ciullo M
    180. Cucca F
    181. Esko T
    182. Franceschini N
    183. Gieger C
    184. Gudnason V
    185. Hayward C
    186. Kraft P
    187. Lawlor DA
    188. Magnusson PKE
    189. Martin NG
    190. Mook-Kanamori DO
    191. Nohr EA
    192. Polasek O
    193. Porteous D
    194. Price AL
    195. Ridker PM
    196. Snieder H
    197. Spector TD
    198. Stöckl D
    199. Toniolo D
    200. Ulivi S
    201. Visser JA
    202. Völzke H
    203. Wareham NJ
    204. Wilson JF
    205. Spurdle AB
    206. Thorsteindottir U
    207. Pollard KS
    208. Easton DF
    209. Tung JY
    210. Chang-Claude J
    211. Hinds D
    212. Murray A
    213. Murabito JM
    214. Stefansson K
    215. Ong KK
    216. Perry JRB
    217. LifeLines Cohort Study
    218. InterAct Consortium
    219. kConFab/AOCS Investigators
    220. Endometrial Cancer Association Consortium
    221. Ovarian Cancer Association Consortium
    222. PRACTICAL consortium
    (2017) Genomic analyses identify hundreds of variants associated with age at menarche and support a role for puberty timing in cancer risk
    Nature Genetics 49:834–841.
    https://doi.org/10.1038/ng.3841
    1. Galton F
    (1886) Regression Towards mediocrity in hereditary stature
    The Journal of the Anthropological Institute of Great Britain and Ireland 15:246–263.
    https://doi.org/10.2307/2841583
    1. Lango Allen H
    2. Estrada K
    3. Lettre G
    4. Berndt SI
    5. Weedon MN
    6. Rivadeneira F
    7. Willer CJ
    8. Jackson AU
    9. Vedantam S
    10. Raychaudhuri S
    11. Ferreira T
    12. Wood AR
    13. Weyant RJ
    14. Segrè AV
    15. Speliotes EK
    16. Wheeler E
    17. Soranzo N
    18. Park JH
    19. Yang J
    20. Gudbjartsson D
    21. Heard-Costa NL
    22. Randall JC
    23. Qi L
    24. Vernon Smith A
    25. Mägi R
    26. Pastinen T
    27. Liang L
    28. Heid IM
    29. Luan J
    30. Thorleifsson G
    31. Winkler TW
    32. Goddard ME
    33. Sin Lo K
    34. Palmer C
    35. Workalemahu T
    36. Aulchenko YS
    37. Johansson A
    38. Zillikens MC
    39. Feitosa MF
    40. Esko T
    41. Johnson T
    42. Ketkar S
    43. Kraft P
    44. Mangino M
    45. Prokopenko I
    46. Absher D
    47. Albrecht E
    48. Ernst F
    49. Glazer NL
    50. Hayward C
    51. Hottenga JJ
    52. Jacobs KB
    53. Knowles JW
    54. Kutalik Z
    55. Monda KL
    56. Polasek O
    57. Preuss M
    58. Rayner NW
    59. Robertson NR
    60. Steinthorsdottir V
    61. Tyrer JP
    62. Voight BF
    63. Wiklund F
    64. Xu J
    65. Zhao JH
    66. Nyholt DR
    67. Pellikka N
    68. Perola M
    69. Perry JR
    70. Surakka I
    71. Tammesoo ML
    72. Altmaier EL
    73. Amin N
    74. Aspelund T
    75. Bhangale T
    76. Boucher G
    77. Chasman DI
    78. Chen C
    79. Coin L
    80. Cooper MN
    81. Dixon AL
    82. Gibson Q
    83. Grundberg E
    84. Hao K
    85. Juhani Junttila M
    86. Kaplan LM
    87. Kettunen J
    88. König IR
    89. Kwan T
    90. Lawrence RW
    91. Levinson DF
    92. Lorentzon M
    93. McKnight B
    94. Morris AP
    95. Müller M
    96. Suh Ngwa J
    97. Purcell S
    98. Rafelt S
    99. Salem RM
    100. Salvi E
    101. Sanna S
    102. Shi J
    103. Sovio U
    104. Thompson JR
    105. Turchin MC
    106. Vandenput L
    107. Verlaan DJ
    108. Vitart V
    109. White CC
    110. Ziegler A
    111. Almgren P
    112. Balmforth AJ
    113. Campbell H
    114. Citterio L
    115. De Grandi A
    116. Dominiczak A
    117. Duan J
    118. Elliott P
    119. Elosua R
    120. Eriksson JG
    121. Freimer NB
    122. Geus EJ
    123. Glorioso N
    124. Haiqing S
    125. Hartikainen AL
    126. Havulinna AS
    127. Hicks AA
    128. Hui J
    129. Igl W
    130. Illig T
    131. Jula A
    132. Kajantie E
    133. Kilpeläinen TO
    134. Koiranen M
    135. Kolcic I
    136. Koskinen S
    137. Kovacs P
    138. Laitinen J
    139. Liu J
    140. Lokki ML
    141. Marusic A
    142. Maschio A
    143. Meitinger T
    144. Mulas A
    145. Paré G
    146. Parker AN
    147. Peden JF
    148. Petersmann A
    149. Pichler I
    150. Pietiläinen KH
    151. Pouta A
    152. Ridderstråle M
    153. Rotter JI
    154. Sambrook JG
    155. Sanders AR
    156. Schmidt CO
    157. Sinisalo J
    158. Smit JH
    159. Stringham HM
    160. Bragi Walters G
    161. Widen E
    162. Wild SH
    163. Willemsen G
    164. Zagato L
    165. Zgaga L
    166. Zitting P
    167. Alavere H
    168. Farrall M
    169. McArdle WL
    170. Nelis M
    171. Peters MJ
    172. Ripatti S
    173. van Meurs JB
    174. Aben KK
    175. Ardlie KG
    176. Beckmann JS
    177. Beilby JP
    178. Bergman RN
    179. Bergmann S
    180. Collins FS
    181. Cusi D
    182. den Heijer M
    183. Eiriksdottir G
    184. Gejman PV
    185. Hall AS
    186. Hamsten A
    187. Huikuri HV
    188. Iribarren C
    189. Kähönen M
    190. Kaprio J
    191. Kathiresan S
    192. Kiemeney L
    193. Kocher T
    194. Launer LJ
    195. Lehtimäki T
    196. Melander O
    197. Mosley TH
    198. Musk AW
    199. Nieminen MS
    200. O'Donnell CJ
    201. Ohlsson C
    202. Oostra B
    203. Palmer LJ
    204. Raitakari O
    205. Ridker PM
    206. Rioux JD
    207. Rissanen A
    208. Rivolta C
    209. Schunkert H
    210. Shuldiner AR
    211. Siscovick DS
    212. Stumvoll M
    213. Tönjes A
    214. Tuomilehto J
    215. van Ommen GJ
    216. Viikari J
    217. Heath AC
    218. Martin NG
    219. Montgomery GW
    220. Province MA
    221. Kayser M
    222. Arnold AM
    223. Atwood LD
    224. Boerwinkle E
    225. Chanock SJ
    226. Deloukas P
    227. Gieger C
    228. Grönberg H
    229. Hall P
    230. Hattersley AT
    231. Hengstenberg C
    232. Hoffman W
    233. Lathrop GM
    234. Salomaa V
    235. Schreiber S
    236. Uda M
    237. Waterworth D
    238. Wright AF
    239. Assimes TL
    240. Barroso I
    241. Hofman A
    242. Mohlke KL
    243. Boomsma DI
    244. Caulfield MJ
    245. Cupples LA
    246. Erdmann J
    247. Fox CS
    248. Gudnason V
    249. Gyllensten U
    250. Harris TB
    251. Hayes RB
    252. Jarvelin MR
    253. Mooser V
    254. Munroe PB
    255. Ouwehand WH
    256. Penninx BW
    257. Pramstaller PP
    258. Quertermous T
    259. Rudan I
    260. Samani NJ
    261. Spector TD
    262. Völzke H
    263. Watkins H
    264. Wilson JF
    265. Groop LC
    266. Haritunians T
    267. Hu FB
    268. Kaplan RC
    269. Metspalu A
    270. North KE
    271. Schlessinger D
    272. Wareham NJ
    273. Hunter DJ
    274. O'Connell JR
    275. Strachan DP
    276. Wichmann HE
    277. Borecki IB
    278. van Duijn CM
    279. Schadt EE
    280. Thorsteinsdottir U
    281. Peltonen L
    282. Uitterlinden AG
    283. Visscher PM
    284. Chatterjee N
    285. Loos RJ
    286. Boehnke M
    287. McCarthy MI
    288. Ingelsson E
    289. Lindgren CM
    290. Abecasis GR
    291. Stefansson K
    292. Frayling TM
    293. Hirschhorn JN
    (2010) Hundreds of variants clustered in genomic loci and biological pathways affect human height
    Nature 467:832–838.
    https://doi.org/10.1038/nature09410
    1. Locke AE
    2. Kahali B
    3. Berndt SI
    4. Justice AE
    5. Pers TH
    6. Day FR
    7. Powell C
    8. Vedantam S
    9. Buchkovich ML
    10. Yang J
    11. Croteau-Chonka DC
    12. Esko T
    13. Fall T
    14. Ferreira T
    15. Gustafsson S
    16. Kutalik Z
    17. Luan J
    18. Mägi R
    19. Randall JC
    20. Winkler TW
    21. Wood AR
    22. Workalemahu T
    23. Faul JD
    24. Smith JA
    25. Zhao JH
    26. Zhao W
    27. Chen J
    28. Fehrmann R
    29. Hedman ÅK
    30. Karjalainen J
    31. Schmidt EM
    32. Absher D
    33. Amin N
    34. Anderson D
    35. Beekman M
    36. Bolton JL
    37. Bragg-Gresham JL
    38. Buyske S
    39. Demirkan A
    40. Deng G
    41. Ehret GB
    42. Feenstra B
    43. Feitosa MF
    44. Fischer K
    45. Goel A
    46. Gong J
    47. Jackson AU
    48. Kanoni S
    49. Kleber ME
    50. Kristiansson K
    51. Lim U
    52. Lotay V
    53. Mangino M
    54. Leach IM
    55. Medina-Gomez C
    56. Medland SE
    57. Nalls MA
    58. Palmer CD
    59. Pasko D
    60. Pechlivanis S
    61. Peters MJ
    62. Prokopenko I
    63. Shungin D
    64. Stančáková A
    65. Strawbridge RJ
    66. Sung YJ
    67. Tanaka T
    68. Teumer A
    69. Trompet S
    70. van der Laan SW
    71. van Setten J
    72. Van Vliet-Ostaptchouk JV
    73. Wang Z
    74. Yengo L
    75. Zhang W
    76. Isaacs A
    77. Albrecht E
    78. Ärnlöv J
    79. Arscott GM
    80. Attwood AP
    81. Bandinelli S
    82. Barrett A
    83. Bas IN
    84. Bellis C
    85. Bennett AJ
    86. Berne C
    87. Blagieva R
    88. Blüher M
    89. Böhringer S
    90. Bonnycastle LL
    91. Böttcher Y
    92. Boyd HA
    93. Bruinenberg M
    94. Caspersen IH
    95. Chen YI
    96. Clarke R
    97. Daw EW
    98. de Craen AJM
    99. Delgado G
    100. Dimitriou M
    101. Doney ASF
    102. Eklund N
    103. Estrada K
    104. Eury E
    105. Folkersen L
    106. Fraser RM
    107. Garcia ME
    108. Geller F
    109. Giedraitis V
    110. Gigante B
    111. Go AS
    112. Golay A
    113. Goodall AH
    114. Gordon SD
    115. Gorski M
    116. Grabe HJ
    117. Grallert H
    118. Grammer TB
    119. Gräßler J
    120. Grönberg H
    121. Groves CJ
    122. Gusto G
    123. Haessler J
    124. Hall P
    125. Haller T
    126. Hallmans G
    127. Hartman CA
    128. Hassinen M
    129. Hayward C
    130. Heard-Costa NL
    131. Helmer Q
    132. Hengstenberg C
    133. Holmen O
    134. Hottenga JJ
    135. James AL
    136. Jeff JM
    137. Johansson Å
    138. Jolley J
    139. Juliusdottir T
    140. Kinnunen L
    141. Koenig W
    142. Koskenvuo M
    143. Kratzer W
    144. Laitinen J
    145. Lamina C
    146. Leander K
    147. Lee NR
    148. Lichtner P
    149. Lind L
    150. Lindström J
    151. Lo KS
    152. Lobbens S
    153. Lorbeer R
    154. Lu Y
    155. Mach F
    156. Magnusson PKE
    157. Mahajan A
    158. McArdle WL
    159. McLachlan S
    160. Menni C
    161. Merger S
    162. Mihailov E
    163. Milani L
    164. Moayyeri A
    165. Monda KL
    166. Morken MA
    167. Mulas A
    168. Müller G
    169. Müller-Nurasyid M
    170. Musk AW
    171. Nagaraja R
    172. Nöthen MM
    173. Nolte IM
    174. Pilz S
    175. Rayner NW
    176. Renstrom F
    177. Rettig R
    178. Ried JS
    179. Ripke S
    180. Robertson NR
    181. Rose LM
    182. Sanna S
    183. Scharnagl H
    184. Scholtens S
    185. Schumacher FR
    186. Scott WR
    187. Seufferlein T
    188. Shi J
    189. Smith AV
    190. Smolonska J
    191. Stanton AV
    192. Steinthorsdottir V
    193. Stirrups K
    194. Stringham HM
    195. Sundström J
    196. Swertz MA
    197. Swift AJ
    198. Syvänen AC
    199. Tan ST
    200. Tayo BO
    201. Thorand B
    202. Thorleifsson G
    203. Tyrer JP
    204. Uh HW
    205. Vandenput L
    206. Verhulst FC
    207. Vermeulen SH
    208. Verweij N
    209. Vonk JM
    210. Waite LL
    211. Warren HR
    212. Waterworth D
    213. Weedon MN
    214. Wilkens LR
    215. Willenborg C
    216. Wilsgaard T
    217. Wojczynski MK
    218. Wong A
    219. Wright AF
    220. Zhang Q
    221. Brennan EP
    222. Choi M
    223. Dastani Z
    224. Drong AW
    225. Eriksson P
    226. Franco-Cereceda A
    227. Gådin JR
    228. Gharavi AG
    229. Goddard ME
    230. Handsaker RE
    231. Huang J
    232. Karpe F
    233. Kathiresan S
    234. Keildson S
    235. Kiryluk K
    236. Kubo M
    237. Lee JY
    238. Liang L
    239. Lifton RP
    240. Ma B
    241. McCarroll SA
    242. McKnight AJ
    243. Min JL
    244. Moffatt MF
    245. Montgomery GW
    246. Murabito JM
    247. Nicholson G
    248. Nyholt DR
    249. Okada Y
    250. Perry JRB
    251. Dorajoo R
    252. Reinmaa E
    253. Salem RM
    254. Sandholm N
    255. Scott RA
    256. Stolk L
    257. Takahashi A
    258. Tanaka T
    259. van 't Hooft FM
    260. Vinkhuyzen AAE
    261. Westra HJ
    262. Zheng W
    263. Zondervan KT
    264. Heath AC
    265. Arveiler D
    266. Bakker SJL
    267. Beilby J
    268. Bergman RN
    269. Blangero J
    270. Bovet P
    271. Campbell H
    272. Caulfield MJ
    273. Cesana G
    274. Chakravarti A
    275. Chasman DI
    276. Chines PS
    277. Collins FS
    278. Crawford DC
    279. Cupples LA
    280. Cusi D
    281. Danesh J
    282. de Faire U
    283. den Ruijter HM
    284. Dominiczak AF
    285. Erbel R
    286. Erdmann J
    287. Eriksson JG
    288. Farrall M
    289. Felix SB
    290. Ferrannini E
    291. Ferrières J
    292. Ford I
    293. Forouhi NG
    294. Forrester T
    295. Franco OH
    296. Gansevoort RT
    297. Gejman PV
    298. Gieger C
    299. Gottesman O
    300. Gudnason V
    301. Gyllensten U
    302. Hall AS
    303. Harris TB
    304. Hattersley AT
    305. Hicks AA
    306. Hindorff LA
    307. Hingorani AD
    308. Hofman A
    309. Homuth G
    310. Hovingh GK
    311. Humphries SE
    312. Hunt SC
    313. Hyppönen E
    314. Illig T
    315. Jacobs KB
    316. Jarvelin MR
    317. Jöckel KH
    318. Johansen B
    319. Jousilahti P
    320. Jukema JW
    321. Jula AM
    322. Kaprio J
    323. Kastelein JJP
    324. Keinanen-Kiukaanniemi SM
    325. Kiemeney LA
    326. Knekt P
    327. Kooner JS
    328. Kooperberg C
    329. Kovacs P
    330. Kraja AT
    331. Kumari M
    332. Kuusisto J
    333. Lakka TA
    334. Langenberg C
    335. Marchand LL
    336. Lehtimäki T
    337. Lyssenko V
    338. Männistö S
    339. Marette A
    340. Matise TC
    341. McKenzie CA
    342. McKnight B
    343. Moll FL
    344. Morris AD
    345. Morris AP
    346. Murray JC
    347. Nelis M
    348. Ohlsson C
    349. Oldehinkel AJ
    350. Ong KK
    351. Madden PAF
    352. Pasterkamp G
    353. Peden JF
    354. Peters A
    355. Postma DS
    356. Pramstaller PP
    357. Price JF
    358. Qi L
    359. Raitakari OT
    360. Rankinen T
    361. Rao DC
    362. Rice TK
    363. Ridker PM
    364. Rioux JD
    365. Ritchie MD
    366. Rudan I
    367. Salomaa V
    368. Samani NJ
    369. Saramies J
    370. Sarzynski MA
    371. Schunkert H
    372. Schwarz PEH
    373. Sever P
    374. Shuldiner AR
    375. Sinisalo J
    376. Stolk RP
    377. Strauch K
    378. Tönjes A
    379. Trégouët DA
    380. Tremblay A
    381. Tremoli E
    382. Virtamo J
    383. Vohl MC
    384. Völker U
    385. Waeber G
    386. Willemsen G
    387. Witteman JC
    388. Zillikens MC
    389. Adair LS
    390. Amouyel P
    391. Asselbergs FW
    392. Assimes TL
    393. Bochud M
    394. Boehm BO
    395. Boerwinkle E
    396. Bornstein SR
    397. Bottinger EP
    398. Bouchard C
    399. Cauchi S
    400. Chambers JC
    401. Chanock SJ
    402. Cooper RS
    403. de Bakker PIW
    404. Dedoussis G
    405. Ferrucci L
    406. Franks PW
    407. Froguel P
    408. Groop LC
    409. Haiman CA
    410. Hamsten A
    411. Hui J
    412. Hunter DJ
    413. Hveem K
    414. Kaplan RC
    415. Kivimaki M
    416. Kuh D
    417. Laakso M
    418. Liu Y
    419. Martin NG
    420. März W
    421. Melbye M
    422. Metspalu A
    423. Moebus S
    424. Munroe PB
    425. Njølstad I
    426. Oostra BA
    427. Palmer CNA
    428. Pedersen NL
    429. Perola M
    430. Pérusse L
    431. Peters U
    432. Power C
    433. Quertermous T
    434. Rauramaa R
    435. Rivadeneira F
    436. Saaristo TE
    437. Saleheen D
    438. Sattar N
    439. Schadt EE
    440. Schlessinger D
    441. Slagboom PE
    442. Snieder H
    443. Spector TD
    444. Thorsteinsdottir U
    445. Stumvoll M
    446. Tuomilehto J
    447. Uitterlinden AG
    448. Uusitupa M
    449. van der Harst P
    450. Walker M
    451. Wallaschofski H
    452. Wareham NJ
    453. Watkins H
    454. Weir DR
    455. Wichmann HE
    456. Wilson JF
    457. Zanen P
    458. Borecki IB
    459. Deloukas P
    460. Fox CS
    461. Heid IM
    462. O'Connell JR
    463. Strachan DP
    464. Stefansson K
    465. van Duijn CM
    466. Abecasis GR
    467. Franke L
    468. Frayling TM
    469. McCarthy MI
    470. Visscher PM
    471. Scherag A
    472. Willer CJ
    473. Boehnke M
    474. Mohlke KL
    475. Lindgren CM
    476. Beckmann JS
    477. Barroso I
    478. North KE
    479. Ingelsson E
    480. Hirschhorn JN
    481. Loos RJF
    482. Speliotes EK
    483. LifeLines Cohort Study
    484. ADIPOGen Consortium
    485. AGEN-BMI Working Group
    486. CARDIOGRAMplusC4D Consortium
    487. CKDGen Consortium
    488. GLGC
    489. ICBP
    490. MAGIC Investigators
    491. MuTHER Consortium
    492. MIGen Consortium
    493. PAGE Consortium
    494. ReproGen Consortium
    495. GENIE Consortium
    496. International Endogene Consortium
    (2015) Genetic studies of body mass index yield new insights for obesity biology
    Nature 518:197–206.
    https://doi.org/10.1038/nature14177
    1. Morris AP
    2. Voight BF
    3. Teslovich TM
    4. Ferreira T
    5. Segrè AV
    6. Steinthorsdottir V
    7. Strawbridge RJ
    8. Khan H
    9. Grallert H
    10. Mahajan A
    11. Prokopenko I
    12. Kang HM
    13. Dina C
    14. Esko T
    15. Fraser RM
    16. Kanoni S
    17. Kumar A
    18. Lagou V
    19. Langenberg C
    20. Luan J
    21. Lindgren CM
    22. Müller-Nurasyid M
    23. Pechlivanis S
    24. Rayner NW
    25. Scott LJ
    26. Wiltshire S
    27. Yengo L
    28. Kinnunen L
    29. Rossin EJ
    30. Raychaudhuri S
    31. Johnson AD
    32. Dimas AS
    33. Loos RJ
    34. Vedantam S
    35. Chen H
    36. Florez JC
    37. Fox C
    38. Liu CT
    39. Rybin D
    40. Couper DJ
    41. Kao WH
    42. Li M
    43. Cornelis MC
    44. Kraft P
    45. Sun Q
    46. van Dam RM
    47. Stringham HM
    48. Chines PS
    49. Fischer K
    50. Fontanillas P
    51. Holmen OL
    52. Hunt SE
    53. Jackson AU
    54. Kong A
    55. Lawrence R
    56. Meyer J
    57. Perry JR
    58. Platou CG
    59. Potter S
    60. Rehnberg E
    61. Robertson N
    62. Sivapalaratnam S
    63. Stančáková A
    64. Stirrups K
    65. Thorleifsson G
    66. Tikkanen E
    67. Wood AR
    68. Almgren P
    69. Atalay M
    70. Benediktsson R
    71. Bonnycastle LL
    72. Burtt N
    73. Carey J
    74. Charpentier G
    75. Crenshaw AT
    76. Doney AS
    77. Dorkhan M
    78. Edkins S
    79. Emilsson V
    80. Eury E
    81. Forsen T
    82. Gertow K
    83. Gigante B
    84. Grant GB
    85. Groves CJ
    86. Guiducci C
    87. Herder C
    88. Hreidarsson AB
    89. Hui J
    90. James A
    91. Jonsson A
    92. Rathmann W
    93. Klopp N
    94. Kravic J
    95. Krjutškov K
    96. Langford C
    97. Leander K
    98. Lindholm E
    99. Lobbens S
    100. Männistö S
    101. Mirza G
    102. Mühleisen TW
    103. Musk B
    104. Parkin M
    105. Rallidis L
    106. Saramies J
    107. Sennblad B
    108. Shah S
    109. Sigurðsson G
    110. Silveira A
    111. Steinbach G
    112. Thorand B
    113. Trakalo J
    114. Veglia F
    115. Wennauer R
    116. Winckler W
    117. Zabaneh D
    118. Campbell H
    119. van Duijn C
    120. Uitterlinden AG
    121. Hofman A
    122. Sijbrands E
    123. Abecasis GR
    124. Owen KR
    125. Zeggini E
    126. Trip MD
    127. Forouhi NG
    128. Syvänen AC
    129. Eriksson JG
    130. Peltonen L
    131. Nöthen MM
    132. Balkau B
    133. Palmer CN
    134. Lyssenko V
    135. Tuomi T
    136. Isomaa B
    137. Hunter DJ
    138. Qi L
    139. Shuldiner AR
    140. Roden M
    141. Barroso I
    142. Wilsgaard T
    143. Beilby J
    144. Hovingh K
    145. Price JF
    146. Wilson JF
    147. Rauramaa R
    148. Lakka TA
    149. Lind L
    150. Dedoussis G
    151. Njølstad I
    152. Pedersen NL
    153. Khaw KT
    154. Wareham NJ
    155. Keinanen-Kiukaanniemi SM
    156. Saaristo TE
    157. Korpi-Hyövälti E
    158. Saltevo J
    159. Laakso M
    160. Kuusisto J
    161. Metspalu A
    162. Collins FS
    163. Mohlke KL
    164. Bergman RN
    165. Tuomilehto J
    166. Boehm BO
    167. Gieger C
    168. Hveem K
    169. Cauchi S
    170. Froguel P
    171. Baldassarre D
    172. Tremoli E
    173. Humphries SE
    174. Saleheen D
    175. Danesh J
    176. Ingelsson E
    177. Ripatti S
    178. Salomaa V
    179. Erbel R
    180. Jöckel KH
    181. Moebus S
    182. Peters A
    183. Illig T
    184. de Faire U
    185. Hamsten A
    186. Morris AD
    187. Donnelly PJ
    188. Frayling TM
    189. Hattersley AT
    190. Boerwinkle E
    191. Melander O
    192. Kathiresan S
    193. Nilsson PM
    194. Deloukas P
    195. Thorsteinsdottir U
    196. Groop LC
    197. Stefansson K
    198. Hu F
    199. Pankow JS
    200. Dupuis J
    201. Meigs JB
    202. Altshuler D
    203. Boehnke M
    204. McCarthy MI
    205. Wellcome Trust Case Control Consortium
    206. Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC) Investigators
    207. Genetic Investigation of ANthropometric Traits (GIANT) Consortium
    208. Asian Genetic Epidemiology Network–Type 2 Diabetes (AGEN-T2D) Consortium
    209. South Asian Type 2 Diabetes (SAT2D) Consortium
    210. DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium
    (2012) Large-scale association analysis provides insights into the genetic architecture and pathophysiology of type 2 diabetes
    Nature Genetics 44:981–990.
    https://doi.org/10.1038/ng.2383
    1. Schunkert H
    2. König IR
    3. Kathiresan S
    4. Reilly MP
    5. Assimes TL
    6. Holm H
    7. Preuss M
    8. Stewart AF
    9. Barbalic M
    10. Gieger C
    11. Absher D
    12. Aherrahrou Z
    13. Allayee H
    14. Altshuler D
    15. Anand SS
    16. Andersen K
    17. Anderson JL
    18. Ardissino D
    19. Ball SG
    20. Balmforth AJ
    21. Barnes TA
    22. Becker DM
    23. Becker LC
    24. Berger K
    25. Bis JC
    26. Boekholdt SM
    27. Boerwinkle E
    28. Braund PS
    29. Brown MJ
    30. Burnett MS
    31. Buysschaert I
    32. Carlquist JF
    33. Chen L
    34. Cichon S
    35. Codd V
    36. Davies RW
    37. Dedoussis G
    38. Dehghan A
    39. Demissie S
    40. Devaney JM
    41. Diemert P
    42. Do R
    43. Doering A
    44. Eifert S
    45. Mokhtari NE
    46. Ellis SG
    47. Elosua R
    48. Engert JC
    49. Epstein SE
    50. de Faire U
    51. Fischer M
    52. Folsom AR
    53. Freyer J
    54. Gigante B
    55. Girelli D
    56. Gretarsdottir S
    57. Gudnason V
    58. Gulcher JR
    59. Halperin E
    60. Hammond N
    61. Hazen SL
    62. Hofman A
    63. Horne BD
    64. Illig T
    65. Iribarren C
    66. Jones GT
    67. Jukema JW
    68. Kaiser MA
    69. Kaplan LM
    70. Kastelein JJ
    71. Khaw KT
    72. Knowles JW
    73. Kolovou G
    74. Kong A
    75. Laaksonen R
    76. Lambrechts D
    77. Leander K
    78. Lettre G
    79. Li M
    80. Lieb W
    81. Loley C
    82. Lotery AJ
    83. Mannucci PM
    84. Maouche S
    85. Martinelli N
    86. McKeown PP
    87. Meisinger C
    88. Meitinger T
    89. Melander O
    90. Merlini PA
    91. Mooser V
    92. Morgan T
    93. Mühleisen TW
    94. Muhlestein JB
    95. Münzel T
    96. Musunuru K
    97. Nahrstaedt J
    98. Nelson CP
    99. Nöthen MM
    100. Olivieri O
    101. Patel RS
    102. Patterson CC
    103. Peters A
    104. Peyvandi F
    105. Qu L
    106. Quyyumi AA
    107. Rader DJ
    108. Rallidis LS
    109. Rice C
    110. Rosendaal FR
    111. Rubin D
    112. Salomaa V
    113. Sampietro ML
    114. Sandhu MS
    115. Schadt E
    116. Schäfer A
    117. Schillert A
    118. Schreiber S
    119. Schrezenmeir J
    120. Schwartz SM
    121. Siscovick DS
    122. Sivananthan M
    123. Sivapalaratnam S
    124. Smith A
    125. Smith TB
    126. Snoep JD
    127. Soranzo N
    128. Spertus JA
    129. Stark K
    130. Stirrups K
    131. Stoll M
    132. Tang WH
    133. Tennstedt S
    134. Thorgeirsson G
    135. Thorleifsson G
    136. Tomaszewski M
    137. Uitterlinden AG
    138. van Rij AM
    139. Voight BF
    140. Wareham NJ
    141. Wells GA
    142. Wichmann HE
    143. Wild PS
    144. Willenborg C
    145. Witteman JC
    146. Wright BJ
    147. Ye S
    148. Zeller T
    149. Ziegler A
    150. Cambien F
    151. Goodall AH
    152. Cupples LA
    153. Quertermous T
    154. März W
    155. Hengstenberg C
    156. Blankenberg S
    157. Ouwehand WH
    158. Hall AS
    159. Deloukas P
    160. Thompson JR
    161. Stefansson K
    162. Roberts R
    163. Thorsteinsdottir U
    164. O'Donnell CJ
    165. McPherson R
    166. Erdmann J
    167. Samani NJ
    168. Cardiogenics
    169. CARDIoGRAM Consortium
    (2011) Large-scale association analysis identifies 13 new susceptibility loci for coronary artery disease
    Nature Genetics 43:333–338.
    https://doi.org/10.1038/ng.784
    1. Selleri L
    2. Depew MJ
    3. Jacobs Y
    4. Chanda SK
    5. Tsang KY
    6. Cheah KS
    7. Rubenstein JL
    8. O'Gorman S
    9. Cleary ML
    (2001)
    Requirement for Pbx1 in skeletal patterning and programming chondrocyte proliferation and differentiation
    Development 128:3543–3557.
    1. Teslovich TM
    2. Musunuru K
    3. Smith AV
    4. Edmondson AC
    5. Stylianou IM
    6. Koseki M
    7. Pirruccello JP
    8. Ripatti S
    9. Chasman DI
    10. Willer CJ
    11. Johansen CT
    12. Fouchier SW
    13. Isaacs A
    14. Peloso GM
    15. Barbalic M
    16. Ricketts SL
    17. Bis JC
    18. Aulchenko YS
    19. Thorleifsson G
    20. Feitosa MF
    21. Chambers J
    22. Orho-Melander M
    23. Melander O
    24. Johnson T
    25. Li X
    26. Guo X
    27. Li M
    28. Shin Cho Y
    29. Jin Go M
    30. Jin Kim Y
    31. Lee JY
    32. Park T
    33. Kim K
    34. Sim X
    35. Twee-Hee Ong R
    36. Croteau-Chonka DC
    37. Lange LA
    38. Smith JD
    39. Song K
    40. Hua Zhao J
    41. Yuan X
    42. Luan J
    43. Lamina C
    44. Ziegler A
    45. Zhang W
    46. Zee RY
    47. Wright AF
    48. Witteman JC
    49. Wilson JF
    50. Willemsen G
    51. Wichmann HE
    52. Whitfield JB
    53. Waterworth DM
    54. Wareham NJ
    55. Waeber G
    56. Vollenweider P
    57. Voight BF
    58. Vitart V
    59. Uitterlinden AG
    60. Uda M
    61. Tuomilehto J
    62. Thompson JR
    63. Tanaka T
    64. Surakka I
    65. Stringham HM
    66. Spector TD
    67. Soranzo N
    68. Smit JH
    69. Sinisalo J
    70. Silander K
    71. Sijbrands EJ
    72. Scuteri A
    73. Scott J
    74. Schlessinger D
    75. Sanna S
    76. Salomaa V
    77. Saharinen J
    78. Sabatti C
    79. Ruokonen A
    80. Rudan I
    81. Rose LM
    82. Roberts R
    83. Rieder M
    84. Psaty BM
    85. Pramstaller PP
    86. Pichler I
    87. Perola M
    88. Penninx BW
    89. Pedersen NL
    90. Pattaro C
    91. Parker AN
    92. Pare G
    93. Oostra BA
    94. O'Donnell CJ
    95. Nieminen MS
    96. Nickerson DA
    97. Montgomery GW
    98. Meitinger T
    99. McPherson R
    100. McCarthy MI
    101. McArdle W
    102. Masson D
    103. Martin NG
    104. Marroni F
    105. Mangino M
    106. Magnusson PK
    107. Lucas G
    108. Luben R
    109. Loos RJ
    110. Lokki ML
    111. Lettre G
    112. Langenberg C
    113. Launer LJ
    114. Lakatta EG
    115. Laaksonen R
    116. Kyvik KO
    117. Kronenberg F
    118. König IR
    119. Khaw KT
    120. Kaprio J
    121. Kaplan LM
    122. Johansson A
    123. Jarvelin MR
    124. Janssens AC
    125. Ingelsson E
    126. Igl W
    127. Kees Hovingh G
    128. Hottenga JJ
    129. Hofman A
    130. Hicks AA
    131. Hengstenberg C
    132. Heid IM
    133. Hayward C
    134. Havulinna AS
    135. Hastie ND
    136. Harris TB
    137. Haritunians T
    138. Hall AS
    139. Gyllensten U
    140. Guiducci C
    141. Groop LC
    142. Gonzalez E
    143. Gieger C
    144. Freimer NB
    145. Ferrucci L
    146. Erdmann J
    147. Elliott P
    148. Ejebe KG
    149. Döring A
    150. Dominiczak AF
    151. Demissie S
    152. Deloukas P
    153. de Geus EJ
    154. de Faire U
    155. Crawford G
    156. Collins FS
    157. Chen YD
    158. Caulfield MJ
    159. Campbell H
    160. Burtt NP
    161. Bonnycastle LL
    162. Boomsma DI
    163. Boekholdt SM
    164. Bergman RN
    165. Barroso I
    166. Bandinelli S
    167. Ballantyne CM
    168. Assimes TL
    169. Quertermous T
    170. Altshuler D
    171. Seielstad M
    172. Wong TY
    173. Tai ES
    174. Feranil AB
    175. Kuzawa CW
    176. Adair LS
    177. Taylor HA
    178. Borecki IB
    179. Gabriel SB
    180. Wilson JG
    181. Holm H
    182. Thorsteinsdottir U
    183. Gudnason V
    184. Krauss RM
    185. Mohlke KL
    186. Ordovas JM
    187. Munroe PB
    188. Kooner JS
    189. Tall AR
    190. Hegele RA
    191. Kastelein JJ
    192. Schadt EE
    193. Rotter JI
    194. Boerwinkle E
    195. Strachan DP
    196. Mooser V
    197. Stefansson K
    198. Reilly MP
    199. Samani NJ
    200. Schunkert H
    201. Cupples LA
    202. Sandhu MS
    203. Ridker PM
    204. Rader DJ
    205. van Duijn CM
    206. Peltonen L
    207. Abecasis GR
    208. Boehnke M
    209. Kathiresan S
    (2010) Biological, clinical and population relevance of 95 loci for blood lipids
    Nature 466:707–713.
    https://doi.org/10.1038/nature09270
    1. van der Harst P
    2. Zhang W
    3. Mateo Leach I
    4. Rendon A
    5. Verweij N
    6. Sehmi J
    7. Paul DS
    8. Elling U
    9. Allayee H
    10. Li X
    11. Radhakrishnan A
    12. Tan ST
    13. Voss K
    14. Weichenberger CX
    15. Albers CA
    16. Al-Hussani A
    17. Asselbergs FW
    18. Ciullo M
    19. Danjou F
    20. Dina C
    21. Esko T
    22. Evans DM
    23. Franke L
    24. Gögele M
    25. Hartiala J
    26. Hersch M
    27. Holm H
    28. Hottenga JJ
    29. Kanoni S
    30. Kleber ME
    31. Lagou V
    32. Langenberg C
    33. Lopez LM
    34. Lyytikäinen LP
    35. Melander O
    36. Murgia F
    37. Nolte IM
    38. O'Reilly PF
    39. Padmanabhan S
    40. Parsa A
    41. Pirastu N
    42. Porcu E
    43. Portas L
    44. Prokopenko I
    45. Ried JS
    46. Shin SY
    47. Tang CS
    48. Teumer A
    49. Traglia M
    50. Ulivi S
    51. Westra HJ
    52. Yang J
    53. Zhao JH
    54. Anni F
    55. Abdellaoui A
    56. Attwood A
    57. Balkau B
    58. Bandinelli S
    59. Bastardot F
    60. Benyamin B
    61. Boehm BO
    62. Cookson WO
    63. Das D
    64. de Bakker PI
    65. de Boer RA
    66. de Geus EJ
    67. de Moor MH
    68. Dimitriou M
    69. Domingues FS
    70. Döring A
    71. Engström G
    72. Eyjolfsson GI
    73. Ferrucci L
    74. Fischer K
    75. Galanello R
    76. Garner SF
    77. Genser B
    78. Gibson QD
    79. Girotto G
    80. Gudbjartsson DF
    81. Harris SE
    82. Hartikainen AL
    83. Hastie CE
    84. Hedblad B
    85. Illig T
    86. Jolley J
    87. Kähönen M
    88. Kema IP
    89. Kemp JP
    90. Liang L
    91. Lloyd-Jones H
    92. Loos RJ
    93. Meacham S
    94. Medland SE
    95. Meisinger C
    96. Memari Y
    97. Mihailov E
    98. Miller K
    99. Moffatt MF
    100. Nauck M
    101. Novatchkova M
    102. Nutile T
    103. Olafsson I
    104. Onundarson PT
    105. Parracciani D
    106. Penninx BW
    107. Perseu L
    108. Piga A
    109. Pistis G
    110. Pouta A
    111. Puc U
    112. Raitakari O
    113. Ring SM
    114. Robino A
    115. Ruggiero D
    116. Ruokonen A
    117. Saint-Pierre A
    118. Sala C
    119. Salumets A
    120. Sambrook J
    121. Schepers H
    122. Schmidt CO
    123. Silljé HH
    124. Sladek R
    125. Smit JH
    126. Starr JM
    127. Stephens J
    128. Sulem P
    129. Tanaka T
    130. Thorsteinsdottir U
    131. Tragante V
    132. van Gilst WH
    133. van Pelt LJ
    134. van Veldhuisen DJ
    135. Völker U
    136. Whitfield JB
    137. Willemsen G
    138. Winkelmann BR
    139. Wirnsberger G
    140. Algra A
    141. Cucca F
    142. d'Adamo AP
    143. Danesh J
    144. Deary IJ
    145. Dominiczak AF
    146. Elliott P
    147. Fortina P
    148. Froguel P
    149. Gasparini P
    150. Greinacher A
    151. Hazen SL
    152. Jarvelin MR
    153. Khaw KT
    154. Lehtimäki T
    155. Maerz W
    156. Martin NG
    157. Metspalu A
    158. Mitchell BD
    159. Montgomery GW
    160. Moore C
    161. Navis G
    162. Pirastu M
    163. Pramstaller PP
    164. Ramirez-Solis R
    165. Schadt E
    166. Scott J
    167. Shuldiner AR
    168. Smith GD
    169. Smith JG
    170. Snieder H
    171. Sorice R
    172. Spector TD
    173. Stefansson K
    174. Stumvoll M
    175. Tang WH
    176. Toniolo D
    177. Tönjes A
    178. Visscher PM
    179. Vollenweider P
    180. Wareham NJ
    181. Wolffenbuttel BH
    182. Boomsma DI
    183. Beckmann JS
    184. Dedoussis GV
    185. Deloukas P
    186. Ferreira MA
    187. Sanna S
    188. Uda M
    189. Hicks AA
    190. Penninger JM
    191. Gieger C
    192. Kooner JS
    193. Ouwehand WH
    194. Soranzo N
    195. Chambers JC
    (2012) Seventy-five genetic loci influencing the human red blood cell
    Nature 492:369–375.
    https://doi.org/10.1038/nature11677
    1. Wood AR
    2. Esko T
    3. Yang J
    4. Vedantam S
    5. Pers TH
    6. Gustafsson S
    7. Chu AY
    8. Estrada K
    9. Luan J
    10. Kutalik Z
    11. Amin N
    12. Buchkovich ML
    13. Croteau-Chonka DC
    14. Day FR
    15. Duan Y
    16. Fall T
    17. Fehrmann R
    18. Ferreira T
    19. Jackson AU
    20. Karjalainen J
    21. Lo KS
    22. Locke AE
    23. Mägi R
    24. Mihailov E
    25. Porcu E
    26. Randall JC
    27. Scherag A
    28. Vinkhuyzen AA
    29. Westra HJ
    30. Winkler TW
    31. Workalemahu T
    32. Zhao JH
    33. Absher D
    34. Albrecht E
    35. Anderson D
    36. Baron J
    37. Beekman M
    38. Demirkan A
    39. Ehret GB
    40. Feenstra B
    41. Feitosa MF
    42. Fischer K
    43. Fraser RM
    44. Goel A
    45. Gong J
    46. Justice AE
    47. Kanoni S
    48. Kleber ME
    49. Kristiansson K
    50. Lim U
    51. Lotay V
    52. Lui JC
    53. Mangino M
    54. Mateo Leach I
    55. Medina-Gomez C
    56. Nalls MA
    57. Nyholt DR
    58. Palmer CD
    59. Pasko D
    60. Pechlivanis S
    61. Prokopenko I
    62. Ried JS
    63. Ripke S
    64. Shungin D
    65. Stancáková A
    66. Strawbridge RJ
    67. Sung YJ
    68. Tanaka T
    69. Teumer A
    70. Trompet S
    71. van der Laan SW
    72. van Setten J
    73. Van Vliet-Ostaptchouk JV
    74. Wang Z
    75. Yengo L
    76. Zhang W
    77. Afzal U
    78. Arnlöv J
    79. Arscott GM
    80. Bandinelli S
    81. Barrett A
    82. Bellis C
    83. Bennett AJ
    84. Berne C
    85. Blüher M
    86. Bolton JL
    87. Böttcher Y
    88. Boyd HA
    89. Bruinenberg M
    90. Buckley BM
    91. Buyske S
    92. Caspersen IH
    93. Chines PS
    94. Clarke R
    95. Claudi-Boehm S
    96. Cooper M
    97. Daw EW
    98. De Jong PA
    99. Deelen J
    100. Delgado G
    101. Denny JC
    102. Dhonukshe-Rutten R
    103. Dimitriou M
    104. Doney AS
    105. Dörr M
    106. Eklund N
    107. Eury E
    108. Folkersen L
    109. Garcia ME
    110. Geller F
    111. Giedraitis V
    112. Go AS
    113. Grallert H
    114. Grammer TB
    115. Gräßler J
    116. Grönberg H
    117. de Groot LC
    118. Groves CJ
    119. Haessler J
    120. Hall P
    121. Haller T
    122. Hallmans G
    123. Hannemann A
    124. Hartman CA
    125. Hassinen M
    126. Hayward C
    127. Heard-Costa NL
    128. Helmer Q
    129. Hemani G
    130. Henders AK
    131. Hillege HL
    132. Hlatky MA
    133. Hoffmann W
    134. Hoffmann P
    135. Holmen O
    136. Houwing-Duistermaat JJ
    137. Illig T
    138. Isaacs A
    139. James AL
    140. Jeff J
    141. Johansen B
    142. Johansson Å
    143. Jolley J
    144. Juliusdottir T
    145. Junttila J
    146. Kho AN
    147. Kinnunen L
    148. Klopp N
    149. Kocher T
    150. Kratzer W
    151. Lichtner P
    152. Lind L
    153. Lindström J
    154. Lobbens S
    155. Lorentzon M
    156. Lu Y
    157. Lyssenko V
    158. Magnusson PK
    159. Mahajan A
    160. Maillard M
    161. McArdle WL
    162. McKenzie CA
    163. McLachlan S
    164. McLaren PJ
    165. Menni C
    166. Merger S
    167. Milani L
    168. Moayyeri A
    169. Monda KL
    170. Morken MA
    171. Müller G
    172. Müller-Nurasyid M
    173. Musk AW
    174. Narisu N
    175. Nauck M
    176. Nolte IM
    177. Nöthen MM
    178. Oozageer L
    179. Pilz S
    180. Rayner NW
    181. Renstrom F
    182. Robertson NR
    183. Rose LM
    184. Roussel R
    185. Sanna S
    186. Scharnagl H
    187. Scholtens S
    188. Schumacher FR
    189. Schunkert H
    190. Scott RA
    191. Sehmi J
    192. Seufferlein T
    193. Shi J
    194. Silventoinen K
    195. Smit JH
    196. Smith AV
    197. Smolonska J
    198. Stanton AV
    199. Stirrups K
    200. Stott DJ
    201. Stringham HM
    202. Sundström J
    203. Swertz MA
    204. Syvänen AC
    205. Tayo BO
    206. Thorleifsson G
    207. Tyrer JP
    208. van Dijk S
    209. van Schoor NM
    210. van der Velde N
    211. van Heemst D
    212. van Oort FV
    213. Vermeulen SH
    214. Verweij N
    215. Vonk JM
    216. Waite LL
    217. Waldenberger M
    218. Wennauer R
    219. Wilkens LR
    220. Willenborg C
    221. Wilsgaard T
    222. Wojczynski MK
    223. Wong A
    224. Wright AF
    225. Zhang Q
    226. Arveiler D
    227. Bakker SJ
    228. Beilby J
    229. Bergman RN
    230. Bergmann S
    231. Biffar R
    232. Blangero J
    233. Boomsma DI
    234. Bornstein SR
    235. Bovet P
    236. Brambilla P
    237. Brown MJ
    238. Campbell H
    239. Caulfield MJ
    240. Chakravarti A
    241. Collins R
    242. Collins FS
    243. Crawford DC
    244. Cupples LA
    245. Danesh J
    246. de Faire U
    247. den Ruijter HM
    248. Erbel R
    249. Erdmann J
    250. Eriksson JG
    251. Farrall M
    252. Ferrannini E
    253. Ferrières J
    254. Ford I
    255. Forouhi NG
    256. Forrester T
    257. Gansevoort RT
    258. Gejman PV
    259. Gieger C
    260. Golay A
    261. Gottesman O
    262. Gudnason V
    263. Gyllensten U
    264. Haas DW
    265. Hall AS
    266. Harris TB
    267. Hattersley AT
    268. Heath AC
    269. Hengstenberg C
    270. Hicks AA
    271. Hindorff LA
    272. Hingorani AD
    273. Hofman A
    274. Hovingh GK
    275. Humphries SE
    276. Hunt SC
    277. Hypponen E
    278. Jacobs KB
    279. Jarvelin MR
    280. Jousilahti P
    281. Jula AM
    282. Kaprio J
    283. Kastelein JJ
    284. Kayser M
    285. Kee F
    286. Keinanen-Kiukaanniemi SM
    287. Kiemeney LA
    288. Kooner JS
    289. Kooperberg C
    290. Koskinen S
    291. Kovacs P
    292. Kraja AT
    293. Kumari M
    294. Kuusisto J
    295. Lakka TA
    296. Langenberg C
    297. Le Marchand L
    298. Lehtimäki T
    299. Lupoli S
    300. Madden PA
    301. Männistö S
    302. Manunta P
    303. Marette A
    304. Matise TC
    305. McKnight B
    306. Meitinger T
    307. Moll FL
    308. Montgomery GW
    309. Morris AD
    310. Morris AP
    311. Murray JC
    312. Nelis M
    313. Ohlsson C
    314. Oldehinkel AJ
    315. Ong KK
    316. Ouwehand WH
    317. Pasterkamp G
    318. Peters A
    319. Pramstaller PP
    320. Price JF
    321. Qi L
    322. Raitakari OT
    323. Rankinen T
    324. Rao DC
    325. Rice TK
    326. Ritchie M
    327. Rudan I
    328. Salomaa V
    329. Samani NJ
    330. Saramies J
    331. Sarzynski MA
    332. Schwarz PE
    333. Sebert S
    334. Sever P
    335. Shuldiner AR
    336. Sinisalo J
    337. Steinthorsdottir V
    338. Stolk RP
    339. Tardif JC
    340. Tönjes A
    341. Tremblay A
    342. Tremoli E
    343. Virtamo J
    344. Vohl MC
    345. Amouyel P
    346. Asselbergs FW
    347. Assimes TL
    348. Bochud M
    349. Boehm BO
    350. Boerwinkle E
    351. Bottinger EP
    352. Bouchard C
    353. Cauchi S
    354. Chambers JC
    355. Chanock SJ
    356. Cooper RS
    357. de Bakker PI
    358. Dedoussis G
    359. Ferrucci L
    360. Franks PW
    361. Froguel P
    362. Groop LC
    363. Haiman CA
    364. Hamsten A
    365. Hayes MG
    366. Hui J
    367. Hunter DJ
    368. Hveem K
    369. Jukema JW
    370. Kaplan RC
    371. Kivimaki M
    372. Kuh D
    373. Laakso M
    374. Liu Y
    375. Martin NG
    376. März W
    377. Melbye M
    378. Moebus S
    379. Munroe PB
    380. Njølstad I
    381. Oostra BA
    382. Palmer CN
    383. Pedersen NL
    384. Perola M
    385. Pérusse L
    386. Peters U
    387. Powell JE
    388. Power C
    389. Quertermous T
    390. Rauramaa R
    391. Reinmaa E
    392. Ridker PM
    393. Rivadeneira F
    394. Rotter JI
    395. Saaristo TE
    396. Saleheen D
    397. Schlessinger D
    398. Slagboom PE
    399. Snieder H
    400. Spector TD
    401. Strauch K
    402. Stumvoll M
    403. Tuomilehto J
    404. Uusitupa M
    405. van der Harst P
    406. Völzke H
    407. Walker M
    408. Wareham NJ
    409. Watkins H
    410. Wichmann HE
    411. Wilson JF
    412. Zanen P
    413. Deloukas P
    414. Heid IM
    415. Lindgren CM
    416. Mohlke KL
    417. Speliotes EK
    418. Thorsteinsdottir U
    419. Barroso I
    420. Fox CS
    421. North KE
    422. Strachan DP
    423. Beckmann JS
    424. Berndt SI
    425. Boehnke M
    426. Borecki IB
    427. McCarthy MI
    428. Metspalu A
    429. Stefansson K
    430. Uitterlinden AG
    431. van Duijn CM
    432. Franke L
    433. Willer CJ
    434. Price AL
    435. Lettre G
    436. Loos RJ
    437. Weedon MN
    438. Ingelsson E
    439. O'Connell JR
    440. Abecasis GR
    441. Chasman DI
    442. Goddard ME
    443. Visscher PM
    444. Hirschhorn JN
    445. Frayling TM
    446. Electronic Medical Records and Genomics (eMEMERGEGE) Consortium
    447. MIGen Consortium
    448. PAGEGE Consortium
    449. LifeLines Cohort Study
    (2014) Defining the role of common variation in the genomic and biological architecture of adult human height
    Nature Genetics 46:1173–1186.
    https://doi.org/10.1038/ng.3097

Decision letter

  1. Clifford J Rosen
    Reviewing Editor; Maine Medical Center Research Institute, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Epigenetic profiling of growth plate chondrocytes sheds insight into regulatory genetic variation influencing height" for consideration by eLife. Your article has been favorably evaluated by Mark McCarthy (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

All the reviewers were enthusiastic about your comprehensive generation of epigenetic data for the growth plate and the potential utilization for future studies on genetic regulatory regions. The ATAC-seq approach and the overlay of GWAS data will provide significant insights into murine and possibly human growth determinants.

There were three major concerns that need to be addressed before a final decision can be made.

First delineation of sex from the E15.5 murine chondrocytes and justification for the E15.5 time point are essential for understanding the framework for growth variants and for future studies by other investigators.

Second, the manuscript still lacks a strong connection to chondrocyte biology as none of the links established by the data directly address (or explain) chondrocyte proliferation or hypertrophy and genes well known to regulate them positively or negatively (PTHrp, Indian hedgehog, FGFs, IGFs, etc.).

Third, although different genes have been identified as possibly being putative regulatory variants there are no experimental or computational evidence implicating the ATAC-seq peaks in regulating the nearby growth plate genes, nor altering transcriptional regulatory activity.

These issues underlie a perception that, independent of your Nature Genetics paper, there is a lack in this manuscript of mechanistic evidence that the epigenomic profiles can be linked to functional correlates. For example, there is a suggestion that rs9920291 disrupted HOXD13 binding, but there is no evidence that HOXD13 actually binds this site at the endogenous locus. As computationally predicted TF motifs are prevalent in the genome, it is hard to assess whether this SNP/motif overlap is biologically meaningful. This requirement for a mechanistic/functional connection may require additional work.

These major points will need to be addressed in the revision along with several other points imbedded in the three reviewers’ comments that have been included below in the comments, and should be considered carefully.

Reviewer #1:

In this paper, Capellini and colleagues utilize a novel method (ATAC-seq) to read the open chromatin regulatory landscape of growth plate chondrocytes in neonatal mice. They were then able to analyze height GWAS data in the context of femur growth to gain new insights into how height GWAS variants act at the growth plate. Briefly they showed that ATAC-seq datasets generated from mouse growth plate chondrocytes overlap with previously identified signals of chondrocyte biology in the mouse, as well as in humans. They identified locations of active regulatory enhancers and Sox9 binding sites as assessed using ChIP seq on mouse rib chondrocytes and provided datasets that overlap with evidence of regulatory control during chondrogenesis and human limb embryogenesis. Their work is a nice supplement to their Nature Genetics paper identifying an ancestral regulatory variant in 'Grow'1 enhancer in mouse limbs that also affects human chondrocyte biology. Finally, the authors detected a variant in the chondroitin sulfate synthase 1 gene which might influence binding of HOXD13, a trans acting regulator of chondrocyte growth.

Strengths: There are several strengths to this paper: First the authors were able to provide very comprehensive data sets on the epigenetic landscape of mouse, and human chondrocytes, an important breakthrough in an otherwise very difficult organ to characterize. Second, they have provided evidence that these regions are critically important for growth of the chondrocyte and thus provide support from previous studies of candidate genes in chondrogenesis such as GDF5, 6, Bmp5 etc. that others have defined. Third they have been able to overlay height GWAS variants that are also significantly enriched in growth plate chondrocyte regulatory regions, particularly for those that are intergenic. Finally, this is state of the art technology in epigenetics; i.e. the use of ATAC methodology, importantly it requires fewer chondrocytes to peruse the chromatin landscape thus breaking new ground.

Limitations: – there are some concerns that need to be fully addressed:

1) The authors used E15.5 embryos from FVB mice to study the epigenetic regions for growth plate regulation. But inbred strains often exhibit numerous nucleotide variants; since C57BL6J mice are the most frequently utilized strain for in vitro and in vivo chondrocyte (and osteoblast) studies, have the authors determined whether the specific nucleotide sequence in the CHYS1 variant, rs9920291, is also present in B6 mouse chondrocytes?

2) Sex is a major determinant of height particularly early in life and during puberty. The early influences of sex steroids likely begins during late embryonic development; i.e. after E11.5 when SRY is expressed in the males, but not females. It would be important to fully understand the epigenetic landscape in respect to sex differences in vivo. Hence, the embryos should have been typed for sex, and both male and female embryos should be examined to more fully understand the transcriptional regulation of chondrocyte differentiation.

3) Can the authors justify why they did not follow up with functional studies on the CHYS1 variant rs9920291 which they identified as being a cis acting site for HOXD13.

4) Figure 6B of the E15.5 embryo showing the lacZ reporter from the ancestral GROW1 enhancer appears nearly identical to their Nature Genetics paper Figure 4B; Can the authors defend its use in this manuscript? Otherwise it should be deleted.

Reviewer #2:

In this study, the authors have performed epigenetic profiling of active chromatin regions in E15.5 mouse embryo femurs to establish possible links with GWAS variants associated with human height previously described by other groups. The regions were found to positively associate with known genes regulating chondrocyte and skeletal biology, particularly those encoding growth plate transcription factors. These and other data lead the authors to conclude that integration of epigenetic information with genetic association data can identify mechanisms important for human height.

The study comes from a research group with recognized expertise in genetics and genomics and for important contributions to the field. It also addresses human height that can be severely affected by congenital conditions, thus introducing an element of translational medicine relevance to the work. The extensive data combining epigenetic profiling using ATAT-seq with several bioinformatics and analytical tools provide a wealth of possible links with specific genes and transcriptome data sets from previous studies. Despite these positive assets, the study raises several concerns.

While the study is technically strong, it is deficient in skeletal cell and developmental biology, raising questions about the major conclusions reached and specific relevance to human height.

The growth plate is certainly the key structure determining long bone and vertebral elongation and thus, body height. It is clearly established that the rates of skeletal growth are determined by chondrocyte proliferation (minor contribution) and chondrocyte hypertrophy (major contribution) in the growth plate. It is also well established that different growth plates within the same organism at any given time "grow" at very different rates, such as proximal and distal growth plates in radius and ulna. None of the links established by the data presented here directly address (or explain) chondrocyte proliferation or hypertrophy and genes known to regulate them positively or negatively (PTHrp, Indian hedgehog, FGFs, IGFs, etc.). The only links established are to general aspects of chondrocyte and skeletal biology, thus making it very difficult to understand how the links have anything to do with "height". The behavior of the growth plate is controlled by intrinsic mechanisms (see above) as well as systemic factors such as growth hormone. Again, none of the links related to this or other systemic cues.

The authors made use of the proximal and distal portions of E15.5 mouse embryo femur and associated perichondrium as starting material for the epigenetic analysis. The tissues were digested into single cell suspension enzymatically before processing for ATAT-seq. There is no explanation of why this developmental stage was chosen and why both ends of the femur were considered. Also, the lengthy enzymatic digestion could have altered gene expression and chromatin configuration.

To establish stronger links, the authors considered the most compelling GWAS variants overlapping ATAT-seq peaks, finding 59 variants distributed over 26 loci. The two examples mentioned are CHSY1 and HOXD13 neither of which has established roles in growth plate function. Indeed, most evidence indicates that HOXD13 regulates autopod skeletal patterning and morphogenesis, but not rates of growth.

The authors make a strong case that the data can provide insights into the meaning and functional relevance of GWAS correlations, but I do not see how they provide any direct insights into what regulates human height. It may be worth considering an overall change in rationale for the study.

Reviewer #3:

The article by Guo et al. describes the epigenetic profiling of mouse growth plate chondrocytes to elucidate potential molecular mechanisms through which genetic variants contribute to human height. The authors perform ATAC-seq profiling to identify open chromatin regions from murine femoral growth plates and demonstrate that femur open chromatin are (i) enriched nearby skeletal development genes, and (ii) enriched for human height GWAS loci, above other epigenomic datasets available from mouse and human. While it has previously been shown that height GWAS loci are enriched nearby musculoskeletal genes (e.g. Figure 3A of Wood et. al 2014), the identification of musculoskeletal regulatory elements possibly affected by height loci is valuable for future studies. The authors proceed to identify height GWAS loci and specific SNPs that overlap femur open chromatin peaks, and provide potential hypotheses for how select SNPs (e.g. rs9920291) may act to influence human height. Unfortunately, despite claiming the identification of "compelling mechanisms for GWAS variants" (Abstract), the authors do not provide sufficient experimental or computational evidence in support, except for rs4911178 that references data in Capellini et al. 2017 (PMID 28671685). In summary, the authors present an important epigenomic dataset that can act as a first step for future study of height GWAS SNPs, but do not provide "strong evidence" (subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, first paragraph) for compelling mechanisms at any height GWAS loci beyond referencing their previously published paper.

1) Substantially fewer GWAS loci overlap ATAC-seq peaks when using the fine-mapped 95% credible variants, instead of all SNPs with r2>0.5 (e.g. 192 loci vs. 317 loci, subsection “Fine-mapped height variants overlap with femoral growth plate open chromatin regions” and subsection “Human height variants are enriched in femoral open chromatin regions”, last paragraph; also 26 loci vs. 46 loci, subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, first paragraph and subsection “A subset of human height variants in femoral open chromatin regions reside near genes differentially expressed in the growth plate”, second paragraph). Given that r2>0.5 is not a strict LD threshold (e.g. Wood et al. 2014 uses r2>0.8 for strict LD), is it still meaningful to report overlap statistics using the r2>0.5 set when the PICS 95% credible set is available?

2) The authors state "44 different genes have been identified as possibly being modulated by these putative regulatory variants". However, in the present manuscript the authors do not provide experimental or computational evidence implicating the ATAC-seq peaks in regulating the nearby growth plate genes, and provide limited evidence that any of these overlapping variants can alter transcriptional regulatory activity.

3) Following point #2, and to support the authors' claim of identifying "prime targets for functional testing", it would be helpful to know how many GWAS SNP overlaps with femur open chromatin regions the authors believe to be biologically meaningful versus coincidence.

4) Similar to point #2, it is premature to state "59 variants thus have strong evidence for being potentially causal variants". For the authors' example of the CHSY1 locus, while there are 4 SNPs in the PICS 95% set that overlap an ATAC-seq peak nearby CHSY1, Supplementary file 6 reveals there are 87 SNPs in the 95% credible set at the CHSY1 locus. Without more evidence, it is difficult to agree with the potential causal variant conclusion.

5) Subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, first paragraph: The authors suggest rs9920291 disrupts HOXD13 binding, but do not provide evidence that HOXD13 actually binds this site at the endogenous locus. As computationally predicted TF motifs are prevalent in the genome, it is hard to assess whether this SNP/motif overlap is biologically meaningful.

6) Wood et al. (2014) found that height GWAS loci were more likely to be eQTLs (in blood) for genes expressed in cartilage. Did the authors examine whether the 468 credible SNPs overlapping femur open chromatin peaks are enriched for eQTL signals? This analysis may strengthen the authors' mechanistic claims.

https://doi.org/10.7554/eLife.29329.036

Author response

[…] There were three major concerns that need to be addressed before a final decision can be made.

First delineation of sex from the E15.5 murine chondrocytes and justification for the E15.5 time point are essential for understanding the framework for growth variants and for future studies by other investigators.

We thank the Reviewing Editor for summarizing these two main concerns and address each below separately.

While we agree that sex is a major determinant of height in humans particularly during puberty, our method of cell collection for ATAC-seq is not yet conducive to examining this aspect of biology. We have optimized our protocol to acquire large numbers of viable chondrocytes that have not been removed from their in vivocontext for extended periods of time.

First, we perform our ATAC-seq on freshly acquired chondrocytes, extracted from the embryonic growth plates, and we must pool embryonic long bone ends across a litter to collect enough for each specific growth plate end (i.e. proximal versus distal femur). In the process, any significant delay in extraction can significantly impact chromatin state, as the cells are no longer fresh and instead are subjected to prolonged exposure to artificial media and in vitro conditions. Genotyping for sex constitutes a major delay in the process.

Second, we perform our ATAC-seq on E15.5 embryos, a stage when sex differences are not readily apparent in the lengths of the skeletal elements (note: embryonic mice in general are not proper for studying sex-specific effects in the skeleton). However, we intentionally chose E15.5 for a number of reasons (see next point below), most notably our ability to easily extract large numbers of viable chondrocytes from the limited extracellular matrix environment using a brief collagenase digestion step (<2 hours). Our current attempts to extract chondrocytes at much later stages (i.e., E18.5 or post-natal stages – when sex differences begin to become apparent in mice) require at least 12 hours of collagenase digestion, but result in an increased ECM content that interferes with ATAC-seq and more importantly a significant increase in cell death (>25%), which results in low quality and biased ATAC-seq libraries. In addition, the prolonged duration of exposure to collagenase and artificial media in vitro that comes hand-in-hand with studying these later stages likely impacts chromatin profiles. In time, when chondrocyte extraction techniques improve, we plan to extract chondrocytes from single early post-natal mice, genotype each embryo for sex while we perform ATAC-seq, and explore sex-specific differences in growth plate chondrocyte chromatin biology.

We chose E15.5 as a stage to explore chondrocyte growth plate biology for a number of reasons:

First, growth plates at this stage show a very limited amount of extracellular matrix, which can interfere with the ATAC-seq protocol. For example, extensive extracellular matrix, characteristic of later fetal and post-natal stages, requires long periods of collagenase digestion (e.g., greater than 12 hours) to isolate chondrocytes, which in turn has a significant impact on the nature of open chromatin profiles in the resultant extracted cell populations. Increased collagenase digestion also increases cell death which also influences chromatin state biology.

Second, growth plates are readily dissected at E15.5 since soft tissues, such as tendon, ligament, and muscle, strip easily off the growth plate using simple mechanical force rather than enzymatic digestions that can alter chromatin state profiles. At later stages, increased enzymatic concentrations and processing times, as well as more mechanical force are necessary to strip soft tissues off of growth plates, and the effects of these treatments on ATAC-seq libraries remain

problematic and unclear.

Third, growth plates at E15.5 display all major chondrocyte differentiated cell types (i.e., spanning resting, proliferative, pre-hypertrophic, hypertrophic, along with perichondrial zones) and at the same time lack evidence in and around the epiphysis of secondary ossification center formation. Thus, we when we extract at this stage, we are less concerned about introducing additional cell types (e.g., osteoblasts).

Fourth, we chose to separate both long bone ends because the growth plates of the proximal and distal femur mature at different stages and under different processes, with the former undergoing complex chondrogenesis and osteogenesis patterns not evident in the latter (Stern et al., PLoS Biol. 2015; Cole et al., Bone.2013), and because proximal and distal growth plates show evidence of differential growth across a number of species, including humans. Thus, by separating growth plates, we are able to identify ATAC-seq peaks unique to each region and likely important for the marked distal femur elongation seen in humans.

Second, the manuscript still lacks a strong connection to chondrocyte biology as none of the links established by the data directly address (or explain) chondrocyte proliferation or hypertrophy and genes well known to regulate them positively or negatively (PTHrp, Indian hedgehog, FGFs, IGFs, etc.).

Thank you for raising this point regarding connections to known chondrocyte biology. In our revised manuscript, we have included ATAC-seq peaks at three well known growth plate genes: Col2a1 (encoding a collagen subunit), Pth1r (parathyroid hormone receptor), and Fgfr1 (fibroblast growth factor)(new Figure 1—figure supplement 1). Our GREAT analysis related to human phenotypes demonstrates many known connections to growth plate-related phenotypes (new Figure 1C). We now include GO Biological Processes from our GREAT analyses, which highlight terms such as “Cartilage condensation”, “Positive regulation of chondrocyte differentiation”, and “Chondrocyte proliferation”. Together, we believe these results highlight how the ATAC-seq peaks are capturing genes with fundamental roles in chondrocyte biology.

Third, although different genes have been identified as possibly being putative regulatory variants there are no experimental or computational evidence implicating the ATAC-seq peaks in regulating the nearby growth plate genes, nor altering transcriptional regulatory activity.

We thank the reviewers for their suggestions of including experimental validation. To this end, we have included a variety of experiments aimed at validating our experimental predictions at the CHSY1 locus.

1) First, using a luciferase reporter assay, we demonstrated that rs9920291 demonstrates allelic regulatory activity in the T/C-28a2 human chondrocyte cell line (new Figure 6B).

2) Using a human cell line (HEK-293FT), we demonstrate that there is allelic skew in CHSY1 expression. We chose to use HEK-293FT because it is heterozygous “C/T” for rs9920291. We had genotyped several human chondrocyte cell lines, and unfortunately they were all homozygous “C” for this variant (new Figure 6C). This is coupled with our analysis of GTEx data showing that rs9920291 is as an eQTL for CHSY1 in human transformed fibroblasts (new Figure 6—figure supplement 1).

3) Next, we used two separate pairs of CRISPR guide RNAs to delete the region surrounding rs9920291 in the T/C-28a2 cell line. We deleted either a large 135 bp region, or a smaller 17 bp region. Deletion of either region enhances CHSY1 expression, suggesting that rs9920291 marks a regulatory element with repressive activity (new Figure 6D).

4) Finally, we demonstrate that overexpression of HOXD13 can increase CHSY1 expression, providing evidence to support our computational predictions that HOXD13 may bind at the open chromatin region marked by rs9920291 and regulate CHSY1 expression (new Figure 6E and new Figure 6—figure supplement 1).

Together, we believe these new results provide strong evidence that our approach can pinpoint plausible candidate causal GWAS variants. Additionally, to our knowledge, this is one of the few height GWAS loci (along with GDF5 from our group) with experimental validation, again highlighting the novelty of our paper.

These issues underlie a perception that, independent of your Nature Genetics paper, there is a lack in this manuscript of mechanistic evidence that the epigenomic profiles can be linked to functional correlates. For example, there is a suggestion that rs9920291 disrupted HOXD13 binding, but there is no evidence that HOXD13 actually binds this site at the endogenous locus. As computationally predicted TF motifs are prevalent in the genome, it is hard to assess whether this SNP/motif overlap is biologically meaningful. This requirement for a mechanistic/functional connection may require additional work.

These major points will need to be addressed in the revision along with several other points imbedded in the three reviewers’ comments that have been included below in the comments, and should be considered carefully.

Reviewer #1:

[…] 1) The authors used E15.5 embryos from FVB mice to study the epigenetic regions for growth plate regulation. But inbred strains often exhibit numerous nucleotide variants; since C57BL6J mice are the most frequently utilized strain for in vitro and in vivo chondrocyte (and osteoblast) studies, have the authors determined whether the specific nucleotide sequence in the CHYS1 variant, rs9920291, is also present in B6 mouse chondrocytes?

Thank you for raising this important point regarding mouse strains and the potential for there to be genetic and epigenetic differences across strains. We chose the FVB strain for this study given the consistently large number of embryos that the FVB strain produces, which allowed us to increase the number of mice from which we could collect difficult to extract tissue. We also wish to clarify that rs9920291 is a human variant that, to our knowledge, is not present in mouse. We note that the reference “C” allele is the ancestral allele that is present on the mouse sequence.

2) Sex is a major determinant of height particularly early in life and during puberty. The early influences of sex steroids likely begins during late embryonic development; i.e. after E11.5 when SRY is expressed in the males, but not females. It would be important to fully understand the epigenetic landscape in respect to sex differences in vivo. Hence, the embryos should have been typed for sex, and both male and female embryos should be examined to more fully understand the transcriptional regulation of chondrocyte differentiation.

Thank you for raising this important point regarding sex of the mice used for analyses. We have addressed this point in our response to the Editor’s first comment (see above).

3) Can the authors justify why they did not follow up with functional studies on the CHYS1 variant rs9920291 which they identified as being a cis acting site for HOXD13.

As suggested by multiple reviewers, we have pursued functional studies of our candidate causal variant rs9920291 at CHSY1. We describe the new experiments we have performed in our response to the Editor’s third comment (see above).

4) Figure 6B of the E15.5 embryo showing the lacZ reporter from the ancestral GROW1 enhancer appears nearly identical to their Nature Genetics paper Figure 4B; Can the authors defend its use in this manuscript? Otherwise it should be deleted.

We thank the reviewer for this comment. We wish to clarify that the lacZ reporter used in this paper is from a stable line, and thus different from the original Nature Genetics paper. Given that the results are very similar, we have decided to remove the GDF5 experiment from our paper. Additionally, we have moved the epigenetic view of the GDF5 locus to the new Figure 6—figure supplement 2.

Reviewer #2:

[…] While the study is technically strong, it is deficient in skeletal cell and developmental biology, raising questions about the major conclusions reached and specific relevance to human height.

Thank you for your careful review of our manuscript. We are pleased that you found our paper “technically strong” and with “positive assets”. Below, we address the criticisms that you raise. We hope that our revised manuscript and new experiments address these gaps and make our manuscript suitable for publication.

The growth plate is certainly the key structure determining long bone and vertebral elongation and thus, body height. It is clearly established that the rates of skeletal growth are determined by chondrocyte proliferation (minor contribution) and chondrocyte hypertrophy (major contribution) in the growth plate. It is also well established that different growth plates within the same organism at any given time "grow" at very different rates, such as proximal and distal growth plates in radius and ulna. None of the links established by the data presented here directly address (or explain) chondrocyte proliferation or hypertrophy and genes known to regulate them positively or negatively (PTHrp, Indian hedgehog, FGFs, IGFs, etc.). The only links established are to general aspects of chondrocyte and skeletal biology, thus making it very difficult to understand how the links have anything to do with "height". The behavior of the growth plate is controlled by intrinsic mechanisms (see above) as well as systemic factors such as growth hormone. Again, none of the links related to this or other systemic cues.

Thank you for raising these points about linking our epigenetic data with biological insights into chondrocyte and skeletal biology. In our revision, we have included several improvements to strengthen the connection between our ATAC-seq data and chondrocyte biology. First, we revised our GREAT analysis by removing ATAC-seq peaks that are also seen in the mouse brain. By doing so, we are limiting ourselves to more chondrocyte-specific regions. When repeating the GREAT analysis, we find many chondrocyte-relevant gene sets, such as “Cartilage condensation”, “Positive regulation of chondrocyte differentiation”, and “Chondrocyte proliferation”. Second, we have included an additional figure (new Figure 1—figure supplement 1) where we show the landscape of ATACseq peaks near several genes with prior connections to chondrocyte biology – Ihh, Pth1r, and Fgfr1. We believe that these revisions help more directly connect our data to chondrocyte biology.

The authors made use of the proximal and distal portions of E15.5 mouse embryo femur and associated perichondrium as starting material for the epigenetic analysis. The tissues were digested into single cell suspension enzymatically before processing for ATAT-seq. There is no explanation of why this developmental stage was chosen and why both ends of the femur were considered. Also, the lengthy enzymatic digestion could have altered gene expression and chromatin configuration.

Thank you for raising this important point regarding our choice of using E15.5 embryos. We have addressed this point in our response to the Editor’s first comment (see above).

To establish stronger links, the authors considered the most compelling GWAS variants overlapping ATAT-seq peaks, finding 59 variants distributed over 26 loci. The two examples mentioned are CHSY1 and HOXD13 neither of which has established roles in growth plate function. Indeed, most evidence indicates that HOXD13 regulates autopod skeletal patterning and morphogenesis, but not rates of growth.

We explored the roles of a number of GWAS loci in the revised edited manuscript, highlighting ATAC-seq regions that overlap with previously published growth plate enhancers, reside near classic genes involved in growth plate biology, and display evidence of overlapping with human height GWAS variants. We also now have demonstrated evidence of their functional involvement in chondrocyte biology. For this manuscript, we focused on CHSY1 as one example relevant to both human height and human health, as it has previously published roles in chondrocyte biology and mutations in CHSY1 cause a variety of skeletal dysplasias (see main text). We identified HOXD13 as a potential transcriptional regulator of CHSY1 at the rs9920291 enhancer variantusing the UniProbe Database (see Materials and methods). While it is well-known that Hoxd13 plays important roles in autopod patterning, it is also expressed in growth plates (Reno et al., 2008) and has been known for some time (along with other Hox genes) to be growth promoting factor (Morgan and Tabin,’ 94; Yokouchi et al.,’ 95; Goff and Tabin,’ 97; Capecchi,’ 97; Papenbrock et al., 2000; Zhao and Potter, 2001). As Hox proteins act redundantly during limb bud development and growth plate function, deciphering the specific functions of Hoxd13 or its targets has been extremely difficult. However, our newly added data reveal that HOXD13 is expressed in chondrocytes and can increase CHSY1 expression when overexpressed in chondrocyte cell lines (new Figure 6E).

The authors make a strong case that the data can provide insights into the meaning and functional relevance of GWAS correlations, but I do not see how they provide any direct insights into what regulates human height. It may be worth considering an overall change in rationale for the study.

We thank the reviewer for raising this point regarding the framing of our work. However, we believe that our work provides many important insights into the regulation of human height. First, we identified 44 genes at height GWAS loci, where the gene is differentially expressed in the growth plate and where there is a height GWAS variant that overlaps ATAC-seq peaks. These 44 genes represent a subset of height-associated genes that likely regulate height at the growth plate. Second, we have identified several transcription factors with binding sites at these ATAC-seq peaks, providing additional insight into the transcriptional controls underlying height. Finally, as detailed above, we now provide functional experiments detailing how common regulatory variation at CHSY1 influences height. These results provide novel insights into the genes, regulatory circuitry, and mechanisms underlying control of height at the growth plate.

Reviewer #3:

[…] The authors proceed to identify height GWAS loci and specific SNPs that overlap femur open chromatin peaks, and provide potential hypotheses for how select SNPs (e.g. rs9920291) may act to influence human height. Unfortunately, despite claiming the identification of "compelling mechanisms for GWAS variants" (Abstract), the authors do not provide sufficient experimental or computational evidence in support, except for rs4911178 that references data in Capellini et al. 2017 (PMID 28671685). In summary, the authors present an important epigenomic dataset that can act as a first step for future study of height GWAS SNPs, but do not provide "strong evidence" (subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, first paragraph) for compelling mechanisms at any height GWAS loci beyond referencing their previously published paper.

We thank the reviewer for his/her careful review our paper and are pleased that the reviewer recognizes the value of our work. To strengthen our connections between ATAC-seq peaks at GWAS loci with biological mechanisms, we now provide extensive functional work at the CHSY1 locus. We have also made many other changes and performed new experiments as suggested by the reviewer.

1) Substantially fewer GWAS loci overlap ATAC-seq peaks when using the fine-mapped 95% credible variants, instead of all SNPs with r2>0.5 (e.g. 192 loci vs. 317 loci, subsection “Fine-mapped height variants overlap with femoral growth plate open chromatin regions” and subsection “Human height variants are enriched in femoral open chromatin regions”, last paragraph; also 26 loci vs. 46 loci, subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, first paragraph and subsection “A subset of human height variants in femoral open chromatin regions reside near genes differentially expressed in the growth plate”, second paragraph). Given that r2>0.5 is not a strict LD threshold (e.g. Wood et al. 2014 uses r2>0.8 for strict LD), is it still meaningful to report overlap statistics using the r2>0.5 set when the PICS 95% credible set is available?

We thank the reviewer for raising these important points regarding proxy SNPs at GWAS loci. We agree with the reviewer that r2>0.5 is not a universally agreed upon LD threshold. We have now included GoShifter analyses at several different r2thresholds ranging from 0.2 to 0.8 (new Figure 3—figure supplement 1). Among the various other LD thresholds tested, an r2 cutoff of 0.5 demonstrated the strongest enrichment (p=0.0059).

The reviewer makes an excellent point regarding whether we should just report PICS credible set overlaps. While the reviewer’s point is valid and well-taken, we have decided to be more inclusive in our approach and report both PICS overlaps and variants overlaps with r2>0.5 in our new Supplementary file 6. This choice was made particularly in light of potential issues with using PICS, as PICS may not be the best fine-mapping tool. Additionally, the height GWAS was imputed to HapMap (as compared to a more dense panel such as 1000 Genomes), limiting our ability to properly fine-map. We do wish to note that we were limited to using PICS for fine-mapping, as to our knowledge, PICS is the only fine-mapping software that doesn’t require densely imputed summary statistics, but only requires a set of lead variants and LD relationships.

2) The authors state "44 different genes have been identified as possibly being modulated by these putative regulatory variants". However, in the present manuscript the authors do not provide experimental or computational evidence implicating the ATAC-seq peaks in regulating the nearby growth plate genes, and provide limited evidence that any of these overlapping variants can alter transcriptional regulatory activity.

Thank you for this comment. We now provide evidence that differentially expressed growth plate genes that are also eQTLs genes in whole blood are enriched at ATAC-seq peaks (see response to point #6). Additionally, as detailed in our response to the Editor’s third comment, we now provide extensive experimental evidence for rs9920291 as a causal variant at the CHSY1 locus.

3) Following point #2, and to support the authors' claim of identifying "prime targets for functional testing", it would be helpful to know how many GWAS SNP overlaps with femur open chromatin regions the authors believe to be biologically meaningful versus coincidence.

We thank the reviewer for raising this important point. Unfortunately, without a “truth” set, it is difficult to evaluate what proportion of our nominated SNPs are likely to be biologically meaningful as opposed to coincidence. However, we do note that using both GoShifter and random matched GWAS loci, we detect very strong enrichments of our height GWAS SNPs at ATAC-seq loci, suggesting that a large proportion of our nominated SNPs are meaningful. We also point out the work of Ulirsch et al., Cell, 2016, which performed MPRA for red blood cell traits and demonstrated that a large proportion of variants nominated by similar approaches as used in our paper demonstrated functional activity in the MPRA assay. As an additional approach to demonstrate that our observed overlaps are not simply coincidence, we also performed GoShifter analyses for enrichments of GWAS associations of several traits (type 2 diabetes, total cholesterol, body mass index, schizophrenia, coronary artery disease, and age of menarche) at our ATAC-seq peaks. For most of these additional traits, we did not detect significant associations (p>0.05) of enrichment (new Figure 3—figure supplement 1). Lastly, as described elsewhere, we now provide new data at the CHSY1 locus to demonstrate the value of our approach.

4) Similar to point #2, it is premature to state "59 variants thus have strong evidence for being potentially causal variants". For the authors' example of the CHSY1 locus, while there are 4 SNPs in the PICS 95% set that overlap an ATAC-seq peak nearby CHSY1, Supplementary file 6 reveals there are 87 SNPs in the 95% credible set at the CHSY1 locus. Without more evidence, it is difficult to agree with the potential causal variant conclusion.

We recognize that it is difficult to determine what proportion of our prioritized variants will turn out to be causal. However, we believe that our ATAC-seq data is bolstering the evidence of these variants being causal. For example, at CHSY1, intersection with ATAC-seq data prioritized the 87 SNPs in the 95% credible set down to just four variants. In our revision, we now provide extensive functional evidence for rs9920291, providing further credence for our claim that our approach can nominate likely causal variants (new Figure 6).

5) Subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, first paragraph: The authors suggest rs9920291 disrupts HOXD13 binding, but do not provide evidence that HOXD13 actually binds this site at the endogenous locus. As computationally predicted TF motifs are prevalent in the genome, it is hard to assess whether this SNP/motif overlap is biologically meaningful.

We agree with the reviewer that predicted motifs are prevalent throughout the genome. While the HOXD13 binding prediction is based on experimental HOXD13 DNA-binding specificities (Hume et al. 2015), there are limited ChIP-seq datasets in relevant tissues to assess whether HOXD13 truly binds at rs9920291, in part because of the known difficulty of generating and using HOXD13-specific antibodies for ChIP-seq. We have toned down the language in the manuscript regarding this HOXD13 connection (subsection “Known and novel targets of human height variants in putative chondrocyte regulatory regions”, second paragraph). Additionally, we do now provide evidence that overexpression of HOXD13 increases expression of CHSY1 (new Figure 6E).

6) Wood et al. (2014) found that height GWAS loci were more likely to be eQTLs (in blood) for genes expressed in cartilage. Did the authors examine whether the 468 credible SNPs overlapping femur open chromatin peaks are enriched for eQTL signals? This analysis may strengthen the authors' mechanistic claims.

We thank the reviewer for this excellent suggestion. We have now included additional analyses, which confirm the observation in Wood et al., that height GWAS loci are more likely to be eQTLs in blood for genes expressed in cartilage. To do this, we used a set of eQTLs from whole blood from Westra et al., 2014 and subsetted for eQTLs genes that are differentially expressed in the growth plate. We then ran GoShifter on this set of corresponding eQTLs variants and indeed detected a significant enrichment when compared with our ATAC-seq peaks (p=0.0004). This is discussed in the third paragraph of the subsection “A subset of human height variants in femoral open chromatin regions reside near genes 261 differentially expressed in the growth plate”.

https://doi.org/10.7554/eLife.29329.037

Article and author information

Author details

  1. Michael Guo

    1. Broad Institute of MIT and Harvard, Cambridge, United States
    2. Division of Endocrinology, Boston Children’s Hospital, Harvard Medical School, Boston, United States
    3. Department of Genetics, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Zun Liu and Jessie Willen
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1357-6389
  2. Zun Liu

    Department of Human Evolutionary Biology, Harvard University, Cambridge, United States
    Contribution
    Formal analysis, Investigation
    Contributed equally with
    Michael Guo and Jessie Willen
    Competing interests
    No competing interests declared
  3. Jessie Willen

    Department of Human Evolutionary Biology, Harvard University, Cambridge, United States
    Contribution
    Data curation, Investigation, Methodology
    Contributed equally with
    Michael Guo and Zun Liu
    Competing interests
    No competing interests declared
  4. Cameron P Shaw

    Department of Human Evolutionary Biology, Harvard University, Cambridge, United States
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
  5. Daniel Richard

    1. Department of Human Evolutionary Biology, Harvard University, Cambridge, United States
    2. Department of Biology, University of Waterloo, Waterloo, Canada
    Contribution
    Formal analysis, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
  6. Evelyn Jagoda

    Department of Human Evolutionary Biology, Harvard University, Cambridge, United States
    Contribution
    Data curation, Investigation, Methodology
    Competing interests
    No competing interests declared
  7. Andrew C Doxey

    Department of Biology, University of Waterloo, Waterloo, Canada
    Contribution
    Methodology, Writing—original draft
    Competing interests
    No competing interests declared
  8. Joel Hirschhorn

    1. Broad Institute of MIT and Harvard, Cambridge, United States
    2. Division of Endocrinology, Boston Children’s Hospital, Harvard Medical School, Boston, United States
    3. Department of Genetics, Harvard Medical School, Boston, United States
    Contribution
    Supervision, Funding acquisition, Writing—review and editing
    Competing interests
    No competing interests declared
  9. Terence D Capellini

    1. Broad Institute of MIT and Harvard, Cambridge, United States
    2. Department of Human Evolutionary Biology, Harvard University, Cambridge, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    tcapellini@fas.harvard.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3842-8478

Funding

National Science Foundation (BCS-1518596)

  • Terence D Capellini

National Institutes of Health (R01 DK075787)

  • Joel Hirschhorn

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors would like to thank the following individuals for assistance with this manuscript: Dr. William Greenleaf and Dr. Jason Buenrostro for advice on ATAC-seq.; Dr. Claire Reardon and the entire Harvard University Bauer Core facility for assistance with next generation sequencing; Dr. Michelle Clamp at Harvard University Research Computing for assistance with ATAC-seq data; Dr. Shinsuke Ohba for providing access to Sox9 ChIP-seq datasets; Dr. Lin Wu for assistance with mouse stable line generation; Dr. Vincenzo Zappavigna for HOXD13 expression constructs; Drs. Li Zeng and Mary Goldring for T/C-28a2 chondrocytes; Dr. Pardis Sabeti for HEK-293FT cells. This work was funded in part by the Harvard University Milton Fund and NSF (BCS-1518596) to TDC, and National Institutes of Health (R01 DK075787) to JNH.

Ethics

Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Heath, and under protocols approved by the Institutional Animal Care and Use Committee (IACUC) at Harvard University (Protocol # 13-04-161). Euthanasia was conducted according to the recommendations and guidelines from the Panel on Euthanasia of the American Veterinary Medical Association and as outlined in Protocol #13-04-161.

Reviewing Editor

  1. Clifford J Rosen, Maine Medical Center Research Institute, United States

Publication history

  1. Received: June 8, 2017
  2. Accepted: November 7, 2017
  3. Version of Record published: December 5, 2017 (version 1)
  4. Version of Record updated: December 14, 2017 (version 2)

Copyright

© 2017, Guo et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 3,512
    Page views
  • 413
    Downloads
  • 16
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links