Puberty is a pivotal life transition marked by complex processes involving physical and emotional changes (Ge, Conger, & Elder, Reference Ge, Conger and Elder2001; Steinberg, Reference Steinberg1987). Although research on the timing of pubertal onset, including age of menarche, reveals it to be heritable (Rowe, Reference Rowe, Rodgers, Rowe and Miller2000), it is also subject to a range of environmental influences, including nutrition, family structure and dynamics, and chronic stress (Ellis, Reference Ellis2004). It is also noteworthy that the timing of puberty has consequences for later life health and reproductive behavior. In the case of females, whose pubertal development is the focus of this report, early age of menarche forecasts varied outcomes across a range of developmental domains, including early sexual debut and pregnancy, greater sexual risk taking (Belsky, Steinberg, Houts, & Halpern-Felsher, Reference Belsky, Steinberg, Houts and Halpern-Felsher2010), early initiation of substance use (Deardorff, Gonzales, Christopher, Roosa, & Millsap, Reference Deardorff, Gonzales, Christopher, Roosa and Millsap2005; Dick, Rose, Kaprio, & Viken, Reference Dick, Rose, Kaprio and Viken2000), mental health problems (Deardorff et al., Reference Deardorff, Hayward, Wilson, Bryson, Hammer and Agras2007; Mendle, Harden, Brooks-Gunn, & Graber, Reference Mendle, Harden, Brooks-Gunn and Graber2010), unhealthy weight gain (Bratberg, Nilsen, Holmen, & Vatten, Reference Bratberg, Nilsen, Holmen and Vatten2007), and an elevated risk of cardiovascular disease and reproductive cancers in later adulthood (Kelsey, Gammon, & John, Reference Kelsey, Gammon and John1993; Lakshman et al., Reference Lakshman, Forouhi, Sharp, Luben, Bingham and Khaw2009; Vo & Carney, Reference Vo and Carney2007).
Psychosocial Acceleration Theory
With respect to family influences, another focus of this report, Belsky, Steinberg, and Draper (Reference Belsky, Steinberg and Draper1991) advanced what has come to be known as psychosocial acceleration theory (Belsky, Reference Belsky2012; Ellis, Reference Ellis2004). This evolutionary theory of socialization stipulated that experiences within the family having to do with psychological, emotional, and relational stress versus support regulate pubertal development in the service of reproductive goals. Psychological and behavioral processes involving attachment, parenting, and marital relations during the first 5–7 years of life were highlighted as being developmentally significant.
More specifically, children have evolved to regard a family environment marked by harsh, inconsistent, or otherwise unsupportive parenting and family relations as a signal or cue that the world they are likely to encounter in the future will be similar and thus that others cannot be trusted; that relationships may be exploitative and pair bonds unstable; and that life itself could be precarious. As a result, development has evolved to respond in an accelerated fashion in the face of such cues in order to reduce the likelihood of an individual dying before reproducing and thus of increasing the chances of passing on genes to the next generation. For these reasons, a child growing up under the conditions just delineated and with the often unconscious understandings just outlined should develop a reproductive strategy characterized by earlier sexual debut, unstable pair bonds, increased number of offspring, and decreased parental investment. In contrast, a child growing up in a stable and supportive family environment should develop in the opposite fashion. In both cases and in accord with evolutionary life-history theory, development is presumed to have been strategically shaped by natural selection to maximize reproductive fitness given the contextual stressors and supports the child encounters. All this is not to say that in modern times the evolved processes just outlined continue to enhance reproductive fitness, only that environmentally sensitive development has evolved in the manner described for these reasons.
Since the emergence of psychosocial acceleration theory, a great deal of research has addressed the puberty prediction unique to this theory of socialization (Belsky, Reference Belsky2012; James & Ellis, in press). The available research indicates that a variety of family stressors, including parenting and parent–child related processes, experienced early in life are associated with earlier age of menarche, such as father absence (Ellis et al., Reference Ellis, McFadyen-Ketchum, Dodge, Pettit and Bates1999), parent–child conflict (Moffitt, Caspi, Belsky, & Silva, Reference Moffitt, Caspi, Belsky and Silva1992), separation from family (Pesonen et al., Reference Pesonen, Räikkönen, Heinonen, Kajantie, Forsén and Eriksson2008), family disruption and father social deviance (Tither & Ellis, Reference Tither and Ellis2008), child maltreatment (Costello, Sung, Worthman, & Angold, Reference Costello, Sung, Worthman and Angold2007), marital conflict (Saxbe & Repetti, Reference Saxbe and Repetti2009), and sexual abuse (Wise, Palmer, Rothman, & Rosenberg, Reference Wise, Palmer, Rothman and Rosenberg2009). No one has shown to date that maternal sensitivity (actually insensitivity) plays a role in accelerating pubertal development. Thus, by focusing on this well-studied parenting construct in the work to be reported, we seek to extend research on environmental influences on pubertal development.
Differential-Susceptibility Theory: Alternative and Conditional Strategies
Further reflection on human evolution and the nature of development led Belsky (Reference Belsky, Rodgers, Rowe and Miller2000) to revise his thinking about whether rearing experiences should regulate the development of reproductive strategies, including pubertal timing, a line of theorizing about variation in susceptibility to environmental influence that has now been applied more generally to the study of human development (Belsky, Reference Belsky, Ellis and Bjorklund2005; Belsky, Bakersman-Kranenburg, & van IJzendoorn, Reference Belsky, Bakermans-Kranenburg and van Ijzendoorn2007; Belsky & Pluess, Reference Belsky and Pluess2009, Reference Belsky and Pluess2013; Ellis, Boyce, Belsky, Bakermans-Kranenburg, & van IJzendoorn, Reference Ellis, Boyce, Belsky, Bakermans-Kranenburg and van Ijzendoorn2011). Coming to appreciate that the future was inherently uncertain and thus that early experience could not always have succeeded in the course of human evolutionary history in matching development to future, within-generation contextual conditions (Belsky, Reference Belsky1997), Belsky (Reference Belsky, Rodgers, Rowe and Miller2000) applied differential-susceptibility theorizing to psychosocial acceleration theory. Central to such theorizing is the notion that there are individual differences in developmental plasticity, and as a result, children are not equally susceptible to the early environmental exposures that have long been thought to shape psychological and behavioral development (e.g., maternal depression, maternal sensitivity, and father involvement).
Such an orientation thus led to the hypothesis that there should be variation in susceptibility to the very socialization influences central to psychosocial acceleration theory. Whereas the rate of development, referring here to the time of sexual maturation, would be faster (occurring earlier) or slower (occurring later) depending on experiences while growing up in the case of “conditional” reproductive strategies, it would be far less a function of developmental experience in the case of “alternative” reproductive strategies (Belsky, Reference Belsky, Rodgers, Rowe and Miller2000). In the parlance of evolutionary, life-history theorists, then, conditional strategists adjust their behavior and/or development depending on the contextual conditions to which they are exposed, whereas alternative strategists do not, adopting a more or less fixed way of functioning regardless of contextual conditions. It is important to appreciate that alternative strategists are not all the same in terms of their rate of development. Thus, some will develop more quickly (i.e., sexual maturation occurring earlier) and some more slowly (i.e., sexual maturation occurring later), but this will not be in response to environmental signals or cues.
Two recent studies provide evidence consistent with this differential-susceptibility related revision of psychosocial acceleration theory. The first, focused on the timing and tempo of pubertal development, was conducted to test Boyce and Ellis's (Reference Boyce and Ellis2005) biological sensitivity to context proposition that children who are more physiologically reactive (e.g., show greater cortisol response to a stressor) would prove more susceptible to environmental influences; the second, focused on age of menarche, examined select Gene × Environment (G × E) interactions. In the first case, Ellis, Shirtcliff, et al. (Reference Ellis, Shirtcliff, Boyce, Deardorff and Essex2011) observed that limited parental supportiveness, measured during preschool, predicted early onset and faster pace of initial pubertal development, with the reverse being true when parenting was highly supportive, but that such rearing effects proved most pronounced and, consistent with psychosocial acceleration theory, in the case of children who scored high in physiological reactivity.
The second study, which directly informs the current work, examined the role of two single nucleotide polymorphisms (SNPs) of the estrogen receptor-α gene (ESR1) in moderating the effect of family cohesion and conflict (Manuck et al., Reference Manuck, Craig, Flory, Halder and Ferrell2011). These genes were selected because of evidence implicating them in the pubertal development process. More specifically, the maturation of the hypothalamic–pituitary–gonadal axis occurs through gradual increases in pulsatile release of pituitary-derived gonadotropins (follicle-stimulating hormone and luteinizing hormone) and heightened tissue exposure to estrogen (DiVall & Radovick, Reference DiVall and Radovick2008, Reference DiVall and Radovick2009). The effects of estrogen are primarily mediated through activation of intracellular estrogen receptors, so variations in estrogen receptor genes likely contribute to pubertal development. The two particular ESR1 polymorphisms, rs9340799 and rs2234693, which are the focus of this report, have been examined previously in studies involving pubertal timing as well as in relation to other estrogen-dependent conditions (e.g., breast cancer; see Manuck et al., Reference Manuck, Craig, Flory, Halder and Ferrell2011).
Consistent with prior research, Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) found that women reporting less cohesion in their families of origin during childhood had their first period at a younger age than did women who experienced more cohesive family environments while growing up. Of special importance to the G × E research reported herein, however, was that these investigators found that such apparent environmental effects were moderated by two polymorphic variations in ESR1 (rs9340799 and rs2234693). For women who were homozygous for the either ESR1 minor alleles, GG for rs9340799 or CC for rs2234693, quality of family environment predicted age of menarche in a manner consistent with psychosocial acceleration theory, but no such effect emerged in the case of women of other genotypes. That is, only when women carried particular “plasticity alleles” did less supportive family environments predict earlier age of menarche and more supportive environments later age of menarche.
A major and explicitly acknowledged limitation of the Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) G × E study was its reliance on retrospective measurements of the critical predictor and outcome variables. This led the investigators to call for further work to determine whether similar results would emerge when prospective data were available, appreciating as they did that memory biases could result in inaccurate reports about experiences earlier in life. Thus, here we seek to determine, using data collected as part of the National Institute of Child Health and Human Development (NICHD) Study of Early Child Care and Youth Development (SECCYD), whether the same two polymorphisms studied by Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) moderate the effect of repeatedly observed maternal sensitivity on age of menarche, in a manner consistent with the Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) results and the differential-susceptibility revision of psychosocial acceleration theory. This issue is especially interesting in light of prior work with a NICHD study that failed to detect a main effect of maternal sensitivity on pubertal development (Belsky, Steinberg, et al., Reference Belsky, Steinberg, Houts, Friedman, DeHart and Cauffman2007). To the extent that genetic moderation of maternal-sensitivity effects emerge in this inquiry, it will serve as reminder that G × E interactions (or any interaction, for that matter) can emerge even in the absence of main effects.
Diathesis–stress versus differential susceptibility
In addition to seeking to replicate the Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) findings using a measure of the childrearing environment obtained through observational assessment well before puberty, we sought to extend research in a second way, by testing alternative and competing models of G × E interaction: diathesis–stress and differential susceptibility. This is an important issue because, as with many G × E findings reported in the literature, Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) conducted no formal statistical test beyond the general test of interaction in the standard regression analysis. Central to our approach developed by Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012; Belsky, Pluess, & Widaman, Reference Belsky, Pluess and Widaman2013) is appreciation that the diathesis–stress and differential-susceptibility frameworks lead to a key differential prediction, which can be evaluated statistically, regarding the nature of the predicted G × E interaction.
Under diathesis–stress theorizing, the predicted interaction should be ordinal in form. Consider a biallelic polymorphism with three possible genotypes, containing zero, one, or two putative risk alleles. (The terminology of “risk” is based on the diathesis–stress notion that some alleles foster “vulnerability” to “adversity.”) According to diathesis–stress, a regression model with a Linear G × Linear E interaction should reveal four outcomes: a small or zero effect of the environment for the (resilient) group with zero risk alleles; a stronger, significant effect of the environment for the group with two risk alleles; a middling outcome by the group with one risk allele; and a crossover point of the linear functions at or near the most positive value for the environment.
Differential susceptibility leads to a contrasting prediction regarding the form of the G × E interaction. The alternate alleles under differential susceptibility are recast as plasticity and nonplasticity alleles, rather than risk and resilience alleles, respectively, because plasticity alleles are presumed to make individuals particularly susceptible to both negative and positive environmental effects (i.e., not just negative ones). The predicted interaction would still have a small (or nil) effect of the environment for the least plastic group, a stronger, significant effect of the environment for the plastic group, and a moderate effect for the moderately plastic group. However, the crossover point of these three linear functions would be near the middle of the distribution of scores on the environmental variable, thus revealing a “for better and for worse” pattern (Belsky, Bakermans-Kranenburg, et al., Reference Belsky, Bakermans-Kranenburg and van Ijzendoorn2007), with “better” outcomes (i.e., later age of menarche) predicted for the most plastic group under more favorable environmental conditions and “worse” outcomes (i.e., earlier age of menarche) for the most plastic group under less favorable ones.
The location of the crossover point for the predicted outcomes is therefore the crucial parameter that distinguishes predictions for the G × E interaction for the competing diathesis–stress and differential-susceptibility positions. Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) proposed a reparameterized regression model that makes the crossover point one of the parameters to be estimated. One major benefit of the reparameterization is that the point estimate of the crossover point is accompanied by a standard error, so that an interval estimate can be calculated. Among other things, the reparameterized model allows model fit under differential-susceptibility and diathesis–stress conditions to be statistically contrasted, with the better fitting model offered as the optimal representation of the data.
Strong versus weak diathesis–stress and differential-susceptibility models
Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012; Belsky et al., Reference Belsky, Pluess and Widaman2013) highlighted four reparameterized models that can provide tests of key parameters consistent with (a) weak and (b) strong differential-susceptibility and (c) weak and (d) strong diathesis–stress predictions without requiring an omnibus test of a G × E interaction before proceeding to examine the form of the interaction. Whereas strong models presume that some individuals are not at all susceptible to environmental effects (i.e., zero-order association between predictor and outcome), weak models presume that all are susceptible but that some are more so than others. Based on results of recent studies favoring strong models (Belsky et al., Reference Belsky, Pluess and Widaman2013; Widaman et al., Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012), we concentrated on distinguishing between the fit of the strong differential-susceptibility and diathesis–stress models, although we conducted additional tests to determine whether the weak versions of each model would improve the fit of models to data. Because they never did, only strong-model findings are reported; the results of the additional analyses are available on request.
Current Study
To summarize, here we sought to generally replicate and extend Manuck et al.'s (Reference Manuck, Craig, Flory, Halder and Ferrell2011) findings that two ESR1 SNPs moderated the effect of family-rearing experience on age of menarche. The replication effort could only be regarded as “general” rather than “specific” because a different measure of the rearing environment was used in this inquiry (i.e., maternal sensitivity) than in the prior one (i.e., conflict/cohesion). The effort to extend the previous work involved using prospective rather than retrospective measures of the rearing environment, relying on reports of a first menstruation obtained within a year of its occurrence rather than decades later, accounting for some of the heritability of menarcheal age by controlling for maternal age of menarche, and using new analytic methods affording evaluation of the relative fit of two alternative models of G × E interaction to the data.
Method
Participants
Participants included a subset of 373 White females enrolled in the multisite NICHD SECCYD on whom data were available on age of menarche. The limited number of minority females with data on menarcheal age limited hypothesis testing within the subgroups and were thus excluded from this study. In addition, males were also excluded because prior findings indicate that their pubertal development is not regulated by family experience (Belsky, Steinberg, et al., Reference Belsky, Steinberg, Houts, Friedman, DeHart and Cauffman2007).
Procedures and measures
This report focuses on three measurements: maternal sensitivity observed between 6 and 54 months of age, genotype, and age of menarche, assessed between 9.5 and 15 years of age.
Maternal sensitivity was assessed using a mother–child semistructured interaction with toys when the child was 6, 15, 24, 36, and 54 months of age (see NICHD Early Child Care Research Network, 2005, for details). Children were included in this report if mothers had data on at least four out of the five measurement points. Using 4- or 7-point rating scales, maternal sensitivity was evaluated by measuring positive, nonintrusive, responsive, and supportive care. A composite maternal sensitivity score, for 6, 15, and 24 months, was obtained by using the sum scores of sensitivity to nondistress, intrusiveness (reverse scored), and positive regard. For later measurements at 36 and 54 months, supportive presence, hostility (reverse scored), and respect for autonomy scores were summed to create the maternal sensitivity composite. Internal consistency (Cronbach α) for these composites ranged from 0.70 to 0.84. The scores at each of these time points (i.e., 6 to 54 months) were averaged to create the final maternal sensitivity composite score. Analyses were based on a mean-centered version of maternal sensitivity, with M = 0.0, SD = 1.30, and a range from −4.35 to 2.55.
Age of menarche was assessed by asking the girls annually between the ages of 9.5 and 15 years whether they had begun menstruating and, if so, their age at their first menstrual period (in years and months). Mothers were also asked to report on their daughter's first menstrual period, and this data was used if information from the girls was missing. In addition, mothers reported on their own age of menarche, in years and months, which was used to create the dependent variable: a residual score of girl's age at menarche when controlling for maternal age of menarche. Retrospective reports of age of menarche have been shown to be highly reliable even over long time intervals (Must et al., Reference Must, Phillips, Naumova, Blum, Harris and Dawson-Hughes2002); however, this should not be read to imply that such retrospective reports are as accurate as those obtained within a year or less of first menses.
Genotyping and DNA reliability
When the children were 15 years of age, they provided a DNA sample using buccal cheek cells. Genotyping was performed for ESR1 rs9340799 and rs2234693. DNA extraction and genotyping for the SECCYD was performed under the direction of Deborah S. Grove, Director for Genetic Analysis, at the Genome Core Facility in the Huck Institutes for Life Sciences, Penn State University. For this study's subsample (n = 373), the frequency distribution of the first genotype, ESR1 rs9340799, differed significantly from Hardy–Weinberg equilibrium (χ² = 6.1, p < .05); however, the second genotype, ESR1 rs2234693, was not significantly different. Despite the first observation, we chose to analyze the first polymorphism given that we were seeking to replicate the Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) G × E result; appreciating that some might question proceeding in this manner, we return to this issue in the discussion.
Of the full sample of male and female participants who provided DNA (N = 695), 93.7% were able to be genotyped for ESR1 alleles (n = 651). Although we do not know for certain, we suspect that the reason some samples could not be genotyped was because of degradation of the biological material. To be noted in this regard is that the assaying of the two SNPs central to this study took place well after most other genotyping had taken place. This was because until the Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) results appeared, there was no basis (in our minds) for focusing on the two ESR1 SNPs. Reliability for each genotype was tested by analyzing the samples twice (11% of the total n = 651; ESR1 rs9340799, n = 72, ESR1 rs2234693, n = 72), with all discrepancies resolved via a third genotyping. For ESR1 rs9340799 (AA = 0, AG = 1, GG = 2), 16.6% of available samples could not be genotyped in this subsample and κ = 0.99, p < .001, 99% agreement. For ESR1 rs2234693 (TT = 0, CT = 1, CC = 2), 4.1% of the available samples could not be genotyped in this subsample and κ = 0.99, p < .001, 99% agreement. For the two genetic variants, the coding on each SNP (0, 1, or 2) reflects the purported number of risk/plasticity alleles, presuming the G and C alleles are the risk/plasticity alleles for the ESR1 rs9340799 and ESR1 rs2234693 SNPs, respectively. For current analyses, individuals with complete data on all key variables (i.e., mother's age of menarche, maternal sensitivity, and girl's age of menarche) and the ESR1 rs9340799 SNP, the distribution of AA, AG, and GG genotypes was 95, 80, and 35, respectively; and the frequencies of TT, TC, and CC genotypes of ESR1 rs2234693 were 71, 110, and 49, respectively.
Data analysis
Because nonlinear predictor–outcome associations can compromise interpretation of G × E findings, we first evaluated whether the environmental predictor, maternal sensitivity, was nonlinearly related to the outcome variable within ESR1 gene allele groups. Nonlinear functions were nonsignificant in all instances, supporting use of linear models throughout.
Next we fit reparameterized regression models adapted from Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012), which had the form
$$Y = B_0 + B_1 X_1 + B_2\lpar X_2 - C\rpar + B_4\lpar \lpar X_2 - C\rpar \times X_3\rpar + E\comma$$where Y is the dependent variable of girl's age of menarche, X 1 is mother's (residualized) age of menarche (mean centered), X 2 is mother sensitivity (mean centered), X 3 is number of risk/plasticity alleles (0, 1, or 2), C is the crossover point, Y is the intercept or predicted value of B 0 at the crossover point, B 1 is the regression weight for mother's (adjusted) age of menarche, B 2is the slope for mother's sensitivity for persons with 0 risk/plasticity alleles, and B 3 is the shift in the slope for mother's sensitivity for persons with 1 or 2 risk/plasticity alleles. The C parameter estimates the value for maternal sensitivity at which regression lines cross for groups with different numbers of risk/plasticity alleles (0, 1, or 2 risk/plasticity alleles). We also evaluated the fit of additive-, dominant-, and recessive-gene models by altering the scoring of the gene variable X 3. In the additive-gene model, X 3 had its original coding of 0, 1, or 2; in the dominant-gene model, X 3 was rescored so that 0 = 0 risk/plasticity alleles, and 1 = 1 or 2 risk/plasticity alleles; and in the recessive-gene model, X 3 was rescored so that 0 = 0 or 1 risk/plasticity alleles, and 1 = 2 resilience alleles.
The model in Equation 1 is the weak differential-susceptibility model, Model 1w. If B 2 in Equation 1 is fixed at 0, this constrains the environmental effect to be exactly 0 in the least plastic group, leading to the strong differential-susceptibility model, Model 1s. Alternatively, if the crossover point C is fixed at the most positive value for the environment, the model is consistent with weak diathesis–stress thinking, or Model 2w. Finally, if both constraints are invoked (B 2 is fixed at 0 and C is fixed at the most positive value for the environment), the model embodies strong diathesis–stress hypotheses, or Model 2s. For additional details, see Widaman et al. (Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012) and Belsky et al. (Reference Belsky, Pluess and Widaman2013). Recall that we restrict reporting to only results of strong models, because these always proved a better fit than weak ones.
We used both SAS PROC NLIN and PROC NLMIXED to ensure accuracy of results. For all models, PROC NLIN and PROC NLMIXED yielded identical findings. In addition, PROC NLMIXED supplied the Akaike and Bayesian information criteria (AIC and BIC, respectively) that are useful in model comparisons. Lower values of AIC and BIC indicate better fit to the data. Both AIC and BIC contain penalties for model complexity, so adding unnecessary parameters will lead to a rise in the index indicating poorer fit to the data. We evaluated relative model fit using the AIC and BIC in connection with statistical significance of model parameters.
Results
Before conducting two sets of analyses examining, in turn, the moderating effect of each ESR1 SNP under investigation on the relation between maternal sensitivity and menarcheal age, each SNP was correlated with maternal sensitivity to check whether any G × E related finding could be an artifact of a gene to environment correlation. This proved not to be the case, because genotype and maternal sensitivity proved to be unrelated.
Predicting age of menarche using the ESR1 rs9340799 SNP and maternal sensitivity
We first fit linear models predicting age of menarche from maternal sensitivity separately for the three allelic groups (i.e., GG, AG, AA). The slopes for maternal sensitivity for the three groups were as follows: AA group, B = −0.05 (SE = 0.07), t (93) = −0.72, p = .47; AG group, B = 0.06 (SE = 0.10), t (78) = 0.65, p = .52; and GG group, B = 0.29 (SE = 0.17), t (33) = 1.70, p = .09. Inspection of Figure 1a clearly shows that the slope for the GG group (with two risk/plasticity alleles) was positive and much stronger than for the AG and AA groups, thereby providing initial indication of genetic moderation of the relation between maternal sensitivity and age of menarche.
Figure 1. Simple slopes predicting age of menarche from maternal sensitivity for girls in different allelic groups: (a) allelic groups (GG, AG, and AA) based on the ESR1 rs9340799 gene and (b) allelic groups (CC, CT, and TT) based on the ESR1 rs2234693 gene.
Next, we fit several sets of reparameterized models (Equation 1). Comparison of additive-, dominant-, and recessive-gene models revealed that the recessive-gene model provided the closest fit, thus replicating Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011). Therefore, Table 1 displays only results of the recessive-gene model. The strong differential-susceptibility model, in which the B 2 coefficient was fixed at zero, had moderately strong fit to the data (R 2 = .201, p < .001). The coefficient for the G × E interaction was of borderline significance, B 3 = 0.25 (SE = 0.13), p = .055, although in the predicted direction. Further, the crossover point estimate was less than 1 SD above the mean of maternal sensitivity, C = 1.40 (SE = 0.97), 95% confidence interval (CI) (−0.51, 3.31), and so was within the range of the environmental predictor, consistent with differential susceptibility. Diathesis–stress could not be rejected, however, because the 95% CI for the crossover point included the most positive value of maternal sensitivity.
Table 1. Results for alternate regression models predicting girl's age of menarche using the ESR1 rs9340799 gene
Note: N = 210. The values are parameter estimates (standard errors). B 0, The girl's age of menarche for the less plastic group; B 1, the linear effect of mother's age of menarche (centered); B 2, the linear effect of the environmental variable (maternal sensitivity, centered) for the less plastic group; B 3, the linear effect of the environmental variable for the more plastic group; C, the crossover point; AIC, Akaike information criterion; BIC, Bayesian information criterion.
a The parameter is fixed at the reported value; the standard error SE is not applicable, so it is listed as (—).
Table 2. Results for alternate regression models predicting girl's age of menarche using the ESR1 rs2234693 gene
Note: N = 230. The values are parameter estimates (standard errors). B 0, The girl's age of menarche for the less plastic group; B 1, the linear effect of mother's age of menarche (centered); B 2, the linear effect of the environmental variable (maternal sensitivity, centered) for the less plastic group; B 3, the linear effect of the environmental variable for the more plastic group; C, the crossover point; AIC, Akaike information criterion; BIC, Bayesian information criterion.
a The parameter is fixed at the reported value; the standard error is not applicable, so it is listed as (—).
The strong diathesis–stress model shown in Table 1 fixed the crossover point at the most positive value of the environmental predictor (C = 2.55). This model had a small and nonsignificant decrease in fit relative to the strong differential-susceptibility model, ΔR 2 = .002, F (1, 206) = 0.62, ns. AIC and BIC values for this model were smaller, thus better, than comparable values for the differential-susceptibility model, thereby indicating that the diathesis–stress model was a more optimal representation. The coefficient for the G × E interaction was statistically significant, B 3 = 0.16 (SE = 0.06), 95% CI (0.04, 0.28), t (226) = 2.57, p = .01.
The results in Table 1 indicate that the allelic groups had an average age of menarche of 12.55 years when maternal sensitivity was most positive (i.e., maternal sensitivity = 2.55). Then, with each one-unit decrease in maternal sensitivity, the model predicts that menarche will occur 0.16 years, or about 2 months, earlier for girls with the GG genotype. From the most positive environment to the average environment is 2.55 units, leading to the prediction that girls with the GG genotype who were reared in an environment with average levels of maternal sensitivity had average ages of menarche about 5 months earlier than girls raised in the most positive environments, where −2.55 × 0.16 = −0.41 years, or −4.9 months. Maternal sensitivity was unrelated to age of menarche for girls with either the AA or AG genotype, so their predicted age of menarche was 12.55 years, regardless of level of maternal sensitivity.
Predicting age of menarche using the ESR1 rs2234693 SNP
Turning to the second SNP, ESR1 rs2234693, fitting linear models predicting age of menarche from maternal sensitivity separately for the three allelic groups (i.e., CC, CT, and TT) resulted in slopes shown in Figure 1b: TT group, B = −0.07 (SE = 0.09), t (69) = −0.85, p = .40; TC group, B = 0.07 (SE = 0.07), t (108) = 0.93, p = .36; and CC group, B = 0.11 (SE = 0.15), t (47) = 0.72, p = .48. These slopes were more variable than those in Figure 1a, indicating that genetic moderation may be less strong.
Once again, the recessive-gene model provided a closer fit to the data than did additive- or dominant-gene models, so only the recessive-gene models are presented. The strong differential-susceptibility model had moderately strong fit to the data (R 2 = .192, p < .001). The coefficient representing the G × E interaction was not significant, B 3 = 0.08 (SE = 0.11), ns, but in the expected direction. Further, the crossover point estimate, C = 5.18 (SE = 8.39), fell well outside the range of the maternal-sensitivity predictor and was poorly identified (i.e., with very large SE), inconsistent with differential susceptibility. The diathesis–stress model had a small and nonsignificant decrease in model fit relative to the differential-susceptibility model, ΔR 2 = .0012, F (1, 226) = 0.34, ns. In addition, its AIC and BIC values were smaller, thus better, than those for the alternative model, thereby indicating that the diathesis–stress model was a more optimal representation of the data. The coefficient for the G × E interaction proved significant, B 3 = 0.14 (SE = 0.05), 95% CI (0.03, 0.24), t (226) = 2.65, p = .007.
The allelic groups had an average age of menarche of 12.55 years when maternal sensitivity was most positive (i.e., maternal sensitivity = 2.55). Then, with each one-unit decrease in maternal sensitivity, the model predicts that menarche occurred 0.14 years, or about 1.7 months, earlier for girls with the CC genotype. From the most positive environment to the average environment is 2.55 units, leading to the prediction that girls with the CC genotype who were reared in an environment with average levels of maternal sensitivity would have average ages of menarche about 4 months earlier than girls raised in the most positive environments, where −2.55 × 0.14 = −0.357 years, or −4.3 months. Maternal sensitivity was unrelated to age of menarche for girls with either the TT or the TC genotype, so their predicted age of menarche was 12.55 years, regardless of level of maternal sensitivity.
Discussion
Psychosocial acceleration theory stipulates that female pubertal timing is regulated, in part, by experiences in the family while growing up (Belsky, Reference Belsky2012; Belsky et al., Reference Belsky, Steinberg and Draper1991). A differential-susceptibility related revision of the theory postulates, however, that such contextual regulation will vary across individuals, with some proving more susceptible to rearing experiences than others (Belsky, Reference Belsky, Rodgers, Rowe and Miller2000, Reference Belsky2012). Empirical support has emerged for both of these theoretical propositions (Ellis, Shirtcliff, et al., Reference Ellis, Shirtcliff, Boyce, Deardorff and Essex2011; Manuck et al., Reference Manuck, Craig, Flory, Halder and Ferrell2011). Here we sought to replicate and extend Manuck et al.'s (Reference Manuck, Craig, Flory, Halder and Ferrell2011) finding that each of two ESR1 SNPs moderate the effect of family-rearing history on age of menarche. Recall that the replication effort could only be regarded as general because a different measure of the rearing environment was used in this inquiry (i.e., maternal sensitivity) than in the prior one (i.e., family cohesion). The effort to extend the previous work involved using prospective rather than retrospective measures of the rearing environment, relying on reports of a first menstruation made within a year of its occurrence rather than decades later, accounting for some of the heritability of menarcheal age by controlling for maternal age of menarche, and using new analytic methods in an attempt to evaluate the relative fit of two alternative models of G × E interaction to the data. Results revealed that both of the ESR1 polymorphisms that Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) found to moderate the effect of the family rearing environment on age of menarche also moderated the effect of maternal sensitivity on age of menarche in the current inquiry; and, just as in Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011), it was the recessive parameterization of the ESR1 alleles that fit the data best, with exactly the same allelic subgroups proving most susceptible to the environmental regulation of pubertal development.
The new method of analysis afforded competitive evaluation of alternative G × E models, specifically enabling us to represent and test key differences in predictions under different models consistent with different forms of G × E interaction. In analyses using the ESR1 rs9340799 SNP, we could not confidently determine whether evidence of genetic moderation took the form of differential susceptibility or diathesis–stress. Recall that evidence of strong differential susceptibility emerged, because the estimate of the crossover point fell within the range of maternal sensitivity in the sample and the less plastic allelic subgroup had a slope for the relation between maternal sensitivity and age of menarche that did not differ from zero. In contrast, this finding could not be confidently embraced because the 95% confidence interval of the crossover point was above the highest observed value of maternal sensitivity, thereby not allowing us to rule out diathesis–stress. Moreover, when information indices (AIC and BIC) were considered, they provided slightly superior support for the fit of the strong diathesis–stress model. It is important to appreciate with regard to this latter point that these information indices penalize for additional parameter estimates; because the strong diathesis–stress model with only three estimates has fewer estimates than the strong differential-susceptibility model, it is regarded as the more parsimonious model, all other things being equal. As such, we concluded that the strong diathesis–stress model was optimal for the data used in this report for the first ESR1 SNP, rs9340799.
For the second ESR1 SNP, rs2234693, the evidence was not as equivocal, and the strong diathesis–stress predictions were more confidently confirmed. Recall that when the differential-susceptibility model was fit to the data, the point estimate of the crossover point fell outside the range of the environmental variable and was poorly defined (i.e., with large standard error), so differential-susceptibility predictions clearly did not hold. Conversely, the diathesis–stress model had levels of explained variance that were virtually identical to those for the differential-susceptibility model; because the diathesis–stress model had fewer parameter estimates, its AIC and BIC values proved superior. Once again, the strong diathesis–stress model had the best fit to the data based on both AIC and BIC, and the practical magnitude of the effect was not trivial and was quite comparable to that found for the ESR1 rs9340799 SNP.
Because our competitive model testing approach did not find clear, unequivocal evidence favoring differential susceptibility or diathesis–stress in our analyses using the first ESR1 SNP, rs9340799, some might be frustrated when faced with a lack of definitive determination regarding which model fit better. However, we believe our model testing approach has merit specifically because it directly addresses competing hypotheses regarding the form that an interaction takes. If the comparative model-fitting approach had not been adopted, an investigator who used more traditional statistical approaches might draw either a differential-susceptibility or a diathesis–stress related conclusion, perhaps depending on his or her theoretical bent, or might even have failed to uncover empirical support replicating the Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) results.
We suspect that, had we available a larger sample, confidence intervals would have been narrower, allowing a more compelling determination of which model fit the data better. Another facet of the sample size issue was the relatively small size of the samples of individuals who were homozygous for the risk/plasticity allele on the ESR1 SNPs. When analyses are performed on small sample sizes, one or two outliers can have a major effect on results. We screened our data for outliers in each group and found no evidence of such. Still, we recommend that future researchers seek to replicate our analyses in additional, hopefully larger samples to confirm or disconfirm our findings. Needless to say, the relatively small sample size is probably the major limitation of the current study, although restricting ourselves to the sample of Caucasian girls also means that the results here need to be replicated with other ethnic/racial groups before they are presumed to be fully generalizable. Another important limitation is that in the first SNP analyzed, the sample available for analysis did not meet criteria for Hardy–Weinberg equilibrium; even though this single departure from equilibrium was minor, it must be acknowledged that it could compromise interpretation of results. The major strength of the current work, of course, is its attempt to replicate and extend theoretically anticipated results reported by Manuck et al. (Reference Manuck, Craig, Flory, Halder and Ferrell2011) while overcoming the primary limit of the data that those investigators had available and clearly acknowledged: its retrospective character.
What should be clear, in any event, is that to some extent the theorized effect of rearing on pubertal timing in girls applies to some girls more than others and that ESR1 SNPs appear to play a role in determining which girls are more and which are less affected by the developmental experiences that have been studied with regard to pubertal timing. Whether either of the ESR1 SNPs under consideration plays a functional role, or simply is associated with other ESR1 SNPs or even other genes and processes that do so, remains to be determined in future research.
