Compared to females, males reach puberty later and develop diverse primary and secondary sexual characteristics (Tanner, Reference Tanner1981). To date, far fewer studies of pubertal timing have focused on males than females, probably due to the difficulty of assessing males’ pubertal development. As menarche is a major event in females’ lives, the timing of menarche is relatively accurately recalled by many females. However, there is no pubertal marker as salient as females’ menarche in males. The most reliable and valid marker of pubertal initiation in males is known to be an increase in testicular volume above 3 mL (Abreu & Kaiser, Reference Abreu and Kaiser2017). However, the assessment of testicular volume requires the presence of a physician or trained health professional, which makes it a suboptimal measure for sensitive populations or large-scale studies (Coleman & Coleman, Reference Coleman and Coleman2002). Recently, it has been recommended that age at voice change (AVC) be used as a proxy marker of male pubertal timing, because AVC is an easy-to-assess and noninvasive measure and because voice change is a characteristic event in males (Hodges-Simeon et al., Reference Hodges-Simeon, Gurven, Cardenas and Gaulin2013). As a late event in puberty, voice change refers to the abrupt decrease in the fundamental voice frequency due to the increased length of vocal cords following the growth spurt of the larynx in response to androgen exposure (Harries et al., Reference Harries, Hawkins, Hacking and Hughes1998). It typically begins during the transition from Tanner pubertal stage 3 to 4, approximately 2 years after testicular enlargement (Harries et al., Reference Harries, Hawkins, Hacking and Hughes1998), and is strongly correlated with other male pubertal milestones, such as testicular growth, gonadarche, sweat odor, axillary hair growth and peak height velocity (Busch et al., Reference Busch, Hollis, Day, Sørensen, Aksglaede, Perry, Ong, Juul and Hagen2019; Yousefi et al., Reference Yousefi, Karmaus, Zhang, Roberts, Matthews, Clayton and Arshad2013).
It is very well documented that the average age of puberty onset in females has declined in most developed countries during the past decades (Euling et al., Reference Euling, Herman-Giddens, Lee, Selevan, Juul, Sørensen, Dunkel, Himes, Teilmann and Swan2008), although there is some indication that the decline has leveled off in recent years (Cole, Reference Cole2000; Papadimitriou et al., Reference Papadimitriou, Fytanidis, Douros, Bakoula, Nicolaidou and Fretzayas2008). To determine whether the same secular trend occurs in male pubertal timing, a panel of researchers reviewed American data on male puberty from 1940 to 1994 and suggested that evidence was insufficient to confirm a downward trend in male pubertal onset (Euling et al., Reference Euling, Herman-Giddens, Lee, Selevan, Juul, Sørensen, Dunkel, Himes, Teilmann and Swan2008). More recently, however, the Copenhagen Puberty Study, a population-based cohort study of healthy Danish children and adolescents, showed a decline of 3 months in age at onset of puberty in males in recent years (1991–2008; Sorensen et al., Reference Sorensen, Aksglaede, Petersen and Juul2010). In another Danish study, Juul et al. (Reference Juul, Magnusdottir, Scheike, Prytz and Skakkebaek2007) have found downward trends in AVC in choirboys during the 10-year period 1994–2003. Ma et al. (Reference Ma, Chen, Chen, Zhu, Xiong, Li, Wang, Liu, Luo, Liu and Du2011) compared ages of spermarche in Chinese boys who participated in regular national surveys since 1979 and found that the spermarcheal age declined over the years. Overall, there has been a tendency toward a decline in age of male pubertal onset in recent years, although more data are necessary to draw a firm conclusion.
Twin studies in Western populations have found that about 45% (Snieder et al., Reference Snieder, MacGregor and Spector1998) to 95% (Loesch et al., Reference Loesch, Huggins, Rogucka, Hoang and Hopper1995) of the variance in age at menarche is attributable to genetic influences. In South Korean females, 72% (95% CI [67%, 76%]) of the variance of menarche was explained by genetic influences (Hur, Jin et al., Reference Hur, Jin and Lee2019). Although molecular genetic studies have shown that quite a few loci involved in pubertal timing are shared between two sexes (Cousminer et al., Reference Cousminer, Stergiakouli, Berry, Ang, Groen-Blokhuis and Körner2014; Day et al., Reference Day, Bulik-Sullivan, Hinds, Finucane, Murabito, Tung, Ong and Perry2015; Ong et al., Reference Ong, Elks, Li, Zhao, Luan, Andersen and Wareham2009), they have also found the presence of sex-specific alleles in pubertal timing (Day et al., Reference Day, Bulik-Sullivan, Hinds, Finucane, Murabito, Tung, Ong and Perry2015). However, as most twin and molecular genetic studies of pubertal timing to date have targeted females, the genetic architecture of male pubertal timing remains relatively poorly understood. Using two national samples of siblings from the Non-Shared Environment of Adolescent Development (NEAD) and the National Longitudinal Study of Adolescent Health (Add Health), Ge et al. (Reference Ge, Natsuaki, Neiderhiser and Reiss2007) have shown that genetic variances associated with male pubertal timing measured by the Pubertal Development Scale (PDS) were 45% (95% CI [15%, 62%]) in the NEAD sample and 40% (95% CI [18%, 55%]) in the Add Health sample, which were lower than the corresponding estimates in girls (50% vs. 46%). Variances due to shared environmental influences were close to zero for both sexes in both samples in the Ge et al. study. Similarly, in a Finnish adolescent twin sample, Mustanski et al. (Reference Mustanski, Viken, Kaprio, Pulkkinen and Rose2004) found that genetic variance accounted for 41% of the variance in male AVC and 59% of the variance in menarche (calculated from twin correlations shown in Table 3). Corley et al. (Reference Corley, Beltz, Wadsworth and Berenbaum2015) also examined scores in the PDS in participants in the Colorado Longitudinal Twin Study and reported that genetic influences on pubertal timing were 62% in males and 77% in females, while shared environmental influences were negligible in both sexes. Eaves et al. (Reference Eaves, Silberg, Foley, Bulik, Maes, Erkanli, Angold, Costello and Worthman2004) reported exceptionally high heritability estimates for the relative timing of pubertal change in participants of the Virginia Twin Study of Adolescent Behavioral Development. Specifically, the authors found that additive genetic factors were 96.3% in girls and 88.0% in boys, and nonshared environmental effects were 0.1% in girls and 12.0% in boys.
Taken together, studies of pubertal timing have consistently suggested that genetic influences in males are substantial (about 40%–88%), but they are somewhat lower than those found in females, and that as with females, nongenetic influences in male pubertal timing largely represent measurement error and nonshared rather than shared environmental influences.
Whereas lean body mass, which primarily reflects muscle mass, begins to increase during early puberty in both boys and girls, fat mass increases during the late stages of puberty only in girls (Wheeler, Reference Wheeler1991). Thus, an inverse relationship between body mass index (BMI) and age at menarche has been consistently reported across studies to date (Bratke et al., Reference Bratke, Bruserud, Brannsether, Aßmus, Bjerknes, Roelants and Júlíusson2017; Wang et al., Reference Wang, Dang, Xing, Li and Xing2016). Using a Finnish twin sample, Kaprio et al. (Reference Kaprio, Rimpelä, Winter, Viken, Rimpelä and Rose1995) demonstrated that common genetic influences were responsible for this relationship. Furthermore, GWAS has shown that earlier age at menarche is associated with the BMI-increasing alleles of several genetic loci (Cousminer et al., Reference Cousminer, Berry, Timpson, Ang, Thiering and Byrne2013). In males, however, the relationship between BMI and pubertal timing remains unclear (Wagner et al., Reference Wagner, Sabin, Pfaffle, Hiemisch, Sergeyev, Korner and Kiess2012). While the majority of phenotypic studies in males have suggested that earlier pubertal timing is associated with higher BMI (Juul et al., Reference Juul, Magnusdottir, Scheike, Prytz and Skakkebaek2007; Ong et al., Reference Ong, Bann, Wills, Ward, Adams, Hardy and Kuh2012), some studies have reported later pubertal onset in obese boys (Wang, Reference Wang2002; Kleber et al., Reference Kleber, Schwarz and Reinehr2011; Lee et al., Reference Lee, Kaciroti, Appugliese, Corwyn, Bradley and Lumeng2010). Using a longitudinal twin cohort from birth to age 18 years, Silventoinen et al. (Reference Silventoinen, Haukka, Dunkel, Tynelius and Rasmussen2008) found that pubertal markers (age at peak height velocity, age at onset of pubertal growth spurt) were negatively associated with childhood BMI in males, and that these associations were explained by common genetic factors. Although the sample size of the study was relatively small (99 monozygotic and 76 dizygotic twin pairs), the results demonstrated genetic regulation of male pubertal growth in the longitudinal context.
Given that BMI and AVC are both highly heritable traits (Eaves et al., Reference Eaves, Silberg, Foley, Bulik, Maes, Erkanli, Angold, Costello and Worthman2004; Silventoinen et al., Reference Silventoinen, Jelenkovic, Yokoyama, Sund, Sugawara, Tanaka and Kaprio2019), it is reasonable to assume that common genetic polymorphisms may influence the association between BMI and AVC. Genetic correlation refers to the amount of overlap in genetic influences on two phenotypes. In a recent GWAS on the basis of over 55,000 males of European ancestry, Day et al. (Reference Day, Bulik-Sullivan, Hinds, Finucane, Murabito, Tung, Ong and Perry2015) reported the genome-wide genetic correlation of −.28 between BMI and AVC. However, Cousminer et al. (Reference Cousminer, Stergiakouli, Berry, Ang, Groen-Blokhuis and Körner2014) showed that some BMI-increasing alleles were associated with earlier and others with delayed pubertal development in boys, suggesting that the relationship between pubertal timing in boys may be complex and needs more studies.
The aims of the present study were to investigate whether or not the mean level of AVC decreases in recent years, estimate heritability for AVC and determine the extent to which common genetic factors influence the relationship between AVC and BMI in South Korean male twins.
Methods
Sample and Measures
Subjects were male twins who participated in telephone surveys conducted in the South Korean Twin Registry (Hur, Kang et al., Reference Hur, Kang, Jeong, Kang and Kim2019). Twins aged 16 years or older at the time of the interview were chosen for the present study, as boys typically experience initiation of voice change by 15 years of age (Buck Louis et al., Reference Buck Louis, Gray, Marcus, Ojeda, Pescovitz and Witchel2008; Hollien, Reference Hollien2012). In total, 955 male twins met the age criteria. Of these, 73 (7.64%) reported that voice change had not occurred yet, while 118 (12.36%) reported that they could not recall their AVC. The total sample consisted of 600 monozygotic (MZ) twins (241 pairs and 118 co-twin missing twins), 214 dizygotic (DZ) twins (82 pairs and 50 co-twin missing twins) and 141 male members of opposite-sex DZ twins. The age of the sample ranged from 16.00 to 29.25 years with a mean of 18.92 (SD = 2.42) years. Twins’ birth years ranged from 1988 to 2001. Twins under 20 years of age were recruited mostly from schools throughout South Korea, while those aged 20 years or older were from Facebook and twin clubs on the Internet and colleges throughout South Korea. During the telephone interview, twins were asked to recall the age of the first onset of voice break. Specifically, the question asked was How old were you when you noticed your voice deepening? AVC was recorded in whole year. Additionally, twins were asked to report their current height (to the nearest centimeter) and weight (to the nearest kilogram). BMI was computed as weight (kg)/height (m)2.
Zygosity of the twins was assessed on the basis of three items asking about twins’ physical similarity and frequency of confusion by family members, close friends and strangers. When compared to DNA analysis, this approach has been shown to achieve over 90% accuracy (Ooki et al., Reference Ooki, Yamada, Asaka and Hayakawa1990). The number of MZ twins was much greater than that of DZ twins in the present sample, reflecting the twin birth rates in the South Korean population for the birth cohorts in the present study (Hur & Kwon, Reference Hur and Kwon2005; Hur & Song, Reference Hur and Song2009).
Statistical Analyses
To examine secular trends in AVC, the total sample was subdivided into two birth cohorts with intervals chosen to ensure that each cohort contained roughly equal numbers: twins born between 1988 and 1993 (n = 464, 48.59%), and those born between 1994 and 2001 (n = 491, 51.41%). The mean AVC and their 95% confidence intervals (CIs) for each cohort as well as for the total sample were estimated by the Kaplan–Meyer survival analysis. The Mantel–Cox rank test was used to compare the mean values of AVC between the two birth cohorts. As twins are not independent samples, the mean differences between two cohorts were tested in the first- and the second-born twins separately as well as in the combined sample.
Maximum likelihood twin correlation and univariate and bivariate model-fitting analyses were conducted to estimate genetic and environmental influences on AVC and BMI and genetic and environmental correlations between the two variables. Using univariate model-fitting analysis, the total variance of each variable was partitioned into four components: additive genetic (A; rMZ = 1.0, rDZ = 0.5), nonadditive genetic (D; rMZ = 1.0, rDZ = 0.25), shared environment (C; rMZ = 1.0, rDZ = 1.0) and nonshared environment plus measurement error (E; rMZ = 0, rDZ = 0). As C and D cannot be estimated simultaneously, ACE and ADE models were fit to the data separately. In bivariate model-fitting analyses, the association between AVC and BMI was decomposed into additive genetic, nonadditive genetic, shared and nonshared environmental factors to estimate correlations associated with each factor.
A raw data approach in Mx (Neale et al., Reference Neale, Boker, Xie and Maes2003) was used to perform model-fitting analyses. To determine the best-fitting most parsimonious model, all submodels were compared to the saturated model using the log-likelihood ratio test (LRT). In addition, Akaike’s information criterion (AIC = −2LL − 2df; Akaike, Reference Akaike1987) was used to evaluate alternative models when the models were not nested to each other. As lower AIC values indicate a better balance between parsimony and goodness-of-fit, the model with the lowest AIC was chosen as the best-fitting model. For twin correlational and model-fitting analyses, AVC and BMI were corrected for age and age2, and the residuals were standardized with a mean of zero and standard deviation of 1.0. Opposite-sex DZ twins were not included in twin correlations or model-fitting analyses.
Results
Secular Trends in AVC
AVC in the present sample was approximately normally distributed with a skewness of −0.07 and a kurtosis of 0.11. The mean value of the AVC in the total sample was 14.19 (95% CI [14.09, 14.29]) years. Table 1 presents the mean values of AVC in the first- and the second-born twins, and the total sample for the two birth cohorts (1988–1993, 1994–2001). As can be seen in the table, means were consistently higher in the older than in the younger cohort across the first- and the second-born twins and the combined sample, confirming a downward trend of mean AVC from 1988 to 2001. Mean differences between the two birth cohorts were significant in all three groups. In the total sample, the mean values of AVC were 14.38 (95% CI [14.25, 14.52]) years for the older cohort and 14.02 (95% CI [13.89, 14.16]) years for the younger cohort.
The mean AVC was not significantly different between MZ and DZ twins (14.22 vs. 14.15 years) or between the first- and the second-born twins (14.19 vs. 14.19 years).
Twin Correlations and Univariate Model-Fitting Analysis
Maximum likelihood MZ and DZ twin correlations were, respectively, .60 (95% CI [.50, .68]) and .25 (95% CI [.03, .45]) for AVC, and .84 (95% CI [.80, .87]) and .48 (95% CI [.31, .62]) for BMI. Table 2 shows the results of univariate model-fitting analysis. As indicated in Table 2, the AE model was best for both AVC and BMI, because differences in chi-square between the AE models and the saturated models were not significant and because AE yielded the lowest AIC among all submodels. In the best-fitting univariate model, additive genetic and nonshared environmental effects plus measurement error were, respectively, .59 (95% CI [.50, .67]) and .41 (95% CI [.33, .50]) for AVC, and .84 (95% CI [.81, .87]) and .16 (95% CI [.13, .19]) for BMI.
Note: −2LL, −2 log-likelihood; A, additive genetic effects; C, shared environmental effects; D, nonadditive genetic effects; E, nonshared environmental effects and measurement error. The fit indices of the submodels were compared to those of the saturated model. The best-fitting model is indicated in bold type.
Bivariate Model-Fitting Analyses
The age-corrected phenotypic correlation between BMI and AVC was modest but significant (r = −.14; 95% CI [−.21, −.07]), suggesting that earlier AVC was associated with higher BMI. Table 3 presents the results of bivariate model-fitting analyses. As expected from univariate model-fitting analysis, a nonsignificant change in −2LL occurred in AE (model 2) when C was removed from the saturated model. Removing either r a or r e individually from AE did not worsen model-fit significantly (models 3 and 4). However, when both r a and r e were eliminated from AE simultaneously, a significant change in −2LL occurred (model 5), indicating that either r a or r e should be maintained in the model. Thus, AIC values were compared between model 3 (without r a) and model 4 (without r e). AIC was much lower in model 4 than in model 3, suggesting that the AE model without r e is better than that without r a. In fact, model 4 had the lowest AIC in all submodels tested. Therefore, it was concluded from both LRT and AIC values that model 4 was the best-fitting one.
Note: −2LL, −2 log-likelihood; A, additive genetic effects; E, nonshared environmental effects and measurement error; r a, additive genetic correlation; r e, nonshared environmental correlation. The fit indices of the submodels were compared to those of the saturated model. The best-fitting model is indicated in bold type.
Figure 1 shows standardized parameter estimates in the best-fitting model. As indicated in Figure 1, the estimates of additive genetic and nonshared environmental influences on AVE and BMI were identical to those found in univariate model-fitting analysis. Additive genetic correlation between AVC and BMI was −.18 (95% CI [−.31, −.06]). Together with nonsignificant, nonshared environmental correlation, these results suggest that the phenotypic association between AVC and BMI, albeit modest, was entirely mediated by common genes.
Discussion
Three main findings emerged from the present study. First, South Korean males showed a downward trend in AVC from 1988 to 2001. Second, heritability for AVC in South Korean males was .59 (95% CI [.50, .67]), whereas shared environmental influences were close to zero. Finally, the phenotypic association between AVC and BMI was modest. However, this association was entirely due to common genetic influences.
The mean AVC in the total sample was 14.19 (95%CI [14.09, 14.29]) years. Consistent with recent studies of male pubertal timing (Juul et al., Reference Juul, Magnusdottir, Scheike, Prytz and Skakkebaek2007; Ma et al., Reference Ma, Chen, Chen, Zhu, Xiong, Li, Wang, Liu, Luo, Liu and Du2011; Sorensen et al., Reference Sorensen, Aksglaede, Petersen and Juul2010), however, the present sample showed a downward trend in AVC from 1988–1993 to 1994–2001 (14.38–14.02 years). Prior studies have shown downward trends in age at menarche in South Korean females as well (Lee et al., Reference Lee, Kim, Oh, Lee and Park2016; Hur, Jin et al., Reference Hur, Jin and Lee2019). Although puberty timing is highly heritable in both sexes, most researchers have speculated that the secular downward trends in pubertal timing are attributable to environmental factors, such as improved health, nutrition and endocrine-disrupting chemicals (Euling et al., Reference Euling, Herman-Giddens, Lee, Selevan, Juul, Sørensen, Dunkel, Himes, Teilmann and Swan2008). Thus, it would be of interest for future studies to explore whether environmental factors influencing the secular downward trends are shared between two sexes or sex-specific.
Heritability estimate for AVC in the present sample was comparable to that found in the Colorado Twin Study (Corley et al., Reference Corley, Beltz, Wadsworth and Berenbaum2015), but somewhat higher than that found in the Ge et al. (Reference Ge, Natsuaki, Neiderhiser and Reiss2007) or the Mustanski et al. (Reference Mustanski, Viken, Kaprio, Pulkkinen and Rose2004) studies, and much lower than that reported in the Eaves et al. (Reference Eaves, Silberg, Foley, Bulik, Maes, Erkanli, Angold, Costello and Worthman2004) study. However, consistent with these studies, heritability was lower as compared to that in female counterparts reported previously (.59 vs. .72; Hur, Jin et al., Reference Hur, Jin and Lee2019). These consistent observations in sex differences in heritability for pubertal timing may be in part due to the fact that measures of pubertal markers are less reliable in males than in females. However, given the presence of sex-specific genes in pubertal timing discovered in previous studies (Day et al., Reference Day, Bulik-Sullivan, Hinds, Finucane, Murabito, Tung, Ong and Perry2015), genetic influences on pubertal timing may be truly lower in males than in females. Additional studies are required to clarify sex differences as well as an overlap between sexes in genetic influences on pubertal timing.
Heritability estimate for BMI was .84 in the present sample, which was very close to the estimates found in previous analyses of South Korean twins (Hur et al., Reference Hur, Jin and Lee2018; Hur, Reference Hur2007). The genetic correlation between BMI and AVC in the present sample was modest (r = −.18), indicating that the amount of overlap in genes influencing BMI and AVC may be relatively small. However, bivariate model-fitting analyses showed that common genes are responsible for the association between BMI and AVC, similar to what has been found in females (Kaprio et al., Reference Kaprio, Rimpelä, Winter, Viken, Rimpelä and Rose1995). The present results are also consistent with findings from recent GWAS. For example, Ong et al. (Reference Ong, Elks, Li, Zhao, Luan, Andersen and Wareham2009) and Widén et al. (Reference Widén, Ripatti, Cousminer, Surakka, Lappalainen, Järvelin, Eriksson, Raitakari, Salomaa, Sovio, Hartikainen, Pouta, McCarthy, Osmond, Kajantie, Lehtimäki, Viikari, Kähönen, Tyler-Smith and Palotie2010) have reported that LIN28B on chromosome 6, a potent, and specific regulator of microRNA processing is significantly associated with both childhood growth and pubertal markers, including AVC.
The relationship between puberty timing and BMI is complex and still controversial among researchers (Day et al., Reference Day, Bulik-Sullivan, Hinds, Finucane, Murabito, Tung, Ong and Perry2015). Although negative genetic correlation has been discovered to explain the relationship between puberty timing and BMI, other biological mechanisms may be involved in the relationship. Several researchers have argued that hormonal mechanisms may mediate the relationship between BMI and puberty timing, because body fat percentage is a central regulator of the pubertal reactivation of the hypothalamic–pituitary axis (Ahmed et al., Reference Ahmed, Ong and Dunger2009). Specifically, Ahmed et al. (Reference Ahmed, Ong and Dunger2009) proposed that rapid weight gain during early infancy and childhood can trigger early pubertal development by increased exposure to sex steroids in prepubertal children. However, environmental agents such as endocrine-disrupting chemicals can also change endogenous hormonal milieu (Greenspan & Lee, Reference Greenspan and Lee2018), environmental influences can be disguised as genetic influences. Thus, further studies are required to explain genetic and environmental causes of the relationship between BMI and pubertal development.
Although the present study is the first to report genetic influences on AVC and genetic correlation between AVC and BMI in South Korean males, there are several limitations in the present study. First, as with many other studies of AVC, this study used self-reported AVC. Self-reported AVC is prone to recall error. However, as there is a relatively short interval between the AVC and the time of assessment in the present sample, the recall error may be relatively small. Second, only AVC was used to measure pubertal timing in the present study. Future studies should include other pubertal markers, such as pubic hair, skin change and testicular volume to enhance our understanding of the genetic architecture of male pubertal timing. Third, the present study examined trends of AVC during a relatively short period of time (1988–2001). In the last century, rapid changes occurred in nutritional and other environments in South Korea. Future studies should extend the time period to better capture secular trends in AVC in South Koreans. Fourth, the sample size in the present study was relatively small for model-fitting analyses. As indicated in the full ADE model in Table 1, there was some hint of nonadditive genetic influence on AVC. However, the present study does not have sufficient power to estimate additive versus nonadditive genetic influences separately. Replication of the present findings with a larger sample is, therefore, necessary. Finally, the present sample comprised only of South Korean males, and thus, caution has to be exercised when the results are generalized to other ethnic groups, especially given that ethnic differences in pubertal timing have been reported (Parent et al., Reference Parent, Teilmann, Juul, Skakkebaek, Toppari and Bourguignon2003).
In summary, the present twin study demonstrated that genetic influences on AVC are significant and substantial in South Korean males, and that common genes are responsible for the association between AVC and BMI.
Acknowledgments
The author would like to thank twins who participated in the South Korean Twin Registry. This study was supported by the National Research Foundation of Korea grant (NRF-371-2011-1 B00047) and by the ‘Development of Health Prediction Technology based on Big Data’ (K18092) funded by the Ministry of Science and ICT(MSIT) of Korea given to the Korea Institute of Oriental Medicine (KIOM).