Introduction

Puberty refers to the re-activation of the hypothalamic-pituitary-gonadal axis, culminating in sexual maturation (Sisk and Foster 2004; Ebling 2005). The timing of the gradual increase in gonadotropin-releasing hormone (GnRH) pulse frequency at the beginning of puberty is complex, with many systems involved (Ojeda et al. 2006). There are substantial individual differences in the timing of puberty and it is likely that genetic factors interact with environmental factors. Several studies have indicated that individual differences in the timing of pubertal development are heritable. Depending on method of assessment and phenotype definition, heritability coefficients range from 50 to 80% (Palmert and Hirschhorn 2003; Meyer et al. 1991; Loesch et al. 1995).

Menarche is closely related to breast development: largely the same genetic effects seem to be involved in the timing of both phenotypes (van den Berg et al. 2006; Pickles et al. 1998). This makes menarche a valid proxy for studying the genetic background of pubertal timing. In a sample of 12-year-old female twins where 14% already had had their first menses, application of a multifactorial threshold model resulted in an estimate of 30% for the heritability of menarche before age 12 (van den Berg et al. 2006). Heritability studies on age at menarche, most of them retrospective and only using adult twin samples, usually find estimates of about 50% or higher (Towne et al. 2005). When controlling for measurement error, Pickles et al. (Pickles et al. 1998) found that nearly all ‘true’ variance in age at menarche could be attributed to additive genetic effects.

A drawback of these studies is that they did not go beyond the application of a relatively simple genetic model, that is, they examined neither assortative mating nor genotype–environment interaction. If there are interactions between genotype and environmental factors that are either shared or non-shared between family members, the heritability coefficient is biased when these interactions are ignored. For instance, the variance attributable to an interaction between additive genetic effects and non-shared environmental effects, when ignored, is attributed to non-shared environmental factors, while it might make more sense to attribute this variance to genetic factors (Falconer and Mackay 1996). Interaction between genetic effects and non-shared environmental effects can be investigated by correlating the difference between the two individuals from a monozygotic (MZ) twin pair with the average age at menarche (or, equivalently, the sum of the ages) of that twin pair (Jinks and Fulker 1970). With a positive correlation, the relative impact of the environment is larger with individuals with a familial predisposition for late menarche.

Another commonly made but untested assumption in studies estimating heritability is that mating is random with regard to the phenotype of interest: one assumes that mates are not genetically correlated regarding age at menarche. Of course, males do not have menses, but they are nevertheless carriers of the genes that influence the age at menarche in their female offspring. If males from families with precocious females have a tendency to mate with relatively precocious females, and males from families with late puberty onset choose to mate with relatively late females, the offspring in one family is expected to be more genetically similar with respect to age at menarche than under random mating (Falconer and Mackay 1996). The assortment need not take place based on the phenotype per se: a genetic correlation is also induced when assortment takes place on a phenotype genetically correlated with age at menarche, for example BMI (Kaprio et al. 1995).

An increase in the genetic similarity of the offspring has important consequences for quantitative genetic modeling. Usually, it is assumed that the additive genetic effects affecting a phenotype correlate 1/2 in full sisters and female dizygotic (DZ) twins. Under positive phenotypic assortment (males mating with females with puberty timing similar to that of their mothers and sisters), the genetic correlation is larger than 1/2. When ignored, heritability is underestimated.

The familial clustering of age at menarche was studied in a large sample of female twins, their mothers, their full sisters and female spouses and first-degree relatives of male twins. All were asked to indicate their age at menarche as part of a longitudinal survey study.

Subjects and methods

Subjects

