Determining the timing of pubertal onset via a multicohort analysis of growth

Objective Growth-based determination of pubertal onset timing would be cheap and practical. We aimed to determine this timing based on pubertal growth markers. Secondary aims were to estimate the differences in growth between cohorts and identify the role of overweight in onset timing. Design This multicohort study includes data from three Finnish cohorts—the Type 1 Diabetes Prediction and Prevention (DIPP, N = 2,825) Study, the Special Turku Coronary Risk Factor Intervention Project (STRIP, N = 711), and the Boy cohort (N = 66). Children were monitored for growth and Tanner staging (except in DIPP). Methods The growth data were analyzed using a Super-Imposition by Translation And Rotation growth curve model, and pubertal onset analyses were run using a time-to-pubertal onset model. Results The time-to-pubertal onset model used age at peak height velocity (aPHV), peak height velocity (PHV), and overweight status as covariates, with interaction between aPHV and overweight status for girls, and succeeded in determining the onset timing. Cross-validation showed a good agreement (71.0% for girls, 77.0% for boys) between the observed and predicted onset timings. Children in STRIP were taller overall (girls: 1.7 [95% CI: 0.9, 2.5] cm, boys: 1.0 [0.3, 2.2] cm) and had higher PHV values (girls: 0.13 [0.02, 0.25] cm/year, boys: 0.35 [0.21, 0.49] cm/year) than those in DIPP. Boys in the Boy cohort were taller (2.3 [0.3, 4.2] cm) compared with DIPP. Overweight girls showed pubertal onset at 1.0 [0.7, 1.4] year earlier compared with other girls. In boys, there was no such difference. Conclusions The novel modeling approach provides an opportunity to evaluate the Tanner breast/genital stage–based pubertal onset timing in cohort studies including longitudinal data on growth but lacking pubertal follow-up.


Introduction
Puberty is a period of rapid growth involving changes in hormonal activity. Most recent studies on pubertal timing have reported a shift of puberty toward earlier ages during the last decades, especially in girls [1,2]. The commonly used marker of pubertal onset in girls is the development of breasts [3]. The emergence of breast buds defines Tanner breast stage 2 (Tanner B2) [4]. In boys, the commonly used marker of pubertal onset is testicular enlargement [3]. An enlargement to the volume of > 3 mL defines Tanner genital stage 2 (Tanner G2) [4], and this volume is consistent with a testis length � 25 mm [5]. However, puberty-induced body changes are a sensitive issue for children and adolescents, and thus, the physical measurement of pubertal onset may be challenging; at the same time, indirect measurement-for example, by hormone levels-is expensive.
One well-known feature of puberty is the adolescent growth spurt. The timing and intensity of the peak of the growth spurt are measured by the peak height velocity (PHV, cm/year) and age at PHV (aPHV, years). Estimation of the PHV and aPHV requires the collection of repeated height measurements and an appropriate analysis methodology. Multiple mathematical models have been proposed to describe adolescent growth [6][7][8][9][10]; the Super-Imposition by Translation And Rotation (SITAR) growth curve model is one of the most recent [11]. The SITAR model has been shown to efficiently summarize growth in adolescence [8] and determine the PHV and aPHV. In addition, it has been suggested that joint modeling of several cohorts with the SITAR model can improve the fit, and that the model can straightforwardly produce cohort-specific growth parameters [12].
In some studies, the start of the growth spurt has been used as the first sign of pubertal maturation [9,10]. In addition, different markers for different stages of pubertal development have been used in the existing literature [13,14], and growth in association with pubertal development has been widely studied [15][16][17][18][19]. However, the literature does not present the relation of the precisely estimated aPHV and PHV with the age at pubertal onset as a continuous variable, with consideration of interval censoring. When researchers are interested in determining the age at pubertal onset as defined by the Tanner breast/genital stages, but the stages have not been measured, determination based on pubertal growth markers becomes an attractive possibility. Since height data are often recorded as a part of routine health surveillance and collected in many clinical follow-up studies, it is generally possible to determine the aPHV and PHV. aPHV is a well-known milestone of puberty, and it correlates with the age at pubertal onset, as does PHV [20,21]. Previous studies have indicated differences in anthropometrics between regions in Finland. For example, adults living in Southwest Finland have been found to be taller than adults living in Northern Ostrobothnia [22]. The present study can shed light on whether this difference can be seen during adolescence.

