Faster ticking rate of the epigenetic clock is associated with faster pubertal development in girls

ABSTRACT

Epigenetic age is an indicator of biological aging, capturing the impact of environmental and behavioral influences across time on cellular function. Deviance between epigenetic age and chronological age (AgeAccel) is a predictor of health. Pubertal timing has similarly been associated with cancer risk and mortality rate among females. We examined the association between AgeAccel and pubertal timing and adolescent breast composition in the longitudinal Growth and Obesity Cohort Study. AgeAccel was estimated in whole blood using the Horvath method at breast Tanner 2 (B2) and 4 (B4). Total breast volume, absolute fibro-glandular volume (FGV), and %FGV were evaluated at B4 using dual X-ray absorptiometry. The impact of AgeAccel (mean: 0; SD: 3.78) across puberty on the time to breast development (thelarche), menarche, and pubertal tempo (thelarche to menarche) was estimated using accelerated failure time models; generalized estimating equations were used to evaluate associations with breast density. A five-year increase in average adolescent AgeAccel was associated with a significant decrease in time to menarche [hazard ratio (HR): 1.37; 95% confidence interval (CI): 1.04, 1.80] adjusting for birth weight, maternal pre-pregnancy body mass index, maternal height, maternal education, B2 height, fat percentage, and cell composition. AgeAccel displayed a stronger inverse association with pubertal tempo (HR: 1.48; 95% CI: 1.10, 1.99). A five-year increase in AgeAccel was associated with 5% greater %FGV, adjusting for B4 percent body fat, and maternal traits (95% CI: 1.01, 1.10). Our study provides unique insight into the influence of AgeAccel on pubertal development in girls, which may have implications for adult health.

KEYWORDS:

• Puberty
• thelarche
• menarche
• breast density
• epigenetic clock
• DNA methylation

Introduction

There is accumulating evidence that DNA methylation (DNAm) can be used to estimate biological age [1,2], capturing the impact of behavioral and environmental influences on gene regulation across the life-course, such as in utero exposures and adult diet [3,4]. Accordingly, epigenetic age has been shown to be a predictor of future health, independent of chronological age. Increased deviance between epigenetic age and chronological age (epigenetic age acceleration; AgeAccel) has been associated with an increased hazard of mortality in several large prospective cohorts after adjusting for established determinants of mortality rate [5,6]. Similarly, greater AgeAccel has been associated with increased risk of overall cancer incidence and death [7,8], and elevated risk of lung and postmenopausal breast cancer in particular [9,10].

As an indicator of biological aging, AgeAccel may also predict the onset of pubertal development. Like AgeAccel, pubertal timing has been associated with cancer risk and all-cause mortality rate among females. Earlier thelarche (age when breast growth begins) and menarche (initiation of menses), as well as increased time between thelarche and menarche (pubertal tempo) have been associated with increased breast cancer incidence [11–13]. Likewise, early menarche has been associated with a greater risk of gynecological cancers (endometrial and ovarian) [14,15]. In addition to this association with cancer incidence, a recent meta-analysis of five international studies reported a 3% decreased hazard of death from all causes for each one-year later age at menarche [16]. Despite these associations with long-term health, little is known about the determinants of pubertal tempo and timing.

For this study, we assessed the influence of AgeAccel on pubertal development in a longitudinal cohort of Chilean girls. Specifically, we analyzed the impact of AgeAccel on the time to thelarche, menarche, and pubertal tempo, as well as its association with adolescent breast density. Similar to pubertal timing, increased breast density is strongly associated with increased breast cancer risk [17–19]. It is during adolescence that the rapidly developing breast tissue is suspected to be most susceptible to carcinogens [20–22], and at which initial peak density is established [23]. However, early life predictors of breast composition are largely unknown.

Results

