- Email: [email protected]
Genome-Wide Association Analysis of Longitudinal Bone Mineral Content Data From the Iowa Bone Development Study
ARTICLE IN PRESS Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health, vol. &, no. &, 111, 2019 Ó 2019 The International Society for Clinical Densitometry. Published by Elsevier Inc. 1094-6950/&:111/$36.00 https://doi.org/10.1016/j.jocd.2019.09.005
Genome-Wide Association Analysis of Longitudinal Bone Mineral Content Data From the Iowa Bone Development Study Camden P. Bay,*,1 Steven M. Levy,2,3 Kathleen F. Janz,4 Brian J. Smith,5,6 John R. Shaffer,7,8 Mary L. Marazita,9,10,11 and Trudy L. Burns12 1
Center for Clinical Investigation, Brigham & Women’s Hospital, Boston, MA, USA; 2 Department of Preventive and Community Dentistry, University of Iowa College of Dentistry, Iowa City, IA, USA; 3 Department of Epidemiology, University of Iowa College of Public Health, Iowa City, IA, USA; 4 Department of Health and Human Physiology, College of Liberal Arts and Sciences, University of Iowa, Iowa City, IA, USA; 5 Holden Comprehensive Cancer Center, University of Iowa Hospitals and Clinics, Iowa City, IA, USA; 6 Department of Biostatistics, University of Iowa College of Public Health, Iowa City, IA, USA; 7 Department of Human Genetics, Graduate School of Public Health, University of Pittsburgh, Pittsburgh, PA, USA; 8 Department of Oral Biology, School of Dental Medicine, University of Pittsburgh, Pittsburgh, PA, USA; 9 Center for Craniofacial and Dental Genetics, Department of Oral Biology, School of Dental Medicine, University of Pittsburgh, Pittsburgh, PA, USA; 10 Department of Human Genetics, Graduate School of Public Health, University of Pittsburgh, Pittsburgh, PA, USA; 11 Clinical and Translational Science, School of Medicine, University of Pittsburgh, Pittsburgh, PA, USA; and 12 Department of Epidemiology, University of Iowa College of Public Health, Iowa City, IA, USA
Abstract The foundation for osteoporosis risk is, in part, established during childhood, adolescence, and young adulthood, all periods of development when bone mass is acquired rapidly. The relative quantity of bone mass accrued is influenced by both lifestyle and genetic factors, although the genetic component is not yet well understood. The purpose of this study was to use a genome-wide association (GWA) analysis to discover single nucleotide polymorphisms (SNPs) associated with: (1) the sex-specific hip bone mineral content at approximately the age of 19 when the amount of bone accrued is near its peak; and (2) the sex-specific rate of hip bone mineral content accrual during the adolescent growth spurt. The Iowa Bone Development Study, a longitudinal cohort study exploring bone health in children, adolescents, and young adults was the source of data. From this cohort, n = 364 (190 females, 174 males) participants were included in GWA analyses to address (1) and n = 258 participants (125 females and 133 males) were included in GWA analyses to address (2). Twenty SNPS were detected having p < 1.0 £ 105. Of most biologic relevance were 2 suggestive SNPs (rs2051756 and rs2866908) detected in an intron of the DKK2 gene through the GWA analysis that explored peak bone mass in females. Key Words: Adolescent bone development; bone mass; DKK2; genome-wide association study; longitudinal.
Introduction A deficiency in the accrual of bone mass during childhood, adolescence, and young adulthood may increase the risk of osteoporosis later in life (1,2). Peak bone mass (the highest bone mass an individual acquires) is generally reached by the early to middle twenties in males and
Received 06/23/19; Revised 09/23/19; Accepted 09/24/19.
*Address correspondence to: Camden P. Bay, Center for Clinical Investigation, Brigham & Women’s Hospital, 221 Longwood Ave, Boston, MA 02115. E-mail: [email protected]
1
ARTICLE IN PRESS 2 by the late teens to early twenties in females, varying by bone site (3,4). Having a higher peak bone mass than one’s physically similar peers may lead to lower osteoporosis risk. The effects of diet and physical activity on bone accrual are well studied (2); however, until recently, study of the genetic component had been restricted to investigations of the genetic predictors of bone pathology in adults (57). In the last several years, studies have been published exploring genes known to be associated with adult bone fragility and childhood bone strength (811). Additionally, one recent study discovered genes specifically associated with childhood bone strength at the one-third distal radius, an area where childhood fractures are common (12). In 2016, a review article by Kemp et al was published, calling for childhood and adolescent bone acquisition research to more readily utilize genome-wide association (GWA) studies to better understand the genetics underlying healthy bone growth and ultimately improve pharmacological treatment of osteoporosis (13). In this study, the potential genetic contributors to bone development were examined for 2 phenotypes using GWA analyses. The first analysis sought to identify genetic variants that explain the peak bone mass of participants (from here on, referred to as the “peak bone mass” analysis) and the second sought to identify genetic variants that describe the rates of bone accrual during the adolescent growth spurt, which is the time of the largest proportion of bone mass accrual (from here on referred to as the “longitudinal” analysis). Additionally, gene £ physical activity interaction effects were assessed using the top SNPs from the peak bone mass GWA analysis (Online Resource 3).
Materials and Methods Data Source Participants in the Iowa Bone Development Study (IBDS) who had both genetic and longitudinal bone measurement data were included. The IBDS is a prospective cohort study based at the University of Iowa which has observed the development of healthy bones in over 600 participants with measurement occasions at the ages of 5, 8, 11, 13, 15, 17, and 19 years (1416). Most participants in this cohort are non-Hispanic White. Many members of this cohort were genotyped through the GENEVA Consortium (Gene Environment Association Studies), providing a total of approx 400 participants with bone measures and genotype data (17). Approval for the study was granted by the University of Iowa Institutional Review Board and written informed consent was given by parents of participants; assent was provided by the children.
Bay et al. and be genotyped for a GWA study exploring the genetic predictors of dental caries (17). Genotyping was performed with the Illumina Human610-Quadv1_B BeadChip (Illumina, San Diego, CA, USA) using DNA samples obtained from blood, mouthwash, buccal swabs, or saliva (17). Genotyping was performed by the Johns Hopkins Center for Inherited Disease Research (CIDR) and cleaning and quality assurance were carried out by the GENEVA Consortium Coordinating Center located at the University of Washington, resulting in 620,901 SNPs for 426 IBDS participants. For this study, participants with missing SNP data for greater than 10% of SNPs were removed (n = 4 participants). Individual SNPS were removed for being missing in greater than 10% of participants (n = 28,120 SNPs), having a HWE (Hardy-Weinberg equilibrium) exact test p value less than 0.0001 (n = 716 SNPs), or having a minor allele frequency (MAF) of less than 5% (n = 61,792 SNPs). Additionally, SNPs on the sex chromosomes were not included in the analyses. Cleaning resulted in 418 unrelated participants; the number of SNPs was 509,842. SNP genotypes were recoded to reflect an additive genetic effect (i.e., taking on the value of 0, 1, or 2 risk alleles, defined as the minor allele) (18,19). Data cleaning was performed by the investigators using PLINK (version 1.9) (20).
Control of Population Stratification Bias due to population stratification was reduced through the use of a subset of principal components describing genome-wide variation (21,22). Among all available SNPs (n = 509,842) a subset of SNPs with little or no linkage disequilibrium was selected (n = 91,752) through scanning windows of 50 SNPs at a time, stepping through the windows 5 SNPs at a time, and removing any SNPs within the window if they had a pairwise coefficient of determination greater than 0.2 by using the indeppairwise function in PLINK (version 1.9). The pca function was then used in PLINK to generate the top 20 principal components; the most informative principal components were then selected using a scree plot (23). A scatter plot of the first 2 principal components, which accounted for most of the variance, was created; participants outside of the central cloud were removed and the principal components were recalculated. This homogenization was to allow for a more fine-grain delineation of genetic variation. Only 2 participants were not White. The top 2 principal components were controlled for in all genetic association analyses. The scree plot, scatter plot of the first 2 principal components before participant removal, and scatter plot of the first 2 principal components after participant removal are displayed in Supplementary Figs. 13 in Online Resource 1.
Bone Measurement and Site Genotyping Individuals participating in the Iowa Fluoride Study (the parent study of the IBDS) were recruited to participate in
Hip bone measurements were obtained for participants during each measurement occasion using Hologic QDR 2000 DXA (ages 5 and 8) or Hologic QDR 4500 (ages
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data 1119) densitometers. QDR 2000 bone measurements were calibrated to QDR 4500 measurements based on an internal calibration study so as to be consistent over time; details of this calibration and quality control have been reported previously (24). The GWAS analysis was focused on hip bone mineral content (BMC). BMC was analyzed instead of the more commonly used bone mineral density (BMD) as it is a better surrogate measure of bone strength and allows for comparison among and within participants—due to the utilization of an areal form of density, BMD is subject to bias when comparing bones of different sizes (2527). In addition, the hip was chosen as the anatomical area of interest as it is one of the sites most likely to fracture when weakened by osteoporosis. In order to account for the influence of body size on BMC, height and weight were controlled for in all analyses. Bone size, which is highly correlated with BMC (28), was not controlled for due to the potential for over-adjustment.
Maturation To account for differences in the timing of maturation among participants, the age of participants was represented on a biological scale rather than a chronological one. This was accomplished by estimating the chronological age of peak height velocity for each participant as a participant-specific baseline in order to create a measure of time since peak maturation. The chronological age at peak height velocity was predicted using the sex-specific equations constructed by Mirwald, Baxter-Jones, Bailey, and Beunen (29). These equations require measurements of height, sitting height, and leg length near peak height velocity, these being the third measurement occasion in the IBDS (approximately age 11 years) for females and the fourth or fifth (approximately ages 13 and 15 years) for males, so participants without data near this required time point were removed (n = 68). As shown in multiple validation studies using data from the Fels Longitudinal Study and the Wroc»aw Growth Study, the use of measurements near peak height velocity is important for obtaining accurate predictions (3032). Using these predictions, chronological ages were transformed to biological ages for all remaining participants. On this scale, 0 represents the predicted time of peak height velocity. It should be noted that extrapolation of biological age to more than a few years beyond peak height velocity may be inappropriate (29).
Genome-Wide Association Analysis A total of 4 GWA analyses were performed, 2 (females only and males only) seeking to discover genetic predictors of approximate peak hip bone mass, and 2 exploring variation in the rate of bone mass accrual during the adolescent growth spurt. Separate GWA analyses were performed for males and females, because the set of SNPs related to bone development may vary by sex (3335). As a
3 supplementary analysis, the results of a GWA analysis pooling over sex are available in Online Resource 2 (Supplementary Table 1). All association testing was performed in R (version 3.0.2) using a cluster computer that allowed for parallel computing. Code is available in Online Resource 1 (Supplementary Code 12). GWA analysis results were assessed with Manhattan plots and QQ plots, constructed using the R package qqman (Supplementary Figs. 411 in Online Resource 1) (36). The UCSC Genome Browser was used for gene annotations (37). Q values, which originate from the false discovery paradigm, are reported, in addition to p values (38,39). Q values were calculated using the p values from within each GWA analysis, rather than p values across all GWA analyses.
Peak Bone Mass Analysis The association between approximate peak hip BMC (age 19) and a SNP genotype was accomplished by creating 2 linear regression models, 1 for females and 1 for males. The models can be seen below, where BMC19 represents the hip BMC at chronological age 19, Genotype can take on 3 values (0, 1, or 2 risk alleles), PC1 and PC2 represent the population stratification principal components, and i indexes the participant. The errors were assumed to be normal and homoscedastic. Missing response data were assumed to be missing completely at random. EðBMC19i Þ ¼ b0 þ b1 Genotypei þ b2 Heighti þ b3 Weighti þ b4 PC1i þ b5 PC2i Studentized residuals were used to check model assumptions and Cook’s distance and likelihood displacement were used to check for influential points. Linearity was assessed graphically and with the addition of polynomial terms in a model excluding Genotype. Testing of the genotype term in the 2 GWA analyses was performed assuming Student’s t-distribution.
Longitudinal Analysis To define the period of the adolescent growth spurt, linear splines with 2 knots were fit with biological age as the metameter for time and hip BMC as the response. Knot locations were varied within a reasonable biological age range; the decision of which knot locations to use was guided by content knowledge, plots of the data, and variation in the likelihood of the model used for fitting the splines. This procedure was performed for males and females separately and no accounting was made for statistical error in the final knot location decision. The region within the 2 knots was then determined to be the best representation of the adolescent growth spurt. A sensitivity analysis was performed to determine the effect of varying knot locations on the slope of the spline. The hypothesis of interest for this GWA analysis was that the rate of hip BMC accrual during the adolescent
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 4
Bay et al.
growth spurt varies by genotype. Testing of the hypothesis was accomplished through the use of a linear mixed model with a random intercept and slope (biological age) fit using restricted maximum likelihood (REML) in R with the lme4 package (40). The response (hip BMC in grams) was treated as continuous and assigned an identity link. Random effects and errors were assumed to be multivariate normal. All longitudinal missing response data were assumed to be missing completely at random. The linear mixed model specification can be seen below, where the BMCijs are the time varying BMC measurements during the adolescent growth spurt, Genotype can take on 3 values (0, 1, or 2 risk alleles), PC1 and PC2 represent the population stratification principal components, Height was specified as a polynomial of degree 2, i indexes the participant, and j is a measurement occasion index. The regression coefficient of interest was the interaction effect between Genotype and Biological Age which was tested assuming Student’s t-distribution with Satterthwaite degrees of freedom.
to determine if the regression coefficients in the above model were appropriate. Model assumptions were checked by using marginal and conditional studentized residuals and empirical best linear unbiased predictors. Mahalanobis distance, likelihood displacement, and Cook’s distance were used to check for influential observations and participants. Lastly, the marginal variancecovariance structure was assessed by using a semivariogram (4143).
Results Six-hundred thirty-one IBDS participants had at least 1 hip BMC measurement (3074 longitudinal measurements total). Sixty-eight were removed due to the inability to estimate peak height velocity, which resulted in a total of 563 participants. After merging with the genetic data and removing 18 participants determined to be genetically distant based on the principal components analysis, 364 participants had hip BMC and genetic data (2058 longitudinal measures). The 18 genetically distant individuals were selected for removal as they were not located within the central cluster of participants in the principal components plot (Supplementary Fig. 2 in Online Resource 1). It was desirable to use principal components for population stratification control that could account for a fine-grain delineation of genetic dissimilarity without being dominated by a small number of genetically distant individuals (Supplementary Fig. 3 in Online Resource 1). Three-hundred sixty-two participants self-reported as “White,” 1 as “Native American,” and 1 as “Hispanic.” There were 190 females and 174 males. Table 1 presents
E BMCij ¼ b0 þ b1 Biological Ageij þ b2 Genotypei þ b3 Genotypei Biological Ageij þ b4 Heightij þ b5 Heightij2 þ b6 Weightij þ b7 PC1i þ b8 PC2i Once fit, the 3 time-varying covariates were decomposed into among- and within-participant effects in order
Table 1 Sample Characteristics by Measurement Occasion for n = 364 Participants with Genetic Data and Bone Measurements. Chronological age (years) Measurement Sex (n) occasion 1 2 3 4 5 6 7
F (150) M (130) F (178) M (158) F (173) M (151) F (173) M (162) F (158) M (143) F (149) M (120) F (120) M (93)
Mean
SD
Median
5.3 5.3 8.7 8.7 11.3 11.3 13.3 13.3 15.3 15.3 17.5 17.5 19.7 19.7
0.4 0.5 0.6 0.7 0.3 0.3 0.4 0.4 0.3 0.3 0.4 0.4 0.6 0.7
5.2 5.2 8.6 8.6 11.2 11.1 13.1 13.2 15.2 15.3 17.5 17.6 19.7 19.6
Hip BMC (g)
Height (cm)
Weight (kg)
Mean SD Median Mean SD Median Mean SD 6.9 7.2 12.3 13.2 18.8 19.1 26.3 27.2 29.4 37.7 30.6 43.5 31.0 44.9
1.3 1.5 2.6 3.0 4.4 4.4 5.1 7.0 5.6 8.2 5.9 9.2 5.8 9.0
6.8 7.0 11.8 12.9 18.5 18.8 25.5 26.2 29.0 36.1 29.8 43.1 29.9 43.6
111.3 112.4 133.0 134.7 149.3 149.3 161.0 162.8 165.0 175.2 166.1 178.9 166.4 180.4
5.3 5.9 6.7 7.4 7.4 7.9 6.5 9.7 6.7 8.3 6.6 7.9 7.1 8.0
110.9 112.3 132.8 134.3 149.6 149.1 161.5 162.7 164.9 174.9 165.7 179.3 166.5 180.8
20.2 20.7 32.2 34.0 45.2 45.9 56.2 58.8 62.6 70.9 68.1 79.4 69.6 83.0
4.0 3.9 9.0 10.4 13.0 13.7 14.7 16.9 14.8 17.0 17.2 18.5 17.9 20.2
Median 19.4 20.0 30.3 31.3 42.4 43.0 52.5 55.9 57.9 67.9 63.0 76.1 64.1 78.1
Measurement occasion: 1 (5 years old), 2 (8), 3 (11), 4 (13), 5 (15), 6 (17), 7 (19). BMC, bone mineral content; SD, standard deviation. Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data
5
Fig. 1. Hip bone mineral content (BMC) by approximate age in female and male Iowa Bone Development Study participants (n = 364). summary statistics of all utilized covariates and hip BMC. Fig. 1 shows the increase in hip BMC over the 7 measurement occasions. Fig. 2 shows the relationship between biological age and chronological age by sex. The mean chronological age at the age 19 measurement occasion was 19.66 years (Standard Deviation = 0.60, Minimum = 18.62, Maximum = 22.10). Twohundred thirteen participants (120 females and 93 males) had measurements available at this time. Table 2 lists the top SNPs (p < 1.0 £ 105) from the peak bone mass GWA analysis by sex, and their estimated regression coefficients. Manhattan plots and QQ plots are available in Online Resource 1 (Supplementary Figs. 411). In the female-specific peak bone mass model, 2 suggestive SNPs, rs2051756 and rs2866908, were found within an intron of the Dickkopf-related protein 2 (DKK2) gene. An increase in the number of minor alleles (C and T,
respectively) at the 2 loci is associated with a decreased peak bone mass (3.18 [95% confidence interval {CI} 4.41 to 1.95] grams and 3.25 [95% CI 4.58 to 1.92] grams per risk allele, respectively—note that the confidence intervals are not adjusted for multiple testing). A LocusZoom plot of the 2 SNPs is displayed in Fig. 3 (44). In the male-specific peak bone mass model, rs2051756 and rs2866908 have raw p values of 0.39 and 0.33 and estimated effects of 1.01 and 1.36, respectively; in the model pooling over males and females, rs2051756 and rs2866908 have raw p values of 0.12 and 0.13 and estimated effects of 1.07 and 1.17, respectively. The adolescent growth spurt linear spline knot locations (beginning and ending of the spurt) were chosen for males to be at the biological ages of 2.0 and 2.0; for females, knot locations of 2.0 and 2.0 were also chosen. Fig. 4 displays the observed growth curves for female and
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 6
Bay et al.
Fig. 2. Scatter plot of chronological age by biological age for each participant (n = 364). male participants, in addition to the mean curve with 2 knots. Two-hundred fifty-nine participants had at least 2 measurements within the adolescent growth spurt (534 measurements total); 125 participants were females and 133 were males. During the GWA analysis, 24 SNPs did not converge in the model with only males; a different optimizer and number of REML iterations was attempted but failed to produce better results. Table 2 lists the top SNPs (p < 1.0 £ 105) by model and their estimated regression coefficients for variation in slope. No SNPs were found to be near genes with biological relevance, and thus, they are not discussed further. As described in Online Resource 2 (Supplementary Table 1), no evidence for association between bone mass and SNPs was detected in the GWA analyses pooling both sexes. Additionally, in the physical activity analyses, no evidence of physical activity £ discovered cross-sectional GWAS SNP interaction effects was found.
Discussion The purpose of this investigation was to discover SNPs associated with 2 representations of bone mass accrual: total hip bone mass measured around the age of 19 when the amount of bone accrued is approximately at its peak, and the rate of hip bone accrual during the adolescent
growth spurt. SNP discovery was performed using sexspecific GWA analyses based on regression models estimating longitudinal and cross-sectional relationships. Twenty SNPs were detected having a p value less than 1.0 £ 105. Gene annotations of the top SNPs from the peak bone mass (age 19) GWA analysis and the longitudinal (bone accrual rate during the adolescent growth spurt) GWA analysis did not lead to the identification of any genes of clear biological relevance other than DKK2 in the female-specific peak bone mass analysis. The DKK2 gene produces a protein that interacts with the Wnt signaling pathway and is associated with embryonic development. In an investigation of the Wnt signaling pathway, Chan et al found that the expression of DKK2 is notably higher in osteoarthritic osteoblasts than in normal ones (45). As rs2051756 alters an RXRA transcription factor binding motif and has chromatin marks suggesting enhancer activity for embryonic and induced pluripotent stem cells lines, there is evidence that it may be a functional variant. Additionally, rs2866908 is in high linkage disequilibrium with rs17584752 (r2 = 0.97, D’ = 0.99) which is an expression quantitative trait locus for the aforementioned DKK2 gene (46). Therefore, it is reasonable to suggest that variation in DKK2 either directly hinders or is a component in a family of genes that hinders peak hip bone mass acquisition in females.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data
7
Table 2 Sex-Specific Top SNPs (p < 1.0 £ 105) From the Peak Hip Bone Mineral Content GWA Analysis and the Longitudinal Hip Bone Mineral Content GWA Analysis. SNP
Chr:Position MA MA frequency Regression coefficient (SE)
Females only Peak BMC
p value
q value
rs2051756 rs7605368 rs12617390 rs2866908 rs2278583
4:108122994 2:42837844 2:42838899 4:108124957 2:42843096
C A G T T
0.19 0.23 0.24 0.17 0.23
3.18 (0.62) 2.76 (0.54) 2.65 (0.53) 3.25 (0.67) 2.55 (0.53)
1.25 £ 106 1.61 £ 106 2.42 £ 106 4.28 £ 106 5.49 £ 106
0.41 0.41 0.41 0.55 0.56
rs7522454 rs2042542
1:60756208 2:218444537
A G
0.30 0.33
0.52 (0.11) 0.54 (0.12)
7.19 £ 106 9.56 £ 106
1 1
rs12322558 rs17106937 rs3937886 rs2023499 rs885210
12:68366420 12:68352547 12:68373206 2:229715103 4:91794764
G G C C C
0.09 0.09 0.09 0.17 0.32
8.34 (1.61) 8.34 (1.61) 8.34 (1.61) 6.31 (1.22) 4.64 (0.92)
1.42 £ 106 1.42 £ 106 1.42 £ 106 1.55 £ 106 2.58 £ 106
0.20 0.20 0.20 0.20 0.26
rs1937920 rs2186172 rs281248 rs281281 rs3371 rs2398203 rs7083869 rs7720553
10:5151955 10:5165450 15:45454125 15:45448723 15:45437398 10:5166783 10:5234295 5:179556243
C G A G T T A A
0.25 0.25 0.32 0.32 0.31 0.24 0.20 0.13
0.98 (0.20) 0.98 (0.20) 0.80 (0.17) 0.80 (0.17) 0.80 (0.17) 0.95 (0.20) 0.99 (0.21) 1.10 (0.23)
2.87 £ 106 2.87 £ 106 5.70 £ 106 5.70 £ 106 5.70 £ 106 5.83 £ 106 6.58 £ 106 7.07 £ 106
0.45 0.45 0.45 0.45 0.45 0.45 0.45 0.45
Longitudinal BMC
Males only Peak BMC
Longitudinal BMC
Regression Coefficient represents the effect on hip BMC (g) for a 1-unit increase in the number of risk alleles in the peak BMC analysis and Regression Coefficient represents the change in the slope of hip BMC (g) over time for a 1-unit increase in the number of risk alleles. The q value for a hypothesis test can be interpreted as the proportion of false positives among all hypothesis test rejections that would be expected if the test’s q value is treated as a threshold for statistical significance. Position is based on the March 2006 HG18 human assembly. Abbr: BMC, bone mineral content; Chr, Chromosome; MA, Minor allele.
Other than this investigation, the DKK2 gene has not been identified in research exploring genetic associations with BMC and BMD in children or adults (712,4750). These studies, some based on single cohorts and some meta-analyses, explored a variety of bone sites and age groups. For example, while observing male and female children (mean age 11.4, standard deviation 4.5) separately in a GWAS design, Chesi et al discovered and replicated associations between variants in 3 genes (CPED1, MIR31HG, and MTAP) and BMC or BMD at the distal radius (12). CPED1 is part of a well-known pathway associated with BMD and the latter 2 genes are novel findings. With total-body BMD in children (mean age 6.2 years, standard deviation 0.3) as an outcome for their GWA, Medina-Gomez et al discovered and replicated associations with variants in the WNT16 (important in bone
formation) and C7orf58 genes (47). Knockout mice provided further evidence for the contribution of WNT16 and C7orf58 to BMD. Timpson et al performed a GWA in children (mean age 9.9 years, standard deviation not reported) and found suggestive evidence that variants in SP7, known to be associated with osteoblast differentiation, may be related to total-body BMD (48). None of the genes of interest detected by these studies or the cited meta-analyses were replicated within the current investigation (Supplementary Table 4 in Online Resource 4). This may be due to this investigation’s exclusive focus on longitudinal change and peak bone mass near the beginning of adulthood and the decision to use BMC instead of the more commonly studied BMD. Additionally, the small study size may have prevented replication of previous results. The detection of the DKK2 gene is also
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 8
Bay et al.
Fig. 3. Plot of base pair location by log10 p value generated in LocusZoom of the top two SNPs (rs2051756 and rs2866908) from the female peak bone mass GWA analysis.
subject to limitations. Given the high q values of the 2 SNPS (0.41 and 0.55) associated with DKK2, and more importantly, the inability to perform a replication study due to a current lack of comparable data, a false positive finding is possible and the biological relevance may be coincidental. Additionally, this gene, if it is truly related, may also be describing body size, even though a body size influence was partially controlled for with the inclusion of carefully specified height and weight covariates in the models. It should also be noted that since the longitudinal GWA p values depended heavily on the correct specification of the correlation structure, p values may have been affected by incorrect specification; however, diagnostic measures were carefully performed to prevent this. Lastly, results may have been affected by measurement error resulting from the use of DXA scans. The design of the IBDS was a major strength of this investigation and allowed for the precise tracking of the entire process of major bone development in children, adolescents, and young adults. Only with this information was it possible to measure when, on average, the adolescent growth spurt was initiated and ended. Additionally, due to the consistent recording of height and weight, biological age could be determined for most participants without sacrificing much precision.
This study specifically focused on the effect that genetic variation has on peak hip BMC and the rate of hip BMC development during the adolescent growth spurt. The focus of future work could be on these specific developmental periods and on this specific bone site, but the importance of the periods of bone development before and after the adolescent growth spurt and the possibility that certain genetic effects may be specific to bone sites is also of great interest.
Acknowledgments We would like to thank the participants of the Iowa Bone Development Study and their families; with their patience and effort, many significant insights into the growth of bone have and will continue to be made. Additionally, we would like to thank Elena Letuchy for preparing and assisting in the interpretation of the bone mineral content and physical activity data. This investigation was supported in part by NIH grants R01-DE09551, R01-DE12101, M01-RR00059, UL1-RR024979, U01DE018903, Dr. Levy’s WrightBush-Shreves Endowed Professor Fund, and Dr. Shaffer’s University of Pittsburgh/UPMC Competitive Medical Research Fund.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data
9
Fig. 4. Observed hip bone mineral content (BMC) spaghetti plots for females and males and mean hip BMC based on the chosen knot locations (2.0 and 2.0 biological age years).
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 10
Supplementary materials Supplementary material associated with this article can be found in the online version at doi:10.1016/j.jocd.2019.09.005.
References 1. Heaney RP, Abrams S, Dawson-Hughes B, et al. 2000 Peak bone mass. Osteoporos Int 11:985–1009. 2. Weaver CM, Gordon CM, Janz KF, et al. 2016 The National Osteoporosis Foundation’s position statement on peak bone mass development and lifestyle factors: a systematic review and implementation recommendations. Osteoporos Int 27: 1281–1386. 3. Baxter-Jones AD, Faulkner RA, Forwood MR, et al. 2011 Bone mineral accrual from 8 to 30 years of age: an estimation of peak bone mass. J Bone Miner Res 26:1729–1739. 4. Nguyen TV, Maynard LM, Towne B, et al. 2001 Sex differences in bone mass acquisition during growth: the Fels Longitudinal Study. J Clin Densitom 4:147–157. 5. Estrada K, Styrkarsdottir U, Evangelou E, et al. 2012 Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat Genet 44:491–501. 6. Ralston SH, Uitterlinden AG. 2010 Genetics of osteoporosis. Endocr Rev 31:629–662. 7. Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, et al. 2008 Multiple genetic loci for bone mineral density and fractures. N Engl J Med 358:2355–2365. 8. Mitchell JA, Chesi A, Elci O, et al. 2016 Genetic risk scores implicated in adult bone fragility associate with pediatric bone density. J Bone Miner Res 31:789–795. 9. Mitchell JA, Chesi A, Elci O, et al. 2015 Genetics of bone mass in childhood and adolescence: effects of sex and maturation interactions. J Bone Miner Res 30:1676–1683. 10. Mitchell JA, Chesi A, Cousminer DL, et al. 2017 Multidimensional bone density phenotyping reveals new insights into genetic regulation of the pediatric skeleton. J Bone Miner Res 33:812–821. 11. Chesi A, Mitchell JA, Kalkwarf HJ, et al. 2017 A genomewide association study identifies two sexspecific loci, at SPTB and IZUMO3, influencing pediatric bone mineral density at multiple skeletal sites. J Bone Miner Res 32: 1274–1281. 12. Chesi A, Mitchell JA, Kalkwarf HJ, et al. 2015 A trans-ethnic genome-wide association study identifies gender-specific loci influencing pediatric aBMD and BMC at the distal radius. Hum Mol Genet 24:5053–5059. 13. Kemp JP, Medina-Gomez C, Tobias JH, et al. 2016 The case for genome-wide association studies of bone acquisition in paediatric and adolescent populations. Bonekey Rep 5:796. 14. Janz KF, Gilmore JM, Levy SM, et al. 2007 Physical activity and femoral neck bone strength during childhood: the Iowa Bone Development Study. Bone 41:216–222. 15. Janz KF, Letuchy EM, Burns TL, et al. 2014 Objectively measured physical activity trajectories predict adolescent bone strength: Iowa Bone Development Study. Br J Sports Med 48:1032–1036. 16. Levy SM, Warren JJ, Davis CS. 2001 Patterns of fluoride intake from birth to 36 months. J Public Health Dent 61: 70–77. 17. Shaffer JR, Wang X, Feingold E, et al. 2011 Genome-wide association scan for childhood caries implicates novel genes. J Dent Res 90:1457–1462.
Bay et al. 18. Bush WS, Moore JH. 2012 Chapter 11: genome-wide association studies. PLoS Comput Biol 8:1–11. 19. Lettre G, Lange C, Hirschhorn JN. 2007 Genetic model testing and statistical power in population-based association studies of quantitative traits. Genet Epidemiol 31:358–362. 20. Chang CC, Chow CC, Tellier LC, et al. 2015 Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4:1–16. 21. Price AL, Patterson NJ, Plenge RM, et al. 2006 Principal components analysis corrects for stratification in genomewide association studies. Nat Genet 38:904–909. 22. Ziegler A, Konig IR. 2010 A statistical approach to genetic epidemiology. Weinheim, Germany: Wiley-Blackwell; 2010. 23. Patterson N, Price AL, Reich D. 2006 Population structure and eigen analysis. PLoS Genet 2:2074–2093. 24. Francis SL, Letuchy EM, Levy SM, Janz KF. 2014 Sustained effects of physical activity on bone health: Iowa Bone Development Study. Bone 63:95–100. 25. Heaney RP. 2003 Bone mineral content, not bone mineral density, is the correct bone measure for growth studies. Am J Clin Nutr 78:350–351. 26. Heaney RP. 2005 BMD: the problem. Osteoporos Int 16: 1013–1015. 27. Carter DR, Bouxsein ML, Marcus R. 1992 New approaches for interpreting projected bone densitometry data. J Bone Miner Res 7:137–145. 28. Budek AZ, Hoppe C, Ingstrup H, et al. 2007 Dietary protein intake and bone mineral content in adolescents-The Copenhagen Cohort Study. Osteoporos Int 18:1661–1667. 29. Mirwald RL, Baxter-Jones AD, Bailey DA, Beunen GP. 2002 An assessment of maturity from anthropometric measurements. Med Sci Sports Exerc 34:689–694. 30. Malina RM, Koziel SM. 2014 Validation of maturity offset in a longitudinal sample of Polish boys. J Sports Sci 32: 424–437. 31. Malina RM, Koziel SM. 2014 Validation of maturity offset in a longitudinal sample of Polish girls. J Sports Sci 32: 1374–1382. 32. Malina RM, Choh AC, Czerwinski SA, Chumlea WC. 2016 Validation of maturity offset in the Fels longitudinal study. Pediatr Exerc Sci 28:439–455. 33. Randall JC, Winkler TW, Kutalik Z, et al. 2013 Sex-stratified genome-wide association studies including 270,000 individuals show sexual dimorphism in genetic loci for anthropometric traits. PLoS Genet 9:1–19. 34. Liu CT, Estrada K, Yerges-Armstrong LM, et al. 2012 Assessment of gene-by-sex interaction effect on bone mineral density. J Bone Miner Res 27:2051–2064. 35. Nielson CM, Klein RF, Orwoll ES. 2012 Sex and the single nucleotide polymorphism: exploring the genetic causes of skeletal sex differences. J Bone Miner Res 27:2047–2050. 36. Turner SD2014qqman: an R package for visualizing GWAS results using Q-Q and Manhattan plots. R statistical software package. 37. Speir ML, Zweig AS, Rosenbloom KR, et al. 2016 The UCSC Genome Browser database: 2016 update. Nucleic Acids Res 44:D717–D725. 38. Storey JD, Tibshirani R. 2003 Statistical significance for genomewide studies. Proc Natl Acad Sci USA 100: 9440–9445. 39. Storey JD, Bass AJ, Dabney A, Robinson D2015qvalue: Qvalue estimation for false discovery rate control. R statistical software package. 40. Bates D, Machler M, Bolker B, Walker S. 2015 Fitting linear mixed-effects models using lme4. J Stat Softw 67:1–48.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data 41. Galecki A, Burzykowski T. 2013 Linear mixed-effects models using R: a step-by-step approach. New York, NY: Springer; 2013. 42. Fitzmaurice GM, Laird NM, Ware JH. 2011 Applied longitudinal analysis. Hoboken, NJ: Wiley; 2011. 43. Nobre JS, da Motta Singer J. 2007 Residual analysis for linear mixed models. Biom J 49:863–875. 44. Pruim RJ, Welch RP, Sanna S, et al. 2010 LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics 26:2336–2337. 45. Chan TF, Couchourel D, Abed E, et al. 2011 Elevated Dickkopf-2 levels contribute to the abnormal phenotype of human osteoarthritic osteoblasts. J Bone Miner Res 26:1399–1410. 46. Ward LD, Kellis M. 2016 HaploReg v4: systematic mining of putative causal variants, cell types, regulators and target genes for human complex traits and disease. Nucleic Acids Res 44:D877–D881.
11 47. Medina-Gomez C, Kemp JP, Estrada K, et al. 2012 Metaanalysis of genome-wide scans for total body BMD in children and adults reveals allelic heterogeneity and agespecific effects at the WNT16 locus. PLoS Genet 8:1–14. 48. Timpson NJ, Tobias JH, Richards JB, et al. 2009 Common variants in the region around Osterix are associated with bone mineral density and growth in childhood. Hum Mol Genet 18: 1510-1507. 49. Kemp JP, Medina-Gomez C, Estrada K, et al. 2014 Phenotypic dissection of bone mineral density reveals skeletal site specificity and facilitates the identification of novel loci in the genetic regulation of bone mass attainment. PLoS Genet 10:1–18. 50. Rivadeneira F, Styrk arsdottir U, Estrada K, et al. 2009 Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat Genet 41:1199–1206.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
Genome-Wide Association Analysis of Longitudinal Bone Mineral Content Data From the Iowa Bone Development Study Camden P. Bay,*,1 Steven M. Levy,2,3 Kathleen F. Janz,4 Brian J. Smith,5,6 John R. Shaffer,7,8 Mary L. Marazita,9,10,11 and Trudy L. Burns12 1
Center for Clinical Investigation, Brigham & Women’s Hospital, Boston, MA, USA; 2 Department of Preventive and Community Dentistry, University of Iowa College of Dentistry, Iowa City, IA, USA; 3 Department of Epidemiology, University of Iowa College of Public Health, Iowa City, IA, USA; 4 Department of Health and Human Physiology, College of Liberal Arts and Sciences, University of Iowa, Iowa City, IA, USA; 5 Holden Comprehensive Cancer Center, University of Iowa Hospitals and Clinics, Iowa City, IA, USA; 6 Department of Biostatistics, University of Iowa College of Public Health, Iowa City, IA, USA; 7 Department of Human Genetics, Graduate School of Public Health, University of Pittsburgh, Pittsburgh, PA, USA; 8 Department of Oral Biology, School of Dental Medicine, University of Pittsburgh, Pittsburgh, PA, USA; 9 Center for Craniofacial and Dental Genetics, Department of Oral Biology, School of Dental Medicine, University of Pittsburgh, Pittsburgh, PA, USA; 10 Department of Human Genetics, Graduate School of Public Health, University of Pittsburgh, Pittsburgh, PA, USA; 11 Clinical and Translational Science, School of Medicine, University of Pittsburgh, Pittsburgh, PA, USA; and 12 Department of Epidemiology, University of Iowa College of Public Health, Iowa City, IA, USA
Abstract The foundation for osteoporosis risk is, in part, established during childhood, adolescence, and young adulthood, all periods of development when bone mass is acquired rapidly. The relative quantity of bone mass accrued is influenced by both lifestyle and genetic factors, although the genetic component is not yet well understood. The purpose of this study was to use a genome-wide association (GWA) analysis to discover single nucleotide polymorphisms (SNPs) associated with: (1) the sex-specific hip bone mineral content at approximately the age of 19 when the amount of bone accrued is near its peak; and (2) the sex-specific rate of hip bone mineral content accrual during the adolescent growth spurt. The Iowa Bone Development Study, a longitudinal cohort study exploring bone health in children, adolescents, and young adults was the source of data. From this cohort, n = 364 (190 females, 174 males) participants were included in GWA analyses to address (1) and n = 258 participants (125 females and 133 males) were included in GWA analyses to address (2). Twenty SNPS were detected having p < 1.0 £ 105. Of most biologic relevance were 2 suggestive SNPs (rs2051756 and rs2866908) detected in an intron of the DKK2 gene through the GWA analysis that explored peak bone mass in females. Key Words: Adolescent bone development; bone mass; DKK2; genome-wide association study; longitudinal.
Introduction A deficiency in the accrual of bone mass during childhood, adolescence, and young adulthood may increase the risk of osteoporosis later in life (1,2). Peak bone mass (the highest bone mass an individual acquires) is generally reached by the early to middle twenties in males and
Received 06/23/19; Revised 09/23/19; Accepted 09/24/19.
*Address correspondence to: Camden P. Bay, Center for Clinical Investigation, Brigham & Women’s Hospital, 221 Longwood Ave, Boston, MA 02115. E-mail: [email protected]
1
ARTICLE IN PRESS 2 by the late teens to early twenties in females, varying by bone site (3,4). Having a higher peak bone mass than one’s physically similar peers may lead to lower osteoporosis risk. The effects of diet and physical activity on bone accrual are well studied (2); however, until recently, study of the genetic component had been restricted to investigations of the genetic predictors of bone pathology in adults (57). In the last several years, studies have been published exploring genes known to be associated with adult bone fragility and childhood bone strength (811). Additionally, one recent study discovered genes specifically associated with childhood bone strength at the one-third distal radius, an area where childhood fractures are common (12). In 2016, a review article by Kemp et al was published, calling for childhood and adolescent bone acquisition research to more readily utilize genome-wide association (GWA) studies to better understand the genetics underlying healthy bone growth and ultimately improve pharmacological treatment of osteoporosis (13). In this study, the potential genetic contributors to bone development were examined for 2 phenotypes using GWA analyses. The first analysis sought to identify genetic variants that explain the peak bone mass of participants (from here on, referred to as the “peak bone mass” analysis) and the second sought to identify genetic variants that describe the rates of bone accrual during the adolescent growth spurt, which is the time of the largest proportion of bone mass accrual (from here on referred to as the “longitudinal” analysis). Additionally, gene £ physical activity interaction effects were assessed using the top SNPs from the peak bone mass GWA analysis (Online Resource 3).
Materials and Methods Data Source Participants in the Iowa Bone Development Study (IBDS) who had both genetic and longitudinal bone measurement data were included. The IBDS is a prospective cohort study based at the University of Iowa which has observed the development of healthy bones in over 600 participants with measurement occasions at the ages of 5, 8, 11, 13, 15, 17, and 19 years (1416). Most participants in this cohort are non-Hispanic White. Many members of this cohort were genotyped through the GENEVA Consortium (Gene Environment Association Studies), providing a total of approx 400 participants with bone measures and genotype data (17). Approval for the study was granted by the University of Iowa Institutional Review Board and written informed consent was given by parents of participants; assent was provided by the children.
Bay et al. and be genotyped for a GWA study exploring the genetic predictors of dental caries (17). Genotyping was performed with the Illumina Human610-Quadv1_B BeadChip (Illumina, San Diego, CA, USA) using DNA samples obtained from blood, mouthwash, buccal swabs, or saliva (17). Genotyping was performed by the Johns Hopkins Center for Inherited Disease Research (CIDR) and cleaning and quality assurance were carried out by the GENEVA Consortium Coordinating Center located at the University of Washington, resulting in 620,901 SNPs for 426 IBDS participants. For this study, participants with missing SNP data for greater than 10% of SNPs were removed (n = 4 participants). Individual SNPS were removed for being missing in greater than 10% of participants (n = 28,120 SNPs), having a HWE (Hardy-Weinberg equilibrium) exact test p value less than 0.0001 (n = 716 SNPs), or having a minor allele frequency (MAF) of less than 5% (n = 61,792 SNPs). Additionally, SNPs on the sex chromosomes were not included in the analyses. Cleaning resulted in 418 unrelated participants; the number of SNPs was 509,842. SNP genotypes were recoded to reflect an additive genetic effect (i.e., taking on the value of 0, 1, or 2 risk alleles, defined as the minor allele) (18,19). Data cleaning was performed by the investigators using PLINK (version 1.9) (20).
Control of Population Stratification Bias due to population stratification was reduced through the use of a subset of principal components describing genome-wide variation (21,22). Among all available SNPs (n = 509,842) a subset of SNPs with little or no linkage disequilibrium was selected (n = 91,752) through scanning windows of 50 SNPs at a time, stepping through the windows 5 SNPs at a time, and removing any SNPs within the window if they had a pairwise coefficient of determination greater than 0.2 by using the indeppairwise function in PLINK (version 1.9). The pca function was then used in PLINK to generate the top 20 principal components; the most informative principal components were then selected using a scree plot (23). A scatter plot of the first 2 principal components, which accounted for most of the variance, was created; participants outside of the central cloud were removed and the principal components were recalculated. This homogenization was to allow for a more fine-grain delineation of genetic variation. Only 2 participants were not White. The top 2 principal components were controlled for in all genetic association analyses. The scree plot, scatter plot of the first 2 principal components before participant removal, and scatter plot of the first 2 principal components after participant removal are displayed in Supplementary Figs. 13 in Online Resource 1.
Bone Measurement and Site Genotyping Individuals participating in the Iowa Fluoride Study (the parent study of the IBDS) were recruited to participate in
Hip bone measurements were obtained for participants during each measurement occasion using Hologic QDR 2000 DXA (ages 5 and 8) or Hologic QDR 4500 (ages
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data 1119) densitometers. QDR 2000 bone measurements were calibrated to QDR 4500 measurements based on an internal calibration study so as to be consistent over time; details of this calibration and quality control have been reported previously (24). The GWAS analysis was focused on hip bone mineral content (BMC). BMC was analyzed instead of the more commonly used bone mineral density (BMD) as it is a better surrogate measure of bone strength and allows for comparison among and within participants—due to the utilization of an areal form of density, BMD is subject to bias when comparing bones of different sizes (2527). In addition, the hip was chosen as the anatomical area of interest as it is one of the sites most likely to fracture when weakened by osteoporosis. In order to account for the influence of body size on BMC, height and weight were controlled for in all analyses. Bone size, which is highly correlated with BMC (28), was not controlled for due to the potential for over-adjustment.
Maturation To account for differences in the timing of maturation among participants, the age of participants was represented on a biological scale rather than a chronological one. This was accomplished by estimating the chronological age of peak height velocity for each participant as a participant-specific baseline in order to create a measure of time since peak maturation. The chronological age at peak height velocity was predicted using the sex-specific equations constructed by Mirwald, Baxter-Jones, Bailey, and Beunen (29). These equations require measurements of height, sitting height, and leg length near peak height velocity, these being the third measurement occasion in the IBDS (approximately age 11 years) for females and the fourth or fifth (approximately ages 13 and 15 years) for males, so participants without data near this required time point were removed (n = 68). As shown in multiple validation studies using data from the Fels Longitudinal Study and the Wroc»aw Growth Study, the use of measurements near peak height velocity is important for obtaining accurate predictions (3032). Using these predictions, chronological ages were transformed to biological ages for all remaining participants. On this scale, 0 represents the predicted time of peak height velocity. It should be noted that extrapolation of biological age to more than a few years beyond peak height velocity may be inappropriate (29).
Genome-Wide Association Analysis A total of 4 GWA analyses were performed, 2 (females only and males only) seeking to discover genetic predictors of approximate peak hip bone mass, and 2 exploring variation in the rate of bone mass accrual during the adolescent growth spurt. Separate GWA analyses were performed for males and females, because the set of SNPs related to bone development may vary by sex (3335). As a
3 supplementary analysis, the results of a GWA analysis pooling over sex are available in Online Resource 2 (Supplementary Table 1). All association testing was performed in R (version 3.0.2) using a cluster computer that allowed for parallel computing. Code is available in Online Resource 1 (Supplementary Code 12). GWA analysis results were assessed with Manhattan plots and QQ plots, constructed using the R package qqman (Supplementary Figs. 411 in Online Resource 1) (36). The UCSC Genome Browser was used for gene annotations (37). Q values, which originate from the false discovery paradigm, are reported, in addition to p values (38,39). Q values were calculated using the p values from within each GWA analysis, rather than p values across all GWA analyses.
Peak Bone Mass Analysis The association between approximate peak hip BMC (age 19) and a SNP genotype was accomplished by creating 2 linear regression models, 1 for females and 1 for males. The models can be seen below, where BMC19 represents the hip BMC at chronological age 19, Genotype can take on 3 values (0, 1, or 2 risk alleles), PC1 and PC2 represent the population stratification principal components, and i indexes the participant. The errors were assumed to be normal and homoscedastic. Missing response data were assumed to be missing completely at random. EðBMC19i Þ ¼ b0 þ b1 Genotypei þ b2 Heighti þ b3 Weighti þ b4 PC1i þ b5 PC2i Studentized residuals were used to check model assumptions and Cook’s distance and likelihood displacement were used to check for influential points. Linearity was assessed graphically and with the addition of polynomial terms in a model excluding Genotype. Testing of the genotype term in the 2 GWA analyses was performed assuming Student’s t-distribution.
Longitudinal Analysis To define the period of the adolescent growth spurt, linear splines with 2 knots were fit with biological age as the metameter for time and hip BMC as the response. Knot locations were varied within a reasonable biological age range; the decision of which knot locations to use was guided by content knowledge, plots of the data, and variation in the likelihood of the model used for fitting the splines. This procedure was performed for males and females separately and no accounting was made for statistical error in the final knot location decision. The region within the 2 knots was then determined to be the best representation of the adolescent growth spurt. A sensitivity analysis was performed to determine the effect of varying knot locations on the slope of the spline. The hypothesis of interest for this GWA analysis was that the rate of hip BMC accrual during the adolescent
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 4
Bay et al.
growth spurt varies by genotype. Testing of the hypothesis was accomplished through the use of a linear mixed model with a random intercept and slope (biological age) fit using restricted maximum likelihood (REML) in R with the lme4 package (40). The response (hip BMC in grams) was treated as continuous and assigned an identity link. Random effects and errors were assumed to be multivariate normal. All longitudinal missing response data were assumed to be missing completely at random. The linear mixed model specification can be seen below, where the BMCijs are the time varying BMC measurements during the adolescent growth spurt, Genotype can take on 3 values (0, 1, or 2 risk alleles), PC1 and PC2 represent the population stratification principal components, Height was specified as a polynomial of degree 2, i indexes the participant, and j is a measurement occasion index. The regression coefficient of interest was the interaction effect between Genotype and Biological Age which was tested assuming Student’s t-distribution with Satterthwaite degrees of freedom.
to determine if the regression coefficients in the above model were appropriate. Model assumptions were checked by using marginal and conditional studentized residuals and empirical best linear unbiased predictors. Mahalanobis distance, likelihood displacement, and Cook’s distance were used to check for influential observations and participants. Lastly, the marginal variancecovariance structure was assessed by using a semivariogram (4143).
Results Six-hundred thirty-one IBDS participants had at least 1 hip BMC measurement (3074 longitudinal measurements total). Sixty-eight were removed due to the inability to estimate peak height velocity, which resulted in a total of 563 participants. After merging with the genetic data and removing 18 participants determined to be genetically distant based on the principal components analysis, 364 participants had hip BMC and genetic data (2058 longitudinal measures). The 18 genetically distant individuals were selected for removal as they were not located within the central cluster of participants in the principal components plot (Supplementary Fig. 2 in Online Resource 1). It was desirable to use principal components for population stratification control that could account for a fine-grain delineation of genetic dissimilarity without being dominated by a small number of genetically distant individuals (Supplementary Fig. 3 in Online Resource 1). Three-hundred sixty-two participants self-reported as “White,” 1 as “Native American,” and 1 as “Hispanic.” There were 190 females and 174 males. Table 1 presents
E BMCij ¼ b0 þ b1 Biological Ageij þ b2 Genotypei þ b3 Genotypei Biological Ageij þ b4 Heightij þ b5 Heightij2 þ b6 Weightij þ b7 PC1i þ b8 PC2i Once fit, the 3 time-varying covariates were decomposed into among- and within-participant effects in order
Table 1 Sample Characteristics by Measurement Occasion for n = 364 Participants with Genetic Data and Bone Measurements. Chronological age (years) Measurement Sex (n) occasion 1 2 3 4 5 6 7
F (150) M (130) F (178) M (158) F (173) M (151) F (173) M (162) F (158) M (143) F (149) M (120) F (120) M (93)
Mean
SD
Median
5.3 5.3 8.7 8.7 11.3 11.3 13.3 13.3 15.3 15.3 17.5 17.5 19.7 19.7
0.4 0.5 0.6 0.7 0.3 0.3 0.4 0.4 0.3 0.3 0.4 0.4 0.6 0.7
5.2 5.2 8.6 8.6 11.2 11.1 13.1 13.2 15.2 15.3 17.5 17.6 19.7 19.6
Hip BMC (g)
Height (cm)
Weight (kg)
Mean SD Median Mean SD Median Mean SD 6.9 7.2 12.3 13.2 18.8 19.1 26.3 27.2 29.4 37.7 30.6 43.5 31.0 44.9
1.3 1.5 2.6 3.0 4.4 4.4 5.1 7.0 5.6 8.2 5.9 9.2 5.8 9.0
6.8 7.0 11.8 12.9 18.5 18.8 25.5 26.2 29.0 36.1 29.8 43.1 29.9 43.6
111.3 112.4 133.0 134.7 149.3 149.3 161.0 162.8 165.0 175.2 166.1 178.9 166.4 180.4
5.3 5.9 6.7 7.4 7.4 7.9 6.5 9.7 6.7 8.3 6.6 7.9 7.1 8.0
110.9 112.3 132.8 134.3 149.6 149.1 161.5 162.7 164.9 174.9 165.7 179.3 166.5 180.8
20.2 20.7 32.2 34.0 45.2 45.9 56.2 58.8 62.6 70.9 68.1 79.4 69.6 83.0
4.0 3.9 9.0 10.4 13.0 13.7 14.7 16.9 14.8 17.0 17.2 18.5 17.9 20.2
Median 19.4 20.0 30.3 31.3 42.4 43.0 52.5 55.9 57.9 67.9 63.0 76.1 64.1 78.1
Measurement occasion: 1 (5 years old), 2 (8), 3 (11), 4 (13), 5 (15), 6 (17), 7 (19). BMC, bone mineral content; SD, standard deviation. Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data
5
Fig. 1. Hip bone mineral content (BMC) by approximate age in female and male Iowa Bone Development Study participants (n = 364). summary statistics of all utilized covariates and hip BMC. Fig. 1 shows the increase in hip BMC over the 7 measurement occasions. Fig. 2 shows the relationship between biological age and chronological age by sex. The mean chronological age at the age 19 measurement occasion was 19.66 years (Standard Deviation = 0.60, Minimum = 18.62, Maximum = 22.10). Twohundred thirteen participants (120 females and 93 males) had measurements available at this time. Table 2 lists the top SNPs (p < 1.0 £ 105) from the peak bone mass GWA analysis by sex, and their estimated regression coefficients. Manhattan plots and QQ plots are available in Online Resource 1 (Supplementary Figs. 411). In the female-specific peak bone mass model, 2 suggestive SNPs, rs2051756 and rs2866908, were found within an intron of the Dickkopf-related protein 2 (DKK2) gene. An increase in the number of minor alleles (C and T,
respectively) at the 2 loci is associated with a decreased peak bone mass (3.18 [95% confidence interval {CI} 4.41 to 1.95] grams and 3.25 [95% CI 4.58 to 1.92] grams per risk allele, respectively—note that the confidence intervals are not adjusted for multiple testing). A LocusZoom plot of the 2 SNPs is displayed in Fig. 3 (44). In the male-specific peak bone mass model, rs2051756 and rs2866908 have raw p values of 0.39 and 0.33 and estimated effects of 1.01 and 1.36, respectively; in the model pooling over males and females, rs2051756 and rs2866908 have raw p values of 0.12 and 0.13 and estimated effects of 1.07 and 1.17, respectively. The adolescent growth spurt linear spline knot locations (beginning and ending of the spurt) were chosen for males to be at the biological ages of 2.0 and 2.0; for females, knot locations of 2.0 and 2.0 were also chosen. Fig. 4 displays the observed growth curves for female and
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 6
Bay et al.
Fig. 2. Scatter plot of chronological age by biological age for each participant (n = 364). male participants, in addition to the mean curve with 2 knots. Two-hundred fifty-nine participants had at least 2 measurements within the adolescent growth spurt (534 measurements total); 125 participants were females and 133 were males. During the GWA analysis, 24 SNPs did not converge in the model with only males; a different optimizer and number of REML iterations was attempted but failed to produce better results. Table 2 lists the top SNPs (p < 1.0 £ 105) by model and their estimated regression coefficients for variation in slope. No SNPs were found to be near genes with biological relevance, and thus, they are not discussed further. As described in Online Resource 2 (Supplementary Table 1), no evidence for association between bone mass and SNPs was detected in the GWA analyses pooling both sexes. Additionally, in the physical activity analyses, no evidence of physical activity £ discovered cross-sectional GWAS SNP interaction effects was found.
Discussion The purpose of this investigation was to discover SNPs associated with 2 representations of bone mass accrual: total hip bone mass measured around the age of 19 when the amount of bone accrued is approximately at its peak, and the rate of hip bone accrual during the adolescent
growth spurt. SNP discovery was performed using sexspecific GWA analyses based on regression models estimating longitudinal and cross-sectional relationships. Twenty SNPs were detected having a p value less than 1.0 £ 105. Gene annotations of the top SNPs from the peak bone mass (age 19) GWA analysis and the longitudinal (bone accrual rate during the adolescent growth spurt) GWA analysis did not lead to the identification of any genes of clear biological relevance other than DKK2 in the female-specific peak bone mass analysis. The DKK2 gene produces a protein that interacts with the Wnt signaling pathway and is associated with embryonic development. In an investigation of the Wnt signaling pathway, Chan et al found that the expression of DKK2 is notably higher in osteoarthritic osteoblasts than in normal ones (45). As rs2051756 alters an RXRA transcription factor binding motif and has chromatin marks suggesting enhancer activity for embryonic and induced pluripotent stem cells lines, there is evidence that it may be a functional variant. Additionally, rs2866908 is in high linkage disequilibrium with rs17584752 (r2 = 0.97, D’ = 0.99) which is an expression quantitative trait locus for the aforementioned DKK2 gene (46). Therefore, it is reasonable to suggest that variation in DKK2 either directly hinders or is a component in a family of genes that hinders peak hip bone mass acquisition in females.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data
7
Table 2 Sex-Specific Top SNPs (p < 1.0 £ 105) From the Peak Hip Bone Mineral Content GWA Analysis and the Longitudinal Hip Bone Mineral Content GWA Analysis. SNP
Chr:Position MA MA frequency Regression coefficient (SE)
Females only Peak BMC
p value
q value
rs2051756 rs7605368 rs12617390 rs2866908 rs2278583
4:108122994 2:42837844 2:42838899 4:108124957 2:42843096
C A G T T
0.19 0.23 0.24 0.17 0.23
3.18 (0.62) 2.76 (0.54) 2.65 (0.53) 3.25 (0.67) 2.55 (0.53)
1.25 £ 106 1.61 £ 106 2.42 £ 106 4.28 £ 106 5.49 £ 106
0.41 0.41 0.41 0.55 0.56
rs7522454 rs2042542
1:60756208 2:218444537
A G
0.30 0.33
0.52 (0.11) 0.54 (0.12)
7.19 £ 106 9.56 £ 106
1 1
rs12322558 rs17106937 rs3937886 rs2023499 rs885210
12:68366420 12:68352547 12:68373206 2:229715103 4:91794764
G G C C C
0.09 0.09 0.09 0.17 0.32
8.34 (1.61) 8.34 (1.61) 8.34 (1.61) 6.31 (1.22) 4.64 (0.92)
1.42 £ 106 1.42 £ 106 1.42 £ 106 1.55 £ 106 2.58 £ 106
0.20 0.20 0.20 0.20 0.26
rs1937920 rs2186172 rs281248 rs281281 rs3371 rs2398203 rs7083869 rs7720553
10:5151955 10:5165450 15:45454125 15:45448723 15:45437398 10:5166783 10:5234295 5:179556243
C G A G T T A A
0.25 0.25 0.32 0.32 0.31 0.24 0.20 0.13
0.98 (0.20) 0.98 (0.20) 0.80 (0.17) 0.80 (0.17) 0.80 (0.17) 0.95 (0.20) 0.99 (0.21) 1.10 (0.23)
2.87 £ 106 2.87 £ 106 5.70 £ 106 5.70 £ 106 5.70 £ 106 5.83 £ 106 6.58 £ 106 7.07 £ 106
0.45 0.45 0.45 0.45 0.45 0.45 0.45 0.45
Longitudinal BMC
Males only Peak BMC
Longitudinal BMC
Regression Coefficient represents the effect on hip BMC (g) for a 1-unit increase in the number of risk alleles in the peak BMC analysis and Regression Coefficient represents the change in the slope of hip BMC (g) over time for a 1-unit increase in the number of risk alleles. The q value for a hypothesis test can be interpreted as the proportion of false positives among all hypothesis test rejections that would be expected if the test’s q value is treated as a threshold for statistical significance. Position is based on the March 2006 HG18 human assembly. Abbr: BMC, bone mineral content; Chr, Chromosome; MA, Minor allele.
Other than this investigation, the DKK2 gene has not been identified in research exploring genetic associations with BMC and BMD in children or adults (712,4750). These studies, some based on single cohorts and some meta-analyses, explored a variety of bone sites and age groups. For example, while observing male and female children (mean age 11.4, standard deviation 4.5) separately in a GWAS design, Chesi et al discovered and replicated associations between variants in 3 genes (CPED1, MIR31HG, and MTAP) and BMC or BMD at the distal radius (12). CPED1 is part of a well-known pathway associated with BMD and the latter 2 genes are novel findings. With total-body BMD in children (mean age 6.2 years, standard deviation 0.3) as an outcome for their GWA, Medina-Gomez et al discovered and replicated associations with variants in the WNT16 (important in bone
formation) and C7orf58 genes (47). Knockout mice provided further evidence for the contribution of WNT16 and C7orf58 to BMD. Timpson et al performed a GWA in children (mean age 9.9 years, standard deviation not reported) and found suggestive evidence that variants in SP7, known to be associated with osteoblast differentiation, may be related to total-body BMD (48). None of the genes of interest detected by these studies or the cited meta-analyses were replicated within the current investigation (Supplementary Table 4 in Online Resource 4). This may be due to this investigation’s exclusive focus on longitudinal change and peak bone mass near the beginning of adulthood and the decision to use BMC instead of the more commonly studied BMD. Additionally, the small study size may have prevented replication of previous results. The detection of the DKK2 gene is also
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 8
Bay et al.
Fig. 3. Plot of base pair location by log10 p value generated in LocusZoom of the top two SNPs (rs2051756 and rs2866908) from the female peak bone mass GWA analysis.
subject to limitations. Given the high q values of the 2 SNPS (0.41 and 0.55) associated with DKK2, and more importantly, the inability to perform a replication study due to a current lack of comparable data, a false positive finding is possible and the biological relevance may be coincidental. Additionally, this gene, if it is truly related, may also be describing body size, even though a body size influence was partially controlled for with the inclusion of carefully specified height and weight covariates in the models. It should also be noted that since the longitudinal GWA p values depended heavily on the correct specification of the correlation structure, p values may have been affected by incorrect specification; however, diagnostic measures were carefully performed to prevent this. Lastly, results may have been affected by measurement error resulting from the use of DXA scans. The design of the IBDS was a major strength of this investigation and allowed for the precise tracking of the entire process of major bone development in children, adolescents, and young adults. Only with this information was it possible to measure when, on average, the adolescent growth spurt was initiated and ended. Additionally, due to the consistent recording of height and weight, biological age could be determined for most participants without sacrificing much precision.
This study specifically focused on the effect that genetic variation has on peak hip BMC and the rate of hip BMC development during the adolescent growth spurt. The focus of future work could be on these specific developmental periods and on this specific bone site, but the importance of the periods of bone development before and after the adolescent growth spurt and the possibility that certain genetic effects may be specific to bone sites is also of great interest.
Acknowledgments We would like to thank the participants of the Iowa Bone Development Study and their families; with their patience and effort, many significant insights into the growth of bone have and will continue to be made. Additionally, we would like to thank Elena Letuchy for preparing and assisting in the interpretation of the bone mineral content and physical activity data. This investigation was supported in part by NIH grants R01-DE09551, R01-DE12101, M01-RR00059, UL1-RR024979, U01DE018903, Dr. Levy’s WrightBush-Shreves Endowed Professor Fund, and Dr. Shaffer’s University of Pittsburgh/UPMC Competitive Medical Research Fund.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data
9
Fig. 4. Observed hip bone mineral content (BMC) spaghetti plots for females and males and mean hip BMC based on the chosen knot locations (2.0 and 2.0 biological age years).
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS 10
Supplementary materials Supplementary material associated with this article can be found in the online version at doi:10.1016/j.jocd.2019.09.005.
References 1. Heaney RP, Abrams S, Dawson-Hughes B, et al. 2000 Peak bone mass. Osteoporos Int 11:985–1009. 2. Weaver CM, Gordon CM, Janz KF, et al. 2016 The National Osteoporosis Foundation’s position statement on peak bone mass development and lifestyle factors: a systematic review and implementation recommendations. Osteoporos Int 27: 1281–1386. 3. Baxter-Jones AD, Faulkner RA, Forwood MR, et al. 2011 Bone mineral accrual from 8 to 30 years of age: an estimation of peak bone mass. J Bone Miner Res 26:1729–1739. 4. Nguyen TV, Maynard LM, Towne B, et al. 2001 Sex differences in bone mass acquisition during growth: the Fels Longitudinal Study. J Clin Densitom 4:147–157. 5. Estrada K, Styrkarsdottir U, Evangelou E, et al. 2012 Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat Genet 44:491–501. 6. Ralston SH, Uitterlinden AG. 2010 Genetics of osteoporosis. Endocr Rev 31:629–662. 7. Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, et al. 2008 Multiple genetic loci for bone mineral density and fractures. N Engl J Med 358:2355–2365. 8. Mitchell JA, Chesi A, Elci O, et al. 2016 Genetic risk scores implicated in adult bone fragility associate with pediatric bone density. J Bone Miner Res 31:789–795. 9. Mitchell JA, Chesi A, Elci O, et al. 2015 Genetics of bone mass in childhood and adolescence: effects of sex and maturation interactions. J Bone Miner Res 30:1676–1683. 10. Mitchell JA, Chesi A, Cousminer DL, et al. 2017 Multidimensional bone density phenotyping reveals new insights into genetic regulation of the pediatric skeleton. J Bone Miner Res 33:812–821. 11. Chesi A, Mitchell JA, Kalkwarf HJ, et al. 2017 A genomewide association study identifies two sexspecific loci, at SPTB and IZUMO3, influencing pediatric bone mineral density at multiple skeletal sites. J Bone Miner Res 32: 1274–1281. 12. Chesi A, Mitchell JA, Kalkwarf HJ, et al. 2015 A trans-ethnic genome-wide association study identifies gender-specific loci influencing pediatric aBMD and BMC at the distal radius. Hum Mol Genet 24:5053–5059. 13. Kemp JP, Medina-Gomez C, Tobias JH, et al. 2016 The case for genome-wide association studies of bone acquisition in paediatric and adolescent populations. Bonekey Rep 5:796. 14. Janz KF, Gilmore JM, Levy SM, et al. 2007 Physical activity and femoral neck bone strength during childhood: the Iowa Bone Development Study. Bone 41:216–222. 15. Janz KF, Letuchy EM, Burns TL, et al. 2014 Objectively measured physical activity trajectories predict adolescent bone strength: Iowa Bone Development Study. Br J Sports Med 48:1032–1036. 16. Levy SM, Warren JJ, Davis CS. 2001 Patterns of fluoride intake from birth to 36 months. J Public Health Dent 61: 70–77. 17. Shaffer JR, Wang X, Feingold E, et al. 2011 Genome-wide association scan for childhood caries implicates novel genes. J Dent Res 90:1457–1462.
Bay et al. 18. Bush WS, Moore JH. 2012 Chapter 11: genome-wide association studies. PLoS Comput Biol 8:1–11. 19. Lettre G, Lange C, Hirschhorn JN. 2007 Genetic model testing and statistical power in population-based association studies of quantitative traits. Genet Epidemiol 31:358–362. 20. Chang CC, Chow CC, Tellier LC, et al. 2015 Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4:1–16. 21. Price AL, Patterson NJ, Plenge RM, et al. 2006 Principal components analysis corrects for stratification in genomewide association studies. Nat Genet 38:904–909. 22. Ziegler A, Konig IR. 2010 A statistical approach to genetic epidemiology. Weinheim, Germany: Wiley-Blackwell; 2010. 23. Patterson N, Price AL, Reich D. 2006 Population structure and eigen analysis. PLoS Genet 2:2074–2093. 24. Francis SL, Letuchy EM, Levy SM, Janz KF. 2014 Sustained effects of physical activity on bone health: Iowa Bone Development Study. Bone 63:95–100. 25. Heaney RP. 2003 Bone mineral content, not bone mineral density, is the correct bone measure for growth studies. Am J Clin Nutr 78:350–351. 26. Heaney RP. 2005 BMD: the problem. Osteoporos Int 16: 1013–1015. 27. Carter DR, Bouxsein ML, Marcus R. 1992 New approaches for interpreting projected bone densitometry data. J Bone Miner Res 7:137–145. 28. Budek AZ, Hoppe C, Ingstrup H, et al. 2007 Dietary protein intake and bone mineral content in adolescents-The Copenhagen Cohort Study. Osteoporos Int 18:1661–1667. 29. Mirwald RL, Baxter-Jones AD, Bailey DA, Beunen GP. 2002 An assessment of maturity from anthropometric measurements. Med Sci Sports Exerc 34:689–694. 30. Malina RM, Koziel SM. 2014 Validation of maturity offset in a longitudinal sample of Polish boys. J Sports Sci 32: 424–437. 31. Malina RM, Koziel SM. 2014 Validation of maturity offset in a longitudinal sample of Polish girls. J Sports Sci 32: 1374–1382. 32. Malina RM, Choh AC, Czerwinski SA, Chumlea WC. 2016 Validation of maturity offset in the Fels longitudinal study. Pediatr Exerc Sci 28:439–455. 33. Randall JC, Winkler TW, Kutalik Z, et al. 2013 Sex-stratified genome-wide association studies including 270,000 individuals show sexual dimorphism in genetic loci for anthropometric traits. PLoS Genet 9:1–19. 34. Liu CT, Estrada K, Yerges-Armstrong LM, et al. 2012 Assessment of gene-by-sex interaction effect on bone mineral density. J Bone Miner Res 27:2051–2064. 35. Nielson CM, Klein RF, Orwoll ES. 2012 Sex and the single nucleotide polymorphism: exploring the genetic causes of skeletal sex differences. J Bone Miner Res 27:2047–2050. 36. Turner SD2014qqman: an R package for visualizing GWAS results using Q-Q and Manhattan plots. R statistical software package. 37. Speir ML, Zweig AS, Rosenbloom KR, et al. 2016 The UCSC Genome Browser database: 2016 update. Nucleic Acids Res 44:D717–D725. 38. Storey JD, Tibshirani R. 2003 Statistical significance for genomewide studies. Proc Natl Acad Sci USA 100: 9440–9445. 39. Storey JD, Bass AJ, Dabney A, Robinson D2015qvalue: Qvalue estimation for false discovery rate control. R statistical software package. 40. Bates D, Machler M, Bolker B, Walker S. 2015 Fitting linear mixed-effects models using lme4. J Stat Softw 67:1–48.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019
ARTICLE IN PRESS Analysis of Longitudinal Bone Mineral Content Data 41. Galecki A, Burzykowski T. 2013 Linear mixed-effects models using R: a step-by-step approach. New York, NY: Springer; 2013. 42. Fitzmaurice GM, Laird NM, Ware JH. 2011 Applied longitudinal analysis. Hoboken, NJ: Wiley; 2011. 43. Nobre JS, da Motta Singer J. 2007 Residual analysis for linear mixed models. Biom J 49:863–875. 44. Pruim RJ, Welch RP, Sanna S, et al. 2010 LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics 26:2336–2337. 45. Chan TF, Couchourel D, Abed E, et al. 2011 Elevated Dickkopf-2 levels contribute to the abnormal phenotype of human osteoarthritic osteoblasts. J Bone Miner Res 26:1399–1410. 46. Ward LD, Kellis M. 2016 HaploReg v4: systematic mining of putative causal variants, cell types, regulators and target genes for human complex traits and disease. Nucleic Acids Res 44:D877–D881.
11 47. Medina-Gomez C, Kemp JP, Estrada K, et al. 2012 Metaanalysis of genome-wide scans for total body BMD in children and adults reveals allelic heterogeneity and agespecific effects at the WNT16 locus. PLoS Genet 8:1–14. 48. Timpson NJ, Tobias JH, Richards JB, et al. 2009 Common variants in the region around Osterix are associated with bone mineral density and growth in childhood. Hum Mol Genet 18: 1510-1507. 49. Kemp JP, Medina-Gomez C, Estrada K, et al. 2014 Phenotypic dissection of bone mineral density reveals skeletal site specificity and facilitates the identification of novel loci in the genetic regulation of bone mass attainment. PLoS Genet 10:1–18. 50. Rivadeneira F, Styrk arsdottir U, Estrada K, et al. 2009 Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat Genet 41:1199–1206.
Journal of Clinical Densitometry: Assessment & Management of Musculoskeletal Health
Volume 00, 2019