PLOS ONE
The association between an increased body mass index (BMI) in childhood and a notable shift to earlier puberty has been widely reported for girls [23,24]. For boys, the evidence of such an association is contradictory [23], with some studies suggesting the association to earlier [24,25] and some studies to later pubertal onset [26]. There is also evidence for different directions of association depending on the magnitude of increase in BMI [27].
The primary aim of the present multicohort study was to estimate the pubertal growth markers, aPHV and PHV, and to determine the timing of pubertal onset based on them. The secondary aims were to investigate whether there are differences in growth during adolescence between three recent Finnish child and adolescent studies and to assess the role of being overweight or obese in the timing of pubertal onset.

Study population and design
The data of the current multicohort study are derived from three Finnish cohort studies-the Type 1 Diabetes Prediction and Prevention (DIPP) Study [28], the Special Turku Coronary Risk Factor Intervention Project (STRIP) [29,30], and the Boy cohort [31]. Multi-cohort study data were collected in two different regions and during different time periods (Table 1).
The DIPP Study is a Finnish prospective population-based birth cohort study. Newborn infants born in Oulu, Tampere and Turku University Hospitals were screened for HLA-DQB1-conferred susceptibility to type 1 diabetes using cord blood samples. Infants carrying increased genetic susceptibility were recruited for regular follow-up. The current study includes 2825 children born in Oulu University Hospital.
The STRIP is a Finnish prospective randomized intervention trial aimed at preventing coronary heart disease and atherosclerosis by a lifestyle-based intervention. One thousand sixtytwo 7-month-old infants born in 1989-1991 were recruited at well-baby clinics in Turku. They were randomized into a life-style based intervention group (n = 540) or a control group (n = 522).
The Boy cohort is a Finnish prospective population-based birth cohort study on the prevalence of congenital cryptorchidism. It included 1494 antenatally recruited infants born 1997-1999 in Turku University Hospital. In addition, 184 boys born in 1997-2002 were included in a case-control study on risk factors of congenital cryptorchidism. Of all the boys, 119 participated in a pubertal follow-up with 52 of them having (cases) and 67 not having (controls) a history of congenital cryptorchidism. The data of the controls in the pubertal follow-up study were used in the current analysis.
Children in all cohorts were monitored for growth; children in the STRIP and Boy cohort were also monitored for pubertal markers (Tanner staging) from the age of 8 years. Tanner staging included the following assessments: • Breast development in girls, evaluated by palpation and measurement of breast tissue with a ruler (Tanner B).

PLOS ONE
• External genitalia in boys, determined by measuring the length of the testes with a ruler (Tanner G), as described in Sadov et al. [31] • Pubic hair for both sexes (Tanner pubic hair stage).
The DIPP Study was included in the analyses based on the further interest to determine the timing of pubertal onset based on growth data in the DIPP children and to improve the fit of the growth curves [12].

Ethics
The DIPP Study (Oulu) was approved by the Ethical Committee of the Faculty of Medicine, University of Oulu, and the Ethical Committee of the city of Oulu. Parents gave their written informed consent first for the genetic testing of the newborn infant and then again for participation in the follow-up.
The STRIP was approved by the Joint Committee on Ethics of the Turku University and the Turku University Central Hospital. Parents of the children gave their written informed consent.
The original birth cohort study, the case-control study and the pubertal follow-up study of the Boy cohort were approved by the Joint Committee on Ethics of the Turku University and the Turku University Central Hospital. Parents of the children gave their written informed consent and the participants of the pubertal follow-up gave their assent for the study.
All the cohorts adhere to the Declaration of Helsinki. Data were pseudonymized before the corresponding author accessed them for the analyses. The written informed consents given within the cohort studies enable the usage of collected data in research, including the current multi-cohort study.