Our study was comprised of a subset of 94 girls from the longitudinal Growth and Obesity Cohort Study (GOCS) in Santiago, Chile with peripheral blood samples collected at Tanner 2 (B2) and Tanner 4 (B4). At the first visit of B2, approximately 59% (55/94) of the girls was categorized as healthy weight [body mass index (BMI) Z-score 5th-85th percentile], 15% (14/94) were overweight (85th-95th percentile), 25% (23/94) were obese (>95th percentile), and 2 girls were underweight (<5th percentile). Age- and sex-adjusted BMI Z-score was relatively consistent between B2 and B4 (R = 0.94, P value<0.001). The median age at thelarche was 9.25 years (95% CI: 9.02, 9.71), and the median age at menarche was 12.2 years (95% CI: 12.0, 12.4). Assuming the age at thelarche to be the midpoint between the last observed visit at Tanner 1 and the first observed visit at B2, each year delay in thelarche was associated with later age at menarche [hazard ratio (HR): 0.68; 95% CI: 0.64, 0.72]. On average, the time between the last observed Tanner 1 visit and the first observed visit at B2 was approximately 6 months (median: 0.64 years; 25th-75th percentile: 0.55 to 0.92 years; range: 0.33 to 3.02 years). Pubertal tempo was estimated as the time between thelarche and menarche; the median lag between these two developmental milestones was 2.86 years (95% CI: 2.69, 3.22). Epigenetic age (DNAmAge) was moderately correlated with chronological age (R = 0.34, P<0.001; Supplemental Figure 1). Median chronological age at the time of Tanner 2 epigenetic assessment was 9.84 years (range: 8.06 to 12.00 years), and 11.40 years at Tanner 4 (range: 9.77 to 13.10 years). Variation in DNAmAge independent of chronological age (AgeAccel) ranged from −7.70 to 8.94 years at B2 [mean: 0.03; standard deviation (SD): 3.54], and from −9.07 to 9.12 years at B4 (mean: −0.03; SD: 4.05). Adjusting for cell-type distribution, AgeAccel was moderately positively correlated between B2 and B4 (R = 0.25, P = 0.02); AgeAccel was slightly less correlated without correcting for cellular heterogeneity (R = 0.23, P = 0.03). The correlation between AgeAccel at B2 and B4 was not significantly modified by the change in fat percentage, height, or age between these two time-points (Supplemental Figure 2).

In this population, AgeAccel averaged across B2 and B4 was not significantly associated with mother's age, education, weight, age at menarche, pre-pregnancy BMI, or smoking status during pregnancy (Table 1). However, average AgeAccel had a slight but significant positive correlation with maternal height (R = 0.21; Table 1). We additionally considered the influence of childhood body size, including weight, height, and BMI, on pubertal epigenetic aging. On average, seven measurements of height, weight, and BMI sex- and age-adjusted Z-score were available for every child prior to B2. None of the associations with age-adjusted childhood body size were significantly modified by the timing of measurement (Table 2). Furthermore, none of these measurements were significantly associated with average peripubertal AgeAccel (Table 2).

Table 1. Influence of maternal and perinatal characteristics on average pubertal AgeAccel in blood (n = 94).

		Association with AgeAccel^a
Maternal Characteristics	N Missing	Unadjusted	Cell-type adjusted
Maternal Education (ref: no post-secondary)	0	1.00	0.90
Maternal Age (years)	5	−0.20	−0.08
Maternal Weight (kg)	1	−0.07	−0.09
Maternal Height (cm)	1	0.21*	0.24*
Maternal Age at menarche (years)	12	0.02	0.01
Pre-Pregnancy BMI (kg/m^2)	1	−0.13	−0.16
Maternal Smoking (ref: smokers)	4	0.02	0.37
Birth weight	1	0.14	0.08

^a Pearson correlation for continuous variables (P value for Pearson correlation), difference in average AgeAccel between categories for dichotomous variables (P value for independent T-test).

*P<0.05, **P<0.01, ***P<0.001.

Table 2. Influence of childhood body size on average cell-type adjusted pubertal AgeAccel in blood.

Childhood Growth	Median Measurements per Child (range)	Absolute Change in AgeAccel (95% CI)	Interaction with Age at Measurement^a
BMI Z-score	7 (3-11)	0.30 (−0.40, 1.00)	0.988
Weight Z-score	7 (3-11)	0.69 (−0.10, 1.48)	0.961
Height Z-score	7 (3-11)	0.74 (−0.06, 1.54)	0.550

^a Multivariable Wald test P value adding a categorical term for age {(1,3], (3,5], (5,7], (7,9], (9,15]} and interaction with childhood growth measure

*P<0.05, **P<0.01, ***P<0.001.

AgeAccel was significantly associated with the age at menarche and pubertal tempo in adjusted models (Table 3). We did not detect a significant association between AgeAccel and the timing of thelarche. A five-year increase in peripheral blood AgeAccel was associated with a significantly earlier age of menarche (HR: 1.29; 95% CI: 1.03, 1.61) adjusting for birth weight (kg), maternal pre-pregnancy BMI (kg/m²), maternal height (cm), maternal education (post-secondary Yes/No), as well as cell-type distribution. Compared to girls with epigenetic age deceleration (AgeAccel: −10 to −2 years), the median age at menarche for girls with epigenetic age acceleration (AgeAccel: 2–10 years) was nearly 5 months earlier (Table 3; Figure 1). This association was consistent after further adjustment for child height Z-score and fat percentage at B2. AgeAccel had an even stronger influence on pubertal tempo (Table 3; Figure. 2). In fully adjusted models, a five-year increase in AgeAccel was associated with decreased time between thelarche and menarche (HR: 1.48; 95% CI: 1.11, 1.98). The median time between these developmental stages was 6.6 months shorter among girls with the greatest AgeAccel (2 to 10 years) compared to the deceleration (ref: −10 to −2 years; HR: 2.05; 95% CI: 1.20, 3.51). Tanner stage at AgeAccel estimation did not significantly modify the association with either age at menarche or the time between thelarche and menarche.