Analyses were based on data from multiples and their family members who are registered with the Netherlands Twin Registry (NTR). In 1991, the NTR started a longitudinal survey study of health, personality and lifestyle. Surveys were sent out in 1991, 1993, 1995, 1997, 2000 and 2002 to adolescent and adult twins and their family members. Twin pairs were asked to participate in all waves, siblings were included since 1995. Parents of twins were asked to participate in 1991, 1993, 1995 and 2002 and spouses (married or unmarried) since 2000. Families with adolescent and young adult twins were recruited through City Councils in 1990–1991 and in 1992–1993. After 1993 an effort was also made to recruit adult and older twins through a variety of approaches. Further details on response rates, response bias and demographic characteristics of the sample can be found elsewhere (Boomsma et al. 2002; Stubbe et al. 2005; Vink et al. 2004; Koopmans et al. 1999). Data from female spouses and first-degree female relatives of male twins were used to examine assortative mating. For the genetic model fitting, data from female-female twin pairs, female twins from opposite sex twins, female full siblings and mothers were used; data from half-siblings and multiples other than twins were excluded from the analyses. We had data on 1,340 families with MZ female twins, 793 families with DZ female twins, 1,078 families with DZ unlike sex twins and 1,141 families without female twins or where zygosity was unknown (16 families). We had data on 4,995 individual twins, 1,296 sisters, 2,946 mothers and 635 female spouses of male twins. Average family size (excluding spouses) was 2.1 (SD = 1.0; median = 2, mode = 1, maximum = 7).

For 993 female twin pairs zygosity was based on DNA polymorphisms. For the other same sex twin pairs, zygosity was based on eight items on physical similarity and the frequency of confusion of the twins by parents, other family members and strangers. Agreement between zygosity based on these items and zygosity based on DNA was 98% (Willemsen et al. 2005). Average age at the time of the first available report on age at menarche was 25 years for twins, 30 years for siblings and 46 years for mothers.

Study design

In 1991, 1993 and 1995 female participants were asked to indicate the age at menarche in years and months. In 2000 and 2002, participants were only asked to indicate the age in years. In addition, in 1993, 1995 and 2002, participants were asked whether they had had their first menstruation prior to indicating the age at which it had occurred. All data concerning age at menarche were rescaled to number of months, where 6 months were added to the data from 2000 and 2002. This way we avoided bias due to the fact that we only had data on the age in years.

Data-analysis

Over the years, the reports on menarche were not entirely consistent. Data were discarded when the reported age at menarche was higher than the age at the time of the questionnaire and when one of more than two reports deviated more than 12 months from the other reports when these were more consistent. Data points were also ignored when participants indicated they had not yet had their first menses but nevertheless reported an age. Total discarded data was less than 1%.

Age at menarche was correlated between spouses of male twins and first-degree relatives of the male twins. Sums and differences in age at menarche in MZ twins were correlated to test for genotype–environment interaction. The covariance structure of age at menarche in mothers, siblings, MZ and DZ twins was assessed and subsequently modeled to estimate heritability in a quantitative genetic analysis. Model fit was assessed by comparing the fit of nested models using likelihood ratio tests. The probability for a Type I error was fixed at 5% for each test.

Results

Consistency

The within-person correlations over time ranged from .71 to .89. The first available report correlated .97 with the average reported age at menarche. The first available report was considered to be the most reliable and used for the analyses. Average age at menarche was 161.7 months (13.5 years, SD = 17.0, minimum = 96, maximum = 246).

Assortative mating

The correlation between the age at menarche in female spouses of male twins is −.08 (N = 88, 95% CI −.28, .12), suggesting that partner choice with regards to pubertal timing is not based on genetic or environmental factors that are shared in (male) twins. Correlations between age at menarche of mothers and that of the spouses of their male twin offspring was estimated at .08 (N = 345, −.02, .18). Correlations between spouses of twins and sisters was −.14 and just significant (N = 229, −.26, −.005) and the correlation between a male twin’s spouse and his female co-twin was −.08 (N = 165, −.22, .06). The inconsistent pattern of correlations between spouses and first-degree relatives—negative for siblings, positive for mothers—and no correlation between the female spouses of twins suggests no assortment for age at menarche or traits related to it. Therefore in the genetic modeling, we assumed that mating is random in the population with regards to age at menarche and that the correlation between any additive genetic effects in female offspring is 1/2.