Defining the age at pubertal onset
For girls, pubertal onset was determined based on Tanner breast stages. The onset of puberty was considered to occur between the last observed Tanner B1 and the first observed Tanner B2. In a time-to-pubertal onset model, the transition between the stages was understood as an interval censored between those points of time. If a girl had missing Tanner B2 measurements, but Tanner B3, B4, and/or B5 measurements were available, the time interval between the last Tanner B1 and the next available Tanner B measurement was used. In the original data provided for boys, Tanner G2 was considered to occur when testicular length was � 20 mm. However, instead of using this old threshold, based on the current recommendation of a testis length of 25 mm as the marker of pubertal onset, 5 ages at pubertal onset were recalculated using testis length measurements. The pubertal onset was considered to occur when the length of the larger testis was � 25 mm at two successive visits [31]. Then, the age at pubertal onset was considered as the interval censored between the last < 25 mm and the first (of the two successive visits) � 25 mm measurement in the time-to-pubertal onset model [32].
If the last measurement was taken before pubertal onset, the age at pubertal onset was handled as right censored: The child had dropped out before pubertal onset, and we only observed that the onset had occurred after that age [32]. If there were no measurements before pubertal onset had taken place, the age at pubertal onset was handled as left censored: We only observed that the age at pubertal onset was less than the age at first measurement [32].

Statistical analyses
The analyses were performed in two stages and conducted separately for boys and girls with R 3.5.3 [33] using the "sitar" function from the "sitar" package [34], the "survreg" function from the "survival" package [35], and the "glmnet" function from the "glmnet" package [36].
Data from all the cohort studies were used in the first stage for SITAR growth curve analyses to investigate the differences in growth between the cohort studies. All growth data from individuals between the ages of 7 and 19 were included. In the second stage, for the time-topubertal onset analyses, only data from the STRIP and Boy cohort were used because they included measurements on pubertal development. In this analysis, children with at least four height measurements were included. This decision was based on the visual inspection of the distributions of the aPHV values by different numbers of height measurements (S1 Fig). The variation in aPHV values was minimal with smaller number of height measurements indicating that the aPHV estimation is based more on the average growth curve than on an individual curve.
First stage: Super-Imposition by Translation And Rotation (SITAR) growth curve model. In the first stage, growth curves were estimated with the SITAR growth curve model, which has been shown to efficiently summarize growth in adolescence using the three dimensions of size, timing, and intensity [8,11]. Size corresponds to the mean height, timing to the aPHV, and intensity to the PHV. The model estimates a mean growth curve through fixedeffect parameters for the dimensions. It also estimates the random-effect parameters, representing the shift from the mean curve for each individual. The random effect for size represents a shift up or down of the height curve of a particular child relative to overall mean height curve (measured in centimeters). The random effect for timing represents a shift earlier or later of the aPHV of a particular child relative to mean APHV (measured in years). The random effect for intensity represents a change of the PHV of a particular child relative to mean PHV (can be viewed as a percentage difference when multiplied by 100). Individual growth curves are estimated through the fixed-effect and random-effect parameters.
The mean aPHV and PHV were derived by differentiating the mean height curve estimated by the SITAR model and by identifying the age at which the derivative curve reached its maximum. The PHV was the value of the derivative at that age. Individual aPHVs and PHVs were derived from the individual height curves in the same way. The differences in growth between cohorts were investigated based on the cohort-specific fixed-effect parameters of the SITAR model. Full technical details of the SITAR modeling are provided in S1 File.
Overweight classification. The overweight status of each child was used as a predictor in the time-to-pubertal onset models to investigate the role of being overweight in the timing of pubertal onset. The overweight status was determined based on the BMI when the children's growth spurt started. The start of the growth spurt was defined as the lowest point in the derivative of the individual growth curve (velocity curve). The corresponding BMI for the child was obtained from the individual BMI-age trajectory fitted by a polynomial spline mixed-effects model. The BMI-age mixed-effects model included the random intercept and random effects for all the polynomial regression coefficients. Children were classified into two classes based on the obtained BMI value as follows: overweight (ISO-BMI > 25, including overweight and obese) and other (including normal weight and underweight) children, based on the Finnish ISO-BMI classification for children [37].
Second stage: Time-to-pubertal onset model. In the second stage, two time-to-pubertal onset models-an extended model and a simple model-were fitted to determine the timing of pubertal onset based on growth markers and to study the role of overweight in the timing. In both models, the outcome was the timing of pubertal onset. The simple model used only aPHV and PHV as predictors. The purpose of the model was to study whether aPHV and PHV alone could determine the timing of pubertal onset. The extended model included all the important variables, which were selected using a lasso regression technique with a time-topubertal onset model. The considered variables were the aPHV, PHV, and size parameter (a i ) obtained from the SITAR model; the overweight status of the child; and all the second-order interactions between the variables. Full technical details of the time-to-pubertal onset models are provided in S1 File.
Determination of the predictive ability of the model. The predictive ability of the time-topubertal onset model was determined using 10-fold cross-validation to obtain the predicted pubertal onsets. Using cross-validation, we could investigate the predictive ability in external data more reliably because the prediction was always done for children who were not included in the estimation of the model. Ten-fold cross-validation was chosen because it has been indicated to be one of the methods with the smallest bias [38].
The estimated individual means at given values of the predictors were used as predicted pubertal onsets. Prediction success was evaluated by comparing the observed pubertal onset age interval to a year length interval centered on the predicted age at pubertal onset. A year length interval was chosen to agree with the prescheduled visits for most of the children. If those intervals overlapped, it was considered to represent an agreement. Only children with an observed interval length of less than 1.5 years were included in the evaluation of success. To visualize the success, prediction intervals were obtained using the appropriate quantiles of the normal distribution centered on the individual prediction and with an estimate of the variance (ŝ 2 ).