Table 3. Relative change in the time to menarche and the time between thelarche and menarche associated with average pubertal AgeAccel in blood.

	HR (95% CI)
	Unadjusted	Adjusted 1^a	Adjusted 2^b
Age at Thelarche (B2)
Continuous (per 5 years)	0.99 (0.85, 1.16)	0.96 (0.80, 1.16)	0.99 (0.82, 1.20)
Categorical {ref: (−10, −2]}	1.00	1.00	1.00
(−2, 2]	1.17 (0.82, 1.69)	1.21 (0.91, 1.62)	1.20 (0.88, 1.63)
(2, 10]	0.97 (0.72, 1.31)	0.98 (0.69, 1.39)	0.92 (0.65, 1.31)
Age at Menarche
Continuous (per 5 years)	1.18 (0.96, 1.46)	1.29* (1.03, 1.61)	1.37* (1.04, 1.80)
Categorical {ref: (−10,−2]}	1.00	1.00	1.00
(−2, 2]	1.17 (0.79, 1.73)	1.42* (1.01, 1.98)	1.58 (1.00, 2.52)
(2, 10]	1.42 (0.96, 2.08)	1.61* (1.04, 2.49)	1.77* (1.03, 3.03)
Pubertal Tempo
Continuous (per 5 years)	1.26 (0.99, 1.59)	1.35* (1.02, 1.79)	1.48** (1.10, 1.99)
Categorical (ref: (−10, −2])	1.00	1.00	1.00
(−2, 2]	0.90 (0.57, 1.41)	1.11 (0.73, 1.69)	1.15 (0.75, 1.76)
(2, 10]	1.66* (1.07, 2.56)	1.71* (1.00, 2.92)	2.05** (1.20, 3.51)

^a Adjusting for birth weight, maternal pre-pregnancy BMI, maternal height, maternal education, and cell-type distribution.

^b Adjusting for variables in Model 1 and child height Z-score and fat percentage at B2

*P<0.05, **P<0.01, ***P<0.001.

Figure 1. Timing of menarche between categories of average pubertal AgeAccel. Plotting a.) the Kaplain Meier cumulative incidence curve for categories of average pubertal AgeAccel; b.) the Kaplain Meier cumulative incidence curve for categories of cell-type adjusted average pubertal AgeAccel; c.) the predicted cumulative incidence curves for categories of cell-type adjusted AgeAccel from accelerated failure time models adjusting for birth weight, maternal pre-pregnancy BMI, maternal height, and maternal education; and d.) the predicted cumulative incidence curves for categories of cell-type adjusted AgeAccel from accelerated failure time models further adjusting for child height Z-score and fat percentage at B2. Color corresponds to AgeAccel category: Red = (−10, −2]; Blue: (−2, 2], and Green = (2, 10]. Dotted lines indicate the median time to menarche for each category.

Figure 2. Timing between thelarche and menarche (pubertal tempo) across categories of average pubertal AgeAccel. Plotting a.) the Turnbull estimated cumulative incidence curve for categories of average pubertal AgeAccel; b.) the Turnbull estimated cumulative incidence curve for categories of cell-type adjusted average pubertal AgeAccel; c.) the predicted cumulative incidence curves for categories of cell-type adjusted AgeAccel from accelerated failure time models adjusting for birth weight, maternal pre-pregnancy BMI, maternal height, and maternal education; and d.) the predicted cumulative incidence curves for categories of cell-type adjusted AgeAccel from accelerated failure time models further adjusting for child height Z-score and fat percentage at B2. Color corresponds to AgeAccel category: Red = (−10, −2]; Blue: (−2, 2], and Green = (2, 10]. Dotted lines indicate the median pubertal tempo for each category.

Given this influence on the rate of pubertal development, we hypothesized AgeAccel may predict developing breast composition measured at B4. A five-year increase in average adolescent AgeAccel was associated with 5% greater geometric mean percent fibroglandular volume (%FGV) adjusting for birth weight, maternal pre-pregnancy BMI, maternal height, and fat percentage at time of breast density measurement (95% CI: 1%, 10%). This association was consistent after additional adjustment for height Z-score at B2 (Table 4). In contrast, pubertal AgeAccel was not significantly associated with total breast volume (BV) or fibroglandular volume (FGV) in multivariable models (Table 4). None of these associations were significantly modified by Tanner stage at AgeAccel estimation. Cellular heterogeneity corrected AgeAccel had a similar positive relationship with %FGV in adjusted models, but these associations did not reach significance (Table 5). Adjusting for age at breast density measurement did not influence the associations observed (results not shown). However, %FGV was associated with estimated ordinal abundance of CD8⁺CD28⁻CD45RA⁻ T cells and naive CD8⁺ T cells (Supplemental Table 1), which increase and decrease with chronological age, respectively. Breast composition was measured prior to menarche among approximately 90% of the girls (n = 82). Sensitivity analyses were conducted excluding those that attained menarche prior to breast measurements to evaluate the potential influence of regular cycling on the observed associations. Both AgeAccel and cell-type corrected AgeAccel were more strongly and significantly associated with %FGV in adjusted models restricted to breast measurements obtained prior to menarche (Supplemental Tables 2 and 3).