Gene–environment interaction

There was a significant correlation between the sums and absolute differences in MZ twins, r = .13, N = 1,122, p < .01, suggesting an interaction between shared genetic or environmental effects and non-shared environmental effects (Jinks and Fulker 1970). The positive sign of the correlation indicates that non-shared environmental effects seem to have more impact in twins that have a familial predisposition for late menarche. Of course, these non-shared environmental effects might simply reflect error variance: possibly, there is more unreliability in the reports of females with late menarche.

When analyzing the MZ twin data from each measurement separately, we see the lowest correlation between MZ sums and absolute differences in 1991 (r = −.01, N = 346) and the highest in 2002 (r = .24, N = 577). The variance of the absolute difference scores in 2002 is 144, while the within-twin absolute differences for 2002 and 2000 reports have variances of 59 and 48, for twin 1 and twin 2, respectively. Thus, the larger part of the environmental variance seems to be true environmental variance, rather than measurement unreliability.

In 2002, we found no correlation between reported age at menarche and absolute difference between the 2000 and 2002 reports for the first twin (r = .07, N = 513, p = .13), but a significant correlation for the second twin (r = .16, N = 524, p < .001). The 2002 twin 1 + twin 2 summed report correlated neither with the within-twin 1 differences, r = .07, p = .19, nor with the within-twin 2 differences, (r = .08, p = .09). But again, the correlation between the absolute difference between twins in 2002 with within-twin 1 2002–2000 difference was significant (r = .17, p = .001), and similarly for twin 2 difference (r = .14, p = .006). Taken together, these results suggest that in 2002, the environmental variance in part reflects error, but whether the error is correlated with the actual age at menarche remains unclear. It seems that the significant correlation between twin sums and differences in 2002 cannot be explained, at least not entirely, by measurement unreliability in females with late menarche.

It should be noted that the GE interaction reflected by the correlation between sum and differences is small, explaining only 1.5% of the variation in the reports. Whatever the extent to which this genotype-dependent environmental variance reflects true variance, this interaction component was ignored in further analyses.

Analysis of means, variances and covariances

A base model was constructed specifying for each of the four different types of families three mean values and three variances for twins, siblings and mother. In addition for each group four covariances were specified for mother–daughter, sibling–sibling, sibling–twin and twin–twin relationships (see Table 1). The twin–twin covariance was not estimated for the opposite sex twin group, nor was the twin–sibling covariance in the male twins/unknown zygosity group as there were no data for this relationship. The -2LL fit statistic was 77233.67 with 9,195 degrees of freedom. In consecutive steps, mean values, variances and covariances were equated across groups, and also within groups, leading to a reference model where all mean values were equal, twin variances were equal, sibling variances were equal and mother variances were equal across twin zygosities. Twin, sibling and mother variance could however not be equated without significant loss of goodness of fit. Since siblings were on average older and the mothers obviously even more, the differences in variance could be due to a memory problem, rendering the reports of siblings and mothers less reliable than in twins. The variance in the 50% oldest mothers (326.09) was significantly larger than the variance in the 50% youngest mothers (290.35), Levene’s test F(1, 2,944) = 3.71, p ∼ .05. Moreover, sibling reports were not available from the 1991 and 1993 surveys and were therefore more often than for twins based on reports from the 2000 and 2002 surveys, where age at menarche was only reported in years. This lack of precision might well have increased measurement error. Indeed, in the twin data for example, the variance in 2000 and 2002 was 307 and 312, respectively, whereas in the preceding waves the variance ranged between 222 and 251. Moreover, MZ twin correlations were .60 and .57 in 2000 and 2002, and between .81 and .87 in the preceding waves.

Table 1 Observed covariance (below diagonal) and correlation matrices (above diagonal) for the different types of families (families with unknown twin zygosity not presented)
Full size table