Determination of the timing of pubertal onset based on growth markers
The mean age at pubertal onset was 10.5 years (SD: 1.25) for STRIP girls and 11.7 (SD: 0.90) and 11.3 (SD: 0.98) years for STRIP and Boy cohort boys, respectively, as estimated from separate cohort-specific time-to-pubertal onset models. Of the 340 girls and 432 boys with pubertal data, 306 girls (90.0%) and 392 boys (90.7%) fulfilled the criterion of four height measurements and were included in the time-to-pubertal onset analyses (Stage 2). Of those children, for 7 (2.3%) girls and 40 (10.2%) boys, pubertal onset occurred after the last follow-up visit (right censored), and for 4 (1.3%) girls, it occurred before the first follow-up visit (left censored).
aPHV had an apparent direct association and PHV an apparent inverse association with age at pubertal onset in both sexes (S2 Fig). Only the size parameter was not associated with the timing of pubertal onset (S2 Fig). For both sexes, the important variables chosen for the extended time-to-pubertal onset model were aPHV, PHV, and overweight status. For girls, the interaction between aPHV and overweight status was added to the model because the association of aPHV to the timing of pubertal onset appeared to differ depending on the overweight status (S3 Fig). aPHV clearly dominated the prediction over PHV for both sexes. All the timeto-pubertal onset model parameters and prediction results are presented in S1 Table. The predictive ability of the extended model was only slightly better than the predictive ability of the simple model for both sexes (S1 Table).
The extended model predictions tended to fall in the center of the data, suggesting good predictive ability for the test data (Fig 1). Most of the observed ages at pubertal onset overlapped with the prediction intervals (Fig 1), indicating that the predicted onset falls close to the observed onset. The agreement for girls was 71.0%, and for boys, it was 77.0%.

The role of being overweight in the timing of pubertal onset
Of the girls, 51 (15.0%) were overweight. Girls with overweight reached puberty earlier than other girls did (−1.0 [95% CI: −1.4, −0.7] years). However, the improvement in overall predictive ability was modest: The agreement increased from 68.1% to 71.0% while adding