Table 4. Relative change in breast volume and composition at B4 associated with average pubertal AgeAccel in blood.

	Percent Change in Geometric Mean (95% CI)
	Model 1^a	Model 2^b	Model 3^c
Total Breast Volume
Continuous (per 5 years)	0.96 (0.89, 1.03)	0.95 (0.88, 1.03)	0.95 (0.88, 1.03)
Categorical {ref: (−10, −2]}	1.00	1.00	1.00
(−2, 2]	1.02 (0.91, 1.14)	1.01 (0.91, 1.12)	1.01 (0.91, 1.12)
(2, 10]	0.92 (0.81, 1.05)	0.92 (0.81, 1.05)	0.92 (0.80, 1.05)
Fibroglandular Volume
Continuous (per 5 years)	1.01 (0.92, 1.11)	1.01 (0.92, 1.11)	1.01 (0.92, 1.10)
Categorical}ref: (−10, −2]}	1.00	1.00	1.00
(−2, 2]	1.07 (0.92, 1.24)	1.06 (0.93, 1.22)	1.06 (0.93, 1.22)
(2, 10]	1.01 (0.87, 1.18)	1.01 (0.87, 1.18)	1.01 (0.87, 1.17)
Percent Fibroglandular Volume
Continuous (per 5 years)	1.05* (1.01, 1.10)	1.05* (1.01, 1.10)	1.05* (1.01, 1.10)
Categorical}ref: (−10,−2]}	1.00	1.00	1.00
(−2, 2]	1.05 (0.96, 1.15)	1.06 (0.97, 1.15)	1.05 (0.97, 1.15)
(2, 10]	1.09* (1.00, 1.19)	1.10* (1.01, 1.20)	1.10* (1.01, 1.19)

^a Adjusting for fat percentage at breast density measurement.

^b Adjusting for variables in Model 1 and birth weight, maternal pre-pregnancy BMI, maternal height, and maternal education.

^c Adjusting for variables in Model 2 and child height Z-score at B2.

*P<0.05, **P<0.01, ***P<0.001.

Table 5. Relative change in breast volume and composition at B4 associated with average pubertal AgeAccel in blood correcting for cellular heterogeneity.

	Percent Change in Geometric Mean (95% CI)
	Model 1^a	Model 2^b	Model 3^c
Total Breast Volume
Continuous (per 5 years)	0.96 (0.90, 1.03)	0.95 (0.89, 1.02)	0.95 (0.89, 1.02)
Categorical {ref: (−10, −2]}	1.00	1.00	1.00
(−2, 2]	1.04 (0.93, 1.16)	1.02 (0.91, 1.15)	1.02 (0.91, 1.15)
(2, 10]	0.96 (0.84, 1.09)	0.94 (0.82, 1.09)	0.94 (0.82, 1.09)
Fibroglandular Volume
Continuous (per 5 years)	0.99 (0.91, 1.09)	0.99 (0.90, 1.09)	0.99 (0.90, 1.09)
Categorical {ref: (−10, −2]}	1.00	1.00	1.00
(−2, 2]	1.08 (0.94, 1.24)	1.06 (0.94, 1.20)	1.06 (0.94, 1.20)
(2, 10]	1.04 (0.88, 1.23)	1.04 (0.88, 1.22)	1.03 (0.87, 1.22)
Percent Fibroglandular Volume
Continuous (per 5 years)	1.03 (0.98, 1.09)	1.04 (0.98, 1.10)	1.04 (0.98, 1.10)
Categorical {ref: (−10,−2]}	1.00	1.00	1.00
(−2, 2]	1.04 (0.95, 1.13)	1.04 (0.95, 1.13)	1.04 (0.95, 1.13)
(2, 10]	1.08 (0.98, 1.20)	1.09 (0.99, 1.21)	1.09 (0.99, 1.20)

^a Adjusting for fat percentage at breast density measurement and cell-type distribution.

^b Adjusting for variables in Model 1 and birth weight, maternal pre-pregnancy BMI, maternal height, and maternal education.

^c Adjusting for variables in Model 2 and child height Z-score at B2.

*P<0.05, **P<0.01, ***P<0.001.