In further analyses, therefore, it was assumed that the surplus sibling and mother variances were due to extra measurement variance. There was no evidence for a systematic telescoping effect where the age at menarche is reported more forward or more backward in time with increasing duration between event and report: means of twins, siblings and mothers could be equated without significant loss of fit. The correlation between age at menarche and age at the time of report was only .06. The correlation between birth year and reported age at menarche was -.06. The distribution in twins and siblings had however significant excess kurtosis, being more leptokurtotic than expected under the assumption of normality. This was not observed in the mothers where the distributions for old and young mothers were slightly platykurtotic but not significantly so.

The sibling–sibling covariance, twin–sibling covariance and DZ twin covariance could all be equated. These were however significantly smaller than the MZ twin covariance, χ 2(1) = 202.16, suggesting a genetic component in the individual differences. The reference model had a −2LL fit statistic of 77267.43 with 9,229 degrees of freedom. The covariance structure estimated using this reference model is given in Table 2. The DZ twin/sibling covariance is a little less than half the MZ twin/sibling covariance, but the mother–daughter covariance is practically half the MZ covariance; taken together this suggests no non-additive gene action.

Table 2 The expected covariance structure based on the reference model (top) and the genetic model (bottom)
Full size table

To put this interpretation to the test, a genetic model was fitted using maximum likelihood estimation on the raw data, hypothesizing additive genetic variance (A), dominance genetic variance (D) and non-shared environmental variance (E), in addition to surplus variance in siblings (S) and mothers (M). The expectation for twin variance was A + D + E, for siblings, A + D + E + S and for mothers A + D + E + M. Since MZ twins are genetically identical, the expectation for the MZ twin covariance was A + D. For DZ twin and siblings covariances the expectation was 1/2A + 1/4D, and for the mother–daughter covariance 1/2A, since parents only transmit additive genetic variance. The model had a −2LL fit statistic of 77268.09 with 9,230 degrees of freedom. The estimates for A, D and E were 189.42, 0.65 and 82.02, respectively, with D not significantly different from 0, Δ − 2LL < 0.01. S and M were 25.88 and 43.91, respectively. Adjusting for the surplus error variance in siblings and mothers, broad heritability is estimated at (A + D)/(A + D + E) = 70% (95% CI 67%, 72%). With the constraint D = 0 the results for narrow heritability were exactly the same.

Discussion

The age at menarche in sisters of male twins correlated significantly with the age at menarche of the spouses of male twins. This was a negative correlation, where a late menarche in the male twin’s sisters predicted early menarche in the male twin’s spouse. However at the same time, the correlation between the age at menarche in the mother and the male twin’s spouse was positive. Thus, there is no consistent evidence that the genetic effects that affect age at menarche are correlated in mating partners. This is important since in quantitative genetic modeling it is often assumed that this is not the case. If a positive genetic correlation between spouses truly exists but is ignored, heritability is underestimated. In the case of age at menarche it can be safely assumed that the additive genetic effects on age at menarche correlate 1/2 in sisters.

Based on the quantitative genetic modeling, significant evidence of additive genetic contribution to individual differences in age at menarche was found, accounting for 70% of the variance in the reports. Non-additive gene action was not detected. If there has been selection for the trait over the course of evolution, it probably was selection for an intermediate optimum.

It should be noted however, that the modeling was based on the assumption that the heritability was constant across generations. Furthermore, it was assumed that there are no environmental influences that are shared among offspring (e.g., maternal effects). Given the present study design it is not possible to disentangle dominant gene action, epistatic effects and environmental effects common to all offspring. In order to disentangle the effects of non-additive genetic effects and environmental effects shared in siblings raised together, data from siblings reared apart could be used. Since results from other studies using only twins suggest shared environmental effects on menarche (van den Berg et al. 2006), dominant and epistatic gene action may therefore have been underestimated in this study. Treloar and Martin (1990) for instance reported evidence of dominant gene action. In the present study, there was no indication of environmental effects shared by mothers and their offspring.