PLOS ONE
The pubertal growth markers and the timing of pubertal onset Table 2  overweight status and its interaction with aPHV to the simple model (S1 Table). The modest improvement was likely due to lower aPHVs among overweight girls compared with other girls.
Of the boys, 77 (17.8%) were overweight. The timing of pubertal onset did not differ between boys with and without overweight (−0.0 ([95% CI: −0.3, 0.2] years). Overweight status was still found to be a necessary variable in the boys' extended model. Nevertheless, the estimated regression coefficient of the overweight status (0.2 [95% CI: 0.0, 0.4] years) remained small in the extended model, and hence, it did not improve the overall prediction results dramatically: The agreement increased from 76.0% to 77.0% (S1 Table). Specifically, overweight appeared to be associated with later pubertal onset among boys with larger aPHV values.

Discussion
We determined the timing of pubertal onset based on growth in terms of individual-specific time intervals representing the period of pubertal onset. In the results, cross-validation showed a good agreement between the observed and predicted timings. The children in Southwest Finland (STRIP and Boy cohort) were taller than the children from Northern Ostrobothnia (DIPP). In addition, STRIP children had higher PHV values compared with DIPP children. Overweight girls had an earlier pubertal onset than normal-weight girls did. In boys, there was no such difference.
Khairullah et al. [16] investigated the relationships between three different height markers, including aPHV, and the pubertal stage at one measurement point, at the age of 13.1-13.8 years. In contrast to the present study, pubertal stage was self-assessed and the association of aPHV with the timing of onset was not considered. Gasser et al. [17] studied the chronology of a variety of pubertal developmental milestones and the associations between them. Associations were studied by determining the proportions of children who reached a specific milestone at some level of the other milestone. Cole et al. [14] determined the most intense stage of pubertal development based on separate applications of different markers, including height and Tanner G/B. However, these researchers did not study the associations between the markers. Zhu et al. [19] determined the pubertal status (prepubertal, pubertal, and post-pubertal) of individuals at single visits based on the height velocity between visits, growth chart review, and clinical adjudication using predefined threshold velocities [18], and they compared their assessments to the clinical Tanner staging visit by visit. However, they did not use modeling to determine the growth-based pubertal status, and they did not determine the age at pubertal onset or aPHV at all. McCormack et al. [20] presented the association between aPHV, PHV, and age at pubertal onset determined using similar definitions as in the present study. They used the age at the first observation of pubertal onset as a timing of onset and classified ages into three groups, whereas our model used continuous values and considered the interval censoring of the outcome, and provided point and interval estimates of the age at pubertal onset. Their results agree with ours: aPHV was directly and PHV was inversely related to pubertal onset age. Our study provides an alternative model to theirs based on Finnish data.
Since the previous literature indicates an association between increased BMI and age at pubertal onset [23,24,26], we considered it important to include the overweight status of the children in the time-to-pubertal onset model. We found that boys with overweight had a similar age at onset of puberty compared with other boys, but they exhibited an earlier aPHV. This finding is similar to that in a recent study from Denmark [25], where testicular volume � 4 mL occurred significantly earlier in obese compared with normal-weight boys but Tanner genital stage-based ages at pubertal onset did not differ between those boys. Together, these findings indicate that one possible explanation for the previous contradictory results in boys could be the variation in the chosen pubertal marker. Compared with other girls, those who were overweight exhibited earlier values in terms of both onset of puberty and aPHV. The observations of earlier aPHV values among girls and boys with overweight are consistent with previous literature [39].
The extended time-to-pubertal onset model we used gave only slightly better prediction results than the simple model that included only aPHV and PHV. Thus, using these two growth markers as the only predictors could be appropriate. We think that these parameters need to be determined sufficiently well, and during the study, we found that the SITAR model was well suited for this purpose. We determined the overweight status at the age of the minimum height velocity (based on the SITAR curves) before the growth spurt. In our data, a more straightforward classification by BMI at 2 years before the aPHV would not have resulted in any major changes because the overweight status at these two alternative times agreed nearly perfectly.
We observed a similar difference in height in adolescents between the Oulu (DIPP) and Turku (STRIP, Boy cohort) regions, as was previously presented with adult data [22]. The finding of the higher PHV values in STRIP children compared with DIPP children could be related to the same phenomenon. Some environmental and genetic factors could explain the difference. In addition, type 1 diabetes is known to adversely affect linear growth [40]. Thus, the genetic susceptibility to type 1 diabetes of DIPP participants may be one factor that explains the differences to the DIPP cohort. However, the risk of developing type 1 diabetes with genetic susceptibility is only 1.7%-8.0% depending on the risk genotype [41], and it is unclear whether genetic susceptibility alone affects growth.
The major strengths of the present study were the large, unique dataset, representing a combination of three carefully followed cohorts. Growth was regularly measured during childhood and adolescence in all the cohorts, and pubertal markers were assessed in two of the cohorts. The SITAR model enabled the straightforward and effective production of cohort-specific and individual-specific growth parameters, and the time-to-pubertal onset model enabled the inclusion of censored pubertal onsets in the model.
The lack of pubertal measurements in the DIPP Study was one limitation of the present research. We had to leave the DIPP Study out of the pubertal analyses, which substantially reduced the amount of data available for the Stage 2 models. Although DIPP children can be a selective sample, the data included children from two other studies, extending the applicability of the results. Cross-validation was used to improve the generalizability of the predictions to external samples.
The estimation of age at pubertal onset is best viewed as an interval estimation: Such estimation provides a prediction interval that can be seen to indicate the range of the most probable ages of pubertal onset. Since real measurements also only provide interval censored data, the interpretation shares similarities with the real measurements (e.g., Tanner stages) taken at fixed times. The determination of pubertal onset timing based on pubertal growth markers is promising method and with replication in other cohorts it could facilitate the determination of the timing in adolescent cohort studies more generally.