We additionally considered the association between the annual change in adolescent AgeAccel and pubertal development. To estimate the change in AgeAccel, we scaled the difference in AgeAccel assessed at Tanner 2 and Tanner 4 by the number of years between the two visits. The median time between these two visits was 1.44 years (range: 0.37 to 3.67 years). The median annual change in AgeAccel was −0.06 years (range: −7.60 to 11.30 years); it was 0.13 years for cell-type corrected AgeAccel (range: −12.20 to 13.40 years). Annual change in adolescent AgeAccel was not significantly associated with mother's age, education, weight, height, age at menarche, pre-pregnancy BMI, or smoking status during pregnancy (Supplemental Table 4). Change in cell-type corrected AgeAccel was also not associated with childhood height, weight, and BMI (Supplemental Table 5). We did not detect a significant association between change in AgeAccel and pubertal timing or breast composition (Supplemental Tables 6-8).

Discussion

In this study, we report that greater AgeAccel was significantly associated with a decreased time to menarche, but not thelarche. Similarly, AgeAccel was not associated with breast or pubic hair Tanner stage among girls in the UK Avon Longitudinal Study of Parents and Children (ALSPAC) [24]. The timing of pubic hair development (pubarche) is considered a consequence of the reappearance of adrenal androgen production (adrenarche). In contrast, the timing of breast budding and menarche are thought to be primarily driven by the reactivation of the hypothalamic-pituitary-gonadal axis (gonadarche), orchestrating an increase in gonadal steroid production. Due to a shared biologic pathway, we might expect concordance between the impact of AgeAccel on breast development and menarcheal age. However, there is evidence that the timing of these two events may have distinct sensitivities to exogenous factors during different critical exposure windows. Factors modifying the timing of breast budding independent of hypothalamic-pituitary-gonadal maturation could contribute to reduced correlation between these pubertal stages. Disparate influences on the age at thelarche and menarche may explain the differing secular trends in the timing of these two events. Compared to the slight trend towards earlier menarche worldwide [25–28], the shift towards earlier thelarche has been much more pronounced [29]. Possible explanations for this disassociation between the timing of thelarche and menarche include distinct associations with childhood obesity [30] or endocrine disrupting chemical exposure [31–33]. In a large Danish cohort, initiation of breast development independent of gonadotropin levels was not explained by childhood obesity, suggesting a potential greater role of environmental exposures on the timing of thelarche [34,35].

Relative to the influence of AgeAccel on menarcheal age, increased AgeAccel was more strongly associated with decreased time between thelarche and menarche (pubertal tempo). Studies of the adult health implications of pubertal development have primarily focused on associations with age at menarche. However, the period between thelarche and menarche has been hypothesized to be a critical exposure window relevant to future breast cancer risk [36]. This time represents the primary phase of ductal growth, characterized by rapid cellular proliferation and differentiation, during which the developing mammary gland may be most susceptible to exogenous exposures. Only one study has investigated the association between pubertal tempo and breast cancer risk. In the large UK Breakthrough Generations Study, decreased time between thelarche and menarche was associated with decreased risk of breast cancer, adjusting for the age at thelarche [12]. Based on the observed association between AgeAccel and pubertal tempo, we might conjecture that increased peripubertal AgeAccel may be an indicator of decreased future cancer risk. However, this theory conflicts with our evidence of the inverse correlation between AgeAccel and menarcheal age, given the increased risk of breast cancer associated with earlier menarche [11,13]. Furthermore, increased adult AgeAccel has been associated with an increased risk of postmenopausal breast cancer [10]. These findings highlight the postulated antagonistic pleiotropy of AgeAccel. Specifically, that increased AgeAccel is a marker of appropriate, rapid early development, but detrimental during adulthood, indicating an increased risk of cancer and mortality rate [5–10].

In addition to influencing the length of this period of rapid ductal development, we hypothesized that AgeAccel may be associated with peak fibro-glandular breast volume, which is suggested to be attained at Tanner 4 [23]. We displayed a positive association between AgeAccel and adolescent %FGV in adjusted models, which was attenuated after correction for cellular heterogeneity. The association between AgeAccel corrected for cellular composition and %FGV remained significant among the subset of participants for whom adolescent breast composition was measured prior to menarche. Previous research has indicated that adolescent %FGV is influenced by regular cycling [23], suggesting a connection between the influence of AgeAccel on age of menarche and breast composition. Similar to pubertal tempo and timing, greater %FGV in adult women is a strong predictor of increased breast cancer risk [17–19]. There is recent evidence to suggest that these breast cancer risk factors may share common variation in disease pathogenesis. Specifically, later age of thelarche and menarche have been associated with greater %FGV and dense breast area in adult women; a longer interval between thelarche and menarche was only associated with increased dense area [37]. Given the observed influence of AgeAccel on pubertal tempo and menarcheal age, we would expect that an increase in AgeAccel would be associated decreased breast density. Paradoxically, we found increased AgeAccel was associated with increased percent dense volume. Many reproductive and lifestyle factors have been associated with adult breast density, including age, parity, age at first birth, breastfeeding, and alcohol intake [38,39]. These changes in the breast composition over time are suspected to capture variation in the susceptibility to carcinogenesis. For example, it is during pregnancy that the breast tissue transitions from primarily the undifferentiated lobular type 1 to the more developed type 4, reducing the proportion of vulnerable epithelial cells [40]. Due to potential changes in both AgeAccel and breast density with age, it is not possible to predict how the relationship between AgeAccel and breast composition may evolve with time. However, given both AgeAccel [10] and breast density have been associated with breast cancer incidence, these findings suggest future studies should investigate whether these risk factors provide shared or distinct information about cancer propensity.