There was a significant correlation between sums and absolute differences in age at menarche in MZ twins. Under the assumption that all variance shared in MZ twins is genetic we may conclude that the impact of environmental factors tends to be larger for those individuals that are genetically predisposed to late menarche. Earlier it has been established that prenatal programming takes place regarding postnatal GnRH secretion under the influence of sex steroids (Sisk and Foster 2004). To the extent that prenatal steroid levels and their effect partly depend on random environmental factors, this possibly reflects the gene–environment interaction effect found here. Late menarche might be the result of the non-additive effects of both a genetic predisposition and random prenatal factors affecting steroid levels. Alternatively, late menarche might result from the interplay between a genetic predisposition and environmental factors affecting prenatal growth and birth weight (Ibanez et al 2000). The interaction may also be due to postnatal factors, where individuals genetically predisposed for late menarche are more under the influence of environmental menarche triggering events such as leptin, glucose and insulin levels and metabolic fuel availability, olfactory cues, diet and stress (Sisk and Foster 2004; Ebling 2005). However, the gene–environment interaction effect was relatively small, accounting for only 1.5% of the variance in twins. Moreover, it cannot be ruled out that at least part of the GE interaction variance reflects error variance associated with reports of late menarche. Whatever the case, the bulk of the variance seems to result from the effects of genes and environment that act independently, with genes explaining most of the variance.

Our results showed an estimate of 70% for the amount of total variance accounted for by additive genetic variance. Since in the analysis gene–environment interaction was ignored, its variance is attributed to non-shared variance. One may therefore conclude that between 70 and 71.5% of the total variance can be attributed to the effects of genotype (Falconer and Mackay 1996).

Total variance however includes variance due to measurement error. Therefore the relative impact of genetic influences may be underestimated. Given the average ages of twins, siblings and the mothers, there was a long interval between the actual menarche and the time of report that may have resulted in substantial error variance. To some extent the error variance was accounted for by allowing larger variances for siblings and mothers, but not all error variance could be accounted for. Although there is indeed evidence of unreliability in retrospective menarche reports, the error seems to be symmetric around the actual age at menarche since the mean values in the three groups were equal (Must et al. 2002; Greif and Ulman 1982; Pillemer et al. 1987). We did however observe leptokurtic distributions in the reports. These might be indicative of a tendency to report menarche more closely to the average age than expected with a normal distribution, possibly related to a memory bias towards the population mean. A leptokurtic distribution for recalled age at menarche has been observed in other studies (Treloar and Martin 1990; Zacharias et al. 1970). However the present results should be interpreted with care, since the test is highly sensitive to the assumption of a continuous variable, which is not the case here where in two surveys one could indicate age at menarche only in years.

We found no systematic telescoping effect, where remembered events are placed more forward or more backward in time with increasing time between actual event and recall (Pickles et al. 1998). However, it should be noted such an effect is confounded with the generally observed decrease in age of puberty onset (Wyshak and Frisch 1982). We found only a very small correlation between reported age at menarche and age at the time of report. Whether this reflects a systematic telescoping effect, an actual decrease in age at menarche in the Dutch population or a combination of both, is uncertain.

Apart from error due to imperfect memory, the type of questions—asking for an age in years rather than years and months—formed a second source of error variance. Although we controlled for surplus error variance in siblings and mothers, being on average older than the twins and more often reported only on the later surveys, the error variance remaining could not be quantified and is included in the non-shared environmental effects. It could well be that nearly all variance is due to genetic effects (Pickles et al. 1998). It has been observed that true age at menarche correlates .79 with the age reported 30 years later (Must et al. 2002). When comparing the 70–71.5% heritability estimate with this reliability measure, one may conclude that actually about 90% of the true variance can be attributed to genetic effects (Schmidt and Hunter 1996).

This study showed high stability of menarche reports and relatively high heritability. This, coupled with results from other studies, suggests that recalled age at menarche is a fairly reliable trait that can be used in linkage and association studies.