Association Between Age at Puberty and Bone Accrual From 10 to 25 Years of Age

Key Points Question Is puberty timing associated with growth-related bone accrual up to adulthood? Findings In this cohort study of 6389 participants who underwent repeated bone density scans from ages 10 to 25 years, later puberty was associated with persistently lower bone mineral density, despite some catch-up during puberty. Meaning People with older pubertal age should be advised on how to maximize bone density and minimize its decrease in later life to help prevent fracture and osteoporosis.


eMethods. Supplemental Methods.
APHV was estimated using Superimposition by Translation and Rotation (SITAR) growth curve analysis (1,2). SITAR is a shape invariant model consisting of a natural cubic spline mean curve of height vs age. SITAR fits this curve with three person-specific random-effects that reflect individual differences from the mean curve in size (random height intercept for person-specific shift up or down in the spline curve), velocity (random age scaling that reflects differences in the duration of the growth spurt in individuals) and tempo (random age intercept reflecting individual differences in the timing of the growth spurt). SITAR models were fitted separately for males and females with 5 degrees of freedom. APHV estimation was restricted to individuals with height measurements between ages 5 and 20 years and at least one height measure after age 9 years. Data cleaning similar to previous ALSPAC analyses was performed to remove data points with velocity exceeding 4 SDs and standardized residuals exceeding 4 in absolute value (2,3,4), leaving a total of 5223 males (with 55 752 height measurements) and 5322 females (with 59 318 height measurements) with an estimate of APHV.
The final SITAR models explained 95.7% and 95.6% of variance in height in males and females respectively. APHV in years was estimated by transforming the random age intercept that reflects individual differences in the tempo (timing) of the growth spurt.
We subsequently investigated the relationship between APHV and whole-body BMD and BMC accrual using two sets of mixed-effects linear spline regression analyses. Linear spline models facilitate the summarising of nonlinear bone accrual through a series of linear splines joined at knot points (5,6). They estimate mean and person-specific bone values at the first assessment and mean and person-specific rates of accrual during each of the different growth periods. Smoothing functions, prior subject knowledge and fit statistics were used to identify the best fitting models (5, 6) and subsequently knot point selection was based on the mean age at DXA scans. Models were fitted using maximum likelihood and included random intercepts and random age slopes. A 1 st order autoregressive correlation structure was used to account for temporal autocorrelation, and a fixed variance structure was used to model heteroskedasticity. In both sets of linear spline models, we examined the relationship between APHV and rates of bone accrual by including APHV-by-age splines interaction terms and included sex-by-age splines interactions terms to test for statistical evidence of different bone trajectories for males and females. All models were adjusted for birth weight, ethnicity, maternal education and childhood BMI and energy intake. The linear spline models provided a good fit to the data, Pearson's correlations between the observed and predicted BMD and BMC values (without random effects) were 0.84 and 0.83.
In the first analyses, age was centred at APHV in order to estimate rates of bone accrual relative to APHV; i.e. from before puberty up to adulthood. Knot points were placed at the mean ages of the 12, 14, and 16-year scans which equated to summarising rate of accrual from 7.7 up to 0.9 years prior to APHV, 0.9 year prior up to 2.1 years following APHV, 2.1 to 3.7 years following APHV, and 3.7 to 16 years after APHV. The mean predicted trajectories relative to APHV were plotted for males and females. For the second analyses, we excluded individuals with accrual from before APHV in order to investigate associations between APHV and subsequent rate of bone accrual up to adulthood. Here, knot points were placed at the mean ages of the 14 (centred at age 14), 16, and 18-year assessments, which then estimated associations with rate of accrual up to age 14, between 14-16, 16-18 and 18-25 years. To assess if difference in bone outcomes persisted into adulthood, predicted bone values at age 25 were regressed on APHV (after adjusting for predicted values at younger ages) and compared to associations with predicted values at age 14. To visualise associations between APHV and bone accrual, predicted mean trajectories were plotted for individuals in the 10 th , 50 th , and 90 th sex-specific APHV percentiles.
In sensitivity analyses, we refitted models after re-estimating APHV from measured heights only (3) and using age at menarche (reported prospectively by the mother in years and months) instead of APHV as an alternative marker of age at puberty in females. For our secondary outcomes, we investigated associations between APHV and subsequent site-specific BMD and BMC accrual from whole-body DXA, and with total hip and femoral neck BMD accrual from hip DXA. Given that hips were only scanned at three different ages (14, 18 and 25), we fitted hip models using a single knot point placed at the mean ages of the 18-year assessment. The predicted mean trajectories were plotted for all secondary outcomes. All analyses were performed in R (R Foundation for Statistical Computing, Vienna); SITAR models were fitted with the 'sitar' package and linear spline models with the 'nlme' package.
The linear spline models reduced bias due to missing bone data by including all measures under the missing at random assumption (i.e. the probability of missing bone data assumed to be related to measured factors but unrelated to unmeasured factors). eTable 1 shows differences between those included in the analysis with complete and partial DXA data, and suggests that the probability of missing at least one bone measure was related to the model covariates, and thus possibly missing at random. Those missing all DXA measures were excluded from the analysis. eTable 2 shows differences between those included in the linear spline models and excluded due to missing all DXA data, and shows that compared with the participants in this study, those eTable APHV-years [mean (SD)] 13.5 (0.9) 13.5 (0.7) 11.6 (0.8) 11.7 (0.8)