In this cohort, childhood body size and several maternal traits were not significantly associated with adolescent AgeAccel. The association between maternal characteristics and offspring AgeAccel has been previously examined in the ALSPAC UK cohort [4]. In concordance with the UK study, we identified a positive association between maternal height and AgeAccel. Similarly, maternal age and smoking were not significantly associated with longitudinal changes in childhood and adolescent AgeAccel in the UK study. Conflictingly, this prior study detected an increase in AgeAccel associated with higher pre-pregnancy BMI, as well as an inverse association between AgeAccel and maternal weight, which was not observed in this population. The estimated influence of perinatal and childhood AgeAccel on subsequent growth patterns was also evaluated in the UK cohort [24]. Early AgeAccel predicted both height and fat mass across adolescence. Among the Chilean girls, we did not detect an association between childhood body size and pubertal AgeAccel. This suggests that the variation in AgeAccel shared across developmental stages is not the component of AgeAccel that predicts growth. Otherwise, we would expect to identify a correlation between AgeAccel and body size, regardless of directionality. In general, adjusting for height and fat percentage at the time of AgeAccel measurement in this cohort strengthened the associations between AgeAccel and pubertal tempo and timing. We hypothesize that early AgeAccel may have a disparate influence on pubertal development relative to the impact of adolescent AgeAccel. These results highlight the necessity of examining the associations between AgeAccel and development across time given the susceptibility of AgeAccel to a variety of exposures, which may change the nature of the association.

Our study has several strengths, including the utilization of a well-characterized longitudinal cohort with frequent assessment of childhood and adolescent growth. Estimation of AgeAccel at two time points enabled a nuanced evaluation of the impact of AgeAccel across puberty on developmental outcomes. A limitation of the timing of these measurements is that we did not assess AgeAccel prior to thelarche. Any association between AgeAccel and thelarche may have been masked by modification of pubertal AgeAccel by the timing of thelarche. However, the lack of influence on the age at thelarche in this cohort echoes a previously reported null association in ALSPAC [24]. Our study may also suffer from potential residual confounding by childhood stress. Childhood adversity has been associated with both pubertal timing and epigenetic aging [41,42]. Although we adjusted for an indicator of socioeconomic status, there may be bias due to other external stressors or trauma not accounted for in this study. Another limitation is that these findings may not be generalizable to other race/ethnicities. The timing of pubertal development has been observed to vary across populations, adjusting for BMI, as well as social and economic indicators [26,43,44]. This suggests that there may be either differing biological susceptibilities to exogenous determinants of pubertal timing, differing frequencies of genetic predictors of timing, or differing distributions of unmeasured confounders. Regardless, future studies in distinct populations will be necessary to confirm the observed associations among Latina girls are consistent across populations. Given Latina girls are suspected to develop earlier [26,43], this study provides unique insight into unstudied susceptible sub-population.

Overall, our study indicates faster adolescent AgeAccel is associated with faster pubertal development. While increased AgeAccel did not significantly affect the onset of breast development, it was associated with decreased time between thelarche and menarche and, correspondingly, an earlier age at menarche. Based on prior research, we would expect a shorter interval from thelarche to menarche and earlier menarche to be associated with decreased fibro-glandular volume [12]. However, increased AgeAccel was also associated with increased breast density. Given the association between pubertal development and adult health, these findings should spur additional research into the future health implications of this accelerated pubertal epigenetic clock ticking.

Materials and methods

Study population

Our study population was the longitudinal Growth and Obesity Cohort Study (GOCS) in Santiago, Chile. Initiated in 2006, children ages 2.6 to 4.0 years were recruited from public nursery schools of six counties in Santiago who met the following inclusion criteria: 1) singletons born in 2002-2003 and with birthweight between 2,500 and 4,500 g and 2) absence of physical (e.g., skin burn), medical (e.g. brain tumor), or endocrine diseases (e.g., hyperthyroidism, hyperprolactinemia) that could alter the growth and/or onset of puberty. Of 1,498 eligible participants, the mothers of 1,195 children (∼80%) accepted to participate in the study. There were no significant differences in age, gender, and birth anthropometry between the final participants and those not enrolled. The GOCS children are representative of the low to middle-income families served by the public nursery schools, which provide free education and food 5 days per week during 11 months of the year. Of the total cohort, 602 are girls with longitudinal data collection. The current study included a subset of 94 randomly selected from the 166 girls with peripheral blood specimens collected at both Tanner 2 (B2) and Tanner 4 (B4) (of 239 with B2 DNA, and 334 with B4 DNA). Informed consent was obtained from all parents or guardians of children before the start of data collection. The study protocol was approved by the Ethics Committee of the Institute of Nutrition and Food Technology, University of Chile.

Breast composition measurements

Breast density was measured at B4, when total fibroglandular volume is assumed to reached its peak [45]. Breast development was assessed by visual inspection using Tanner's rating scale approximately every 6 months [46]. After confirming the absence of pregnancy by a urine test, dual X-ray absorptiometry (DXA) was used to measure the volume of dense tissue (absolute fibroglandular volume) of the breast. The UCSF DXA breast scanning protocol was developed by Shepherd et al. in the Department of Radiology and Biomedical Imaging, University of California, San Francisco (version 5) [47]. Each breast was scanned with the use of the GE iDXA system software (version 13.6, GE Healthcare, Madison, WI, USA). A quality control phantom containing reference breast density materials was scanned throughout the study to assure a stable calibration. This approach has been shown to have high validity and precision for measuring breast density in girls at different Tanner stages [45]. The radiation dose from DXA is exceedingly small; DXA is commonly used to measure whole body percent fat and body mass in pediatric studies [48]. Total projected breast area was manually delineated on each image and breast fibroglandular volume (FGV; cm³) and total volume (BV; cm³) were estimated using a two-compartment model of adipose and fibro-glandular tissue with a software developed by Shepherd et al. at the University of California [49]. Percent FGV (%FGV) was defined as the proportion FGV relative to BV times 100. We averaged the values of the left and right breast for all analyses.

Tanner stage

Starting in 2009, breast development was assessed by two trained dietitians (kappa with pediatric endocrinologist = 0.85) by visual inspection using Tanner's rating scale approximately every 6 months [46]. Palpation was additionally used to distinguish breast Tanner 1 from Tanner 2. For statistical analysis, the timing of breast development (thelarche) was assumed to be between the age at last Tanner 1 (B1) visit and the age at the first Tanner 2 (B2) visit. None of the girls included in this study population reached B2 prior to entry into the cohort, and all girls reached B2 during the course of the study.

Anthropometry

Maternal pre-pregnancy body mass index (BMI; kg/m²) was self-reported at study entry. Birth weight (kg) of the daughters, as well as annual height and weight measurements prior to 2006, were abstracted from medical charts. After study initiation, body fat percentage was estimated at each visit using Tanita-BC-418 MA bioelectrical impedance measurements (Tanita-Corporation, Tokyo, Japan), according to the manufacture's guidelines and at measurement frequency of 50 kHz (accuracy 0.1 kg) [50]. Sex- and age-adjusted height, weight, and BMI Z-scores at each visit were calculated based on the World Health Organization (WHO) growth charts.

Age at menarche

Prior to the onset of B4, girls were asked to report their first menstrual bleeding at each 6-month visit. After achieving B4, girls were contacted by study dieticians every 3 months to survey whether the girl had reached menarche. During this phone interview, a questionnaire was used to differentiate menarche from other potential causes of vaginal bleeding, such as vaginal infection, urinary infection or trauma. Longitudinal follow-up of participants enabled the confirmation of menarche onset. Models for the time to menarche include incident cases, as well as right censored individuals (individuals that did not reach menarche during the course of the study).

Epigenetic age acceleration

Peripheral blood buffy coat DNA methylation (DNAm) at B2 and B4 was measured in bisulfite-treated DNA using the Illumina Infinium MethylationEPIC array at the USC Epigenome Center (Los Angeles, CA). Data pre-processing was performed utilizing a modified beta-mixture quantile normalization method that was originally developed by Teschendorff et. al [1,51]. This normalization method both reduces the influence of technical artifacts across batches and standardizes the data to the training data of the epigenetic clock. Epigenetic age was estimated using the Horvath method, which predicts age based on the DNAm of 353 CpG loci on the Illumina microarray using previously reported coefficient values [1]. Estimated proportions of CD8⁺ T cells, CD4⁺ T cells, natural killer cells, B cells, monocytes, and granulocytes data were calculated using the Houseman approach [52], and estimated abundance measures of plasma B cells, CD8⁺CD28⁻CD45RA⁻ T cells, naive CD8⁺ T cells, and naive CD4⁺ T cells were calculated using the advanced analysis option of the epigenetic clock software [1]. Epigenetic age was regressed on chronological age, and residuals from this model, termed epigenetic age acceleration (AgeAccel), were used for subsequent analyses. By design, these estimates are independent of chronological age, with positive values indicating epigenetic age acceleration relative to actual age, and negative values indicating deceleration. Epigenetic age was also regressed on age as well as estimated cell-type composition (age + naive CD8⁺ T cells + CD8⁺CD28⁻CD45RA⁻ T cells + plasma B cells + CD4⁺ T + natural killer + monocytes + granulocytes) to provide a measure of AgeAccel independent of cellular heterogeneity.

Statistical methods

Accelerated failure time (AFT) models were used to assess the influence of AgeAccel on time to thelarche and menarche as well as pubertal tempo, assuming a Weibull distribution. A cluster statement was used to account for within-subject correlation between B2 and B4 AgeAccel estimates. Accordingly, inference was based on robust standard errors estimated using the Huber sandwich estimator. Use of AFT models allowed us to account for interval censoring of thelarche onset between the age at last B1 visit and the age at first B2 (years). For the time to menarche analysis, survival time for incident cases was the age at menarche (years), estimated based on time between the self-reported date of first menses and date of birth. Survival time for right censored individuals was the age at last clinic visit, based on the time between date of last visit and date of birth. Pubertal tempo was defined as the time between thelarche and menarche, accounting for interval censoring of thelarche onset. For incident cases of menarche, intervals were defined as the time (years) from the date of last B1 visit to menarche, and the time (years) from first B2 visit to menarche. If menarche was right censored, we estimated the time to the last clinic visit, and including an indicator for censored outcome. Time-varying associations between AgeAccel and the age at thelarche and menarche were investigated by adding Tanner stage at blood sample, as well as an interaction between AgeAccel and Tanner stage, to our models. We similarly investigated effect modification by Tanner stage at AgeAccel measurement on pubertal tempo. If Tanner stage at AgeAccel estimation did not significantly modify the association with pubertal endpoint (Wald test; P<0.05), we reported associations with average AgeAccel across time points. Generalized estimating equations (GEEs) were used to jointly estimate the influence of AgeAccel at B2 and B4 on breast density measurements, including: total breast volume (cm³), absolute FGV (cm³), and percent FGV (%). This non-conventional application of GEE models allows for the simultaneous modeling of the associations with B2 and B4 AgeAccel, with a working independence assumption and model based standard errors [53,54]. Breast density measurements were log-transformed prior to analysis to improve normality. Time-varying associations between AgeAccel and breast density were investigated by evaluating the significance of the interaction between AgeAccel and Tanner stage on breast density. Reported associations were stratified by Tanner stage if significant effect modification was detected (Wald test; P<0.05). Associations between breast density measurements and AgeAccel were reported as the relative change in geometric mean and 95% confidence interval (CI) by exponentiating the associations with log transformed breast density measurement. When modeled continuously, AgeAccel was scaled so that a one-unit change represented a 5-year shift in AgeAccel. Accordingly, we report the relative change in hazard associated with a 5-year increase in AgeAccel. We additionally modeled categories of AgeAccel to identify potential non-monotonic relationships and for ease of interpretation. All models were adjusted for estimated cell-type distribution, maternal pre-pregnancy BMI, maternal height, maternal education, as an indicator of socioeconomic status, and birth weight, which have been previously associated with AgeAccel [4,55]. Maternal BMI and birth weight represent possible confounders, given the reported impact of both maternal BMI and birth weight on age at menarche [56,57], as well the association between birth weight and adult breast density [58]. Socioeconomic status has similarly been associated with both pubertal timing and breast density [43,59]. There is suggestion that maternal height may also be a predictor of breast cancer risk, but the literature is limited [60]. We considered models further adjusted for child's height Z-score and fat percentage at B2, which have been previously associated with perinatal or childhood AgeAccel [24]. Breast density models were additionally adjusted for fat percentage at measurement to account for the influence of adiposity on breast volume. Missing adjustment variables were mean imputed. GEEs were also used to assess the influence of childhood growth, including BMI, height, and weight Z-score, on average AgeAccel across puberty. We restricted to body size measurements that occurred before the first visit of B2. Effect modification by categories of age at body size measurement {(1,3 years], (3,5 years], (5,7 years], (7,9 years], (9,15 years)} was evaluated via multivariate Wald test (P<0.05). Similar models were run to evaluate the association between annual change in adolescent AgeAccel and pubertal development. Annual change in adolescent AgeAccel was estimated by the difference in AgeAccel estimated at Tanner 2 and Tanner 4, scaled by the number of years between the two visits. The only differences in the analytic approach were the use of standard AFT models to estimate the association between change in AgeAccel and pubertal timing, and the use of generalized linear models to evaluate the relation with breast composition. All statistical analysis was performed in R Version 3.4.1 and figures were generated using ggplot2 [61].

Disclosure of potential conflicts of interest

No potential conflict of interest was reported by the authors. NIH grants supporting this work are listed under funding details.

Additional information

Funding

This work was supported by Public Health Service grant R01CA158313 from the National Cancer Institute, National Institutes of Health, US Department of Health and Human Services (to KBM), and by the Breast Cancer and the Environment Research Program (BCERP) award grant U01ES026130 from the National Institute of Environmental Health Sciences and the National Cancer Institute, National Institutes of Health, Department of Health and Human Services (to KBM).