Introduction

Previous studies in twins and adoptees showed a substantial heritability in internalizing and externalizing problems over childhood [6, 9, 16, 20, 21]. Bartels et al. [6] showed that the heritability of internalizing problems decreased during childhood from 59% at age 3 to 37% at age 12, and the contribution of shared environmental influences increased over time from 13% at age 3 to 37% at age 12. For externalizing problems, small increases in heritability were observed for boys from age 3 (57%) to age 12 (64%), whereas for girls the influences of genetic (50%) and environmental factors remained stable over this period. Significant genetic and shared environmental influences have also been reported for anxious-depressive and withdrawn behavior at age 12 [24].

In addition to differences in heritability as a function of age and sex, genetic effects may also be conditional on environmental exposures, a process often referred to as gene–environment interaction (G × E) [15]. For the current study we were interested in the moderating role of parental divorce on the genetic architecture of internalizing and externalizing problems in children. The experience of parental divorce is often associated with higher mean levels of internalizing and externalizing problems in children [3, 17, 34]. So far, little is known about a possible interaction effect of parental divorce on sources of variation in internalizing and externalizing problems.

In disadvantageous environments, heritabilities may differ from those in advantageous environments. So far, there is empirical support for both higher and lower heritabilities of child problem behavior in the context of disadvantageous environments. Lower heritabilities in disadvantageous environments can result from the suppression of genetic effects by environmental exposures that push or predispose a child to show behavioral problems [10, 32, 35]. Several studies reported that heritabilities are lower in disadvantageous environments. For instance, Hicks et al. [18] showed that the heritability of internalizing problems in 17-year-olds was lower in the context of greater environmental adversity, which included parental divorce. Also, it has been reported that the heritability of adolescent depressive symptoms was lower at higher levels of parental negativity (e.g., frequency of conflicts, punitiveness, and parent–child disagreement) [13]. Another study showed that the heritability of conduct problems in children aged 5–17 years dropped from 77 to 0% at high levels of family dysfunction [11]. The heritability of boys’ antisocial behavior was lower in adolescents from families with lower socioeconomic status [37]. A lower heritability does, however, not necessarily imply less genetic variance. Heritability is a ratio (genetic variance divided by total phenotypic variance), so enlarged environmental variance may result in a lower heritability, even when the absolute amount of genetic variance remains the same. For instance, Feinberg et al. [13] showed that the amount of non-shared environmental variance was higher at high levels of parental negativity, whereas the absolute amount of genetic variance remained stable.

Alternatively, when genetic effects are triggered or amplified by disadvantageous environments, heritabilities will increase [27, 35]. South and Krueger [36] reported a higher heritability of internalizing problems when marital quality was lower, suggesting that persons with a genetic predisposition to internalizing problems may be more likely to express this predisposition in the context of a problematic marriage. Similarly, the heritability of child and adolescent depressive symptoms was higher when the levels of family conflict were high [33]. The heritabilities of adolescent depressive symptoms and externalizing problems were higher when adolescents had experienced more stressful life events [19, 25], and Feinberg et al. [13] reported that the heritability of antisocial behavior was greater at low levels of parental warmth. Regarding these previous studies that report higher as well as lower heritabilities in the context of disadvantageous environments, it is obvious that G × E findings with regard to child and adolescent psychopathology are mixed.

We previously reported in this journal that parental divorce was preceded and predicted by higher levels of externalizing problems in girls at age 3 [34]. Furthermore, parental reports indicated that children at age 12 years with divorced parents showed more internalizing and externalizing problems than children with married parents. The aim of the current study was to investigate as to whether genetic and environmental influences on internalizing and externalizing problems are modified by parental divorce. Therefore, we compared the genetic architecture of internalizing and externalizing problems at ages 3 and 12 between children from divorced and non-divorced families. We expected to find the genetic architecture of pre- and post-divorce internalizing and externalizing problems to be influenced by parental divorce. However, since earlier G × E findings are mixed we could not formulate a specific hypothesis about the direction of effects.

A unique aspect of this study is that all families were still intact at age 3, and a subgroup of children experienced parental divorce between ages 3 and 12 years. This way, we were able to investigate variation in pre-divorce and post-divorce internalizing and externalizing problems in children. Another unique aspect is the large sample size of 4,592 twin pairs. As the genetic architecture of child problem behaviors may be different for boys and girls [6], and because the moderating effect of parental divorce may depend on sex [38], we took divorce, sex, age, and the interaction of sex by divorce into account as moderators.

Methods

Sample and procedure

All participating families are members of the Netherlands Twin Registry (NTR) [7, 8]. For the present study, data from surveys mailed and collected within 3 months of the third and 12th birthdays of twins born between 1986 and 1996 were analyzed. In general, reminders were sent to non-responders 2–3 months after the first mailing. The response rate was 73% at age 3 and 61% at age 12 [7]. For 5,017 twin pairs data were available at age 12 (i.e., survey returned by at least one parent). From this dataset, twin pairs were excluded if one or both of the children suffered from a severe handicap (N = 175), if information on divorce or date of divorce could not be determined (N = 25 twin pairs), if there was a maternal or paternal death (N = 66 twin pairs), if there was an irregular family situation (e.g., the twins did not live with either of their biological parents; N = 39 twin pairs), if parents were already divorced when the child was 3 years old (N = 116 twin pairs), and when information on zygosity was missing (N = 4 twin pairs).

The study sample consisted of 4,592 twin pairs (52% girls), of whom 367 (8%) had experienced a parental divorce between the ages 3 and 12. There were 1,722 monozygotic (MZ) twin pairs, and 2,870 dizygotic (DZ) twin pairs, from which there were 1,386 same-sex DZ twin pairs. For 3,895 twin pairs mother ratings on internalizing and externalizing problems were available at both target ages. For 52 twin pairs, mother ratings were missing at age 12, but information on divorce status was available. For 645 twin pairs, mother ratings on internalizing and externalizing problems were missing at age 3. Twin pairs were not excluded from the study when data on internalizing or externalizing problems were missing on one of the target ages. Attrition analyses showed that drop-out after age 3 was not related to problem behaviors at age 3 [7].

Zygosity was determined for 852 (27%) same-sex twin pairs by DNA analyses or blood group polymorphisms. For all other same-sex twin pairs, zygosity was determined by discriminant analyses using longitudinal questionnaire items on physical resemblance.

Measures

Divorce status

Divorce status for birth cohorts 1986–1991 was assessed in surveys mailed when the children were aged 12 years. Divorce status for birth cohorts 1992 and 1993 was assessed when the children were aged 10 and 12 years. For birth cohorts 1994–1996 we assessed divorce status at ages 7, 10, and 12 years. There were two items in the surveys to assess divorce status. The first item was Which of the following statements corresponds best with the twin’s current family situation? Possible answers were (1) the biological mother and the biological father are married or are living together, (2) the biological mother and the biological father have been divorced since …, (3) the biological mother has been deceased since…, (4) the biological father has been deceased since…, and (5) different family situation, namely… The second item was Which of the following statements corresponds best with the twin’s current living situation? Possible answers were (1) the twins live with both biological parents, (2) the twins live with the biological mother only, (3) the twins live with the biological father only, (4) the twins live with the biological mother and her partner (not the biological father), (5) the twins live with the biological father and his partner (not the biological mother), and (6) different living situation, namely… The information from both items was combined to classify divorce status at ages 3 and 12 years. When parents indicated to be divorced, we deduced divorce status at age 3 from date of divorce.

Internalizing and externalizing problems

At age 3, parental ratings were collected with the Dutch translation of the Child Behavior Checklist (CBCL/2–3) [1], which is a reliable and validated instrument to measure internalizing and externalizing problems in young children [22, 23]. It consists of 100 items that are scored on a 3-point scale that is based on the occurrence of the behavior during the 2 months prior to the survey: 0 if the problem item was not true, 1 if the item was somewhat or sometimes true, and 2 if it was very true or often true. The Internalizing scale (INT3) covers anxious, withdrawn, and depressed behaviors, and consisted of 19 items; the Externalizing scale (EXT3) measures oppositional, aggressive, and overactive behaviors, and consisted of 31 items.

Parental ratings at age 12 were collected with the Dutch translation of the Child Behavior Checklist (CBCL/4–18) [2, 39], that consists of 120 items with the same three response categories as the CBCL/2–3. Parents were instructed to rate the child’s behavior over the preceding 6 months. The Internalizing scale (INT12) measures anxious, depressed, and withdrawn behaviors, as well as functional somatic symptoms, and consists of 31 items. The Externalizing scale (EXT12) consists of 30 items and measures aggressive and rule-breaking behaviors. Sum-scores were used in the analyses.

The good reliability and validity of the CBCL were confirmed in our sample: Cronbach’s alphas were 0.79 (INT3), 0.92 (EXT3), 0.87 (INT12), and 0.90 (EXT12). The reliability estimates were not different between boys and girls. For the present study, only mother ratings were used. Missing data were not replaced by father ratings.

Strategy of analysis

The models for G × E analyses are based on an extension of the classical twin design, in which the relative contributions of genetic and environmental factors to individual differences in internalizing and externalizing problems can be inferred from the different levels of genetic relatedness between MZ and DZ twins. Monozygotic twins are genetically (nearly) identical, while DZ twins share on average 50% of their segregating genes. We can use this difference in similarity to estimate different portions of variance. Individual differences may be due to additive genetic (A), shared environmental (C), and non-shared environmental (E) factors [29]. Shared environmental factors are responsible for resemblance between members of a twin pair, whereas non-shared environmental factors contribute to differences between members of a twin pair. Additive genetic factors are correlated 1.0 in MZ twins, because it is assumed that MZ twins share (nearly) 100% of their genes. On average, additive genetic factors are correlated 0.5 in the DZ twins. By definition, shared environmental factors are correlated 1.0 in both MZ and DZ twins, assuming that the environment shared by a twin pair does not depend on zygosity. This is referred to as the equal environments assumption. Non-shared (or unique) environmental factors are uncorrelated, and may also include measurement error.

Genetic structural equation modeling in Mx [28] was used with the raw data maximum likelihood procedure for estimation of parameters. First, means, variances, and twin correlations for INT3, INT12, EXT3, and EXT12 were estimated. Sex, divorce status, and their interaction, were coded as dummy variables (i.e., 0 and 1), and included as fixed effects in the models. To assess whether they significantly influenced means and variances of INT3, INT12, EXT3, and EXT12, we tested whether constraining each regression weight (i.e., βD, βS, βS×D) at zero led to a significant deterioration of model fit. By including the mean effects in the genetic models, gene–environment correlations were controlled for [31]. Gene–environment correlations are present when exposure to parental divorce depends on a child’s genotype.

The pattern of twin correlations provides a first indication of the magnitude of genetic and environmental influences. If genes contribute to individual differences, the MZ correlations are higher than the DZ correlations. More specifically, there must be an effect of genes when the MZ correlation is about twice as large as the DZ correlation, since MZ twins share twice as many genes as DZ twins. Shared environmental influences are present if the MZ correlations are less than twice as high as the DZ correlations, because this means that the DZ twins are more alike than we would expect on the basis of their shared genes. There are non-shared environmental effects when the MZ correlation is smaller than 1. An indication for G × E is provided when the difference between the MZ and DZ correlations differs between the divorced and non-divorced groups.

Next, sex, divorce status, and their interaction were included as moderators on the genetic architecture of INT3, INT12, EXT3, and EXT12, according to the method as proposed by Purcell [31]. There were nine moderation effects to be tested for each outcome (Fig. 1): the effect of divorce status (D) on A (αD), on C (γD), and on E (ηD), the effect of sex (S) on A (αS), on C (γS), and on E (ηS), and the effect of sex × divorce status (S × D) on A (αS×D), C (γS×D), and E (ηS×D). In all models the effects of sex, divorce status, and their interaction were included as fixed effects. We assessed the significance of the moderation effects by testing whether fixing them at zero resulted in a significant deterioration of model fit. Each moderation effect was separately tested. Assuming that the latent factors A, C, and E have unit variance, the expectation for the phenotype equals:

Fig. 1
figure 1

Genetic model for pre- and post-divorce internalizing and externalizing problems with moderating effects of sex and divorce status. Effect of divorce status (βD), of sex (βS), and their interaction (βS*D) on internalizing and externalizing problems, the triangle (M) represents the INT/EXT residual score, A additive genetic score, C shared environmental factor, E unique environmental factor, a genetic path coefficient, c shared environmental path coefficient, e unique environmental path coefficient, D divorce status, S t1 sex twin 1, S t2 sex twin 2, α D moderation effect of divorce status on A, γ D moderation effect of divorce status on C, η D moderation effect of divorce status on E, α S moderation effect of sex on A, γ S moderation effect of sex on C, η S moderation effect of sex on E, α S×D moderation effect of sex × divorce status on A, γ S×D moderation effect of sex × divorce status on C, η S×D  = moderation effect of sex × divorce status on E. Additive genetic factors are correlated 1.0 in MZ twins and 0.5 in DZ twins. Shared environmental factors are correlated 1.0 in both MZ and DZ twins

Full size image
$$ \begin{gathered} {\text{P }} = \, ({\text{a }} + \, \alpha_{\text{D}} \times {\text{D }} + \alpha_{\text{S}} \times {\text{S }} + \alpha_{{{\text{S}} \times {\text{D}}}} \times {\text{SD}}) \, + \, ({\text{c }} + \gamma_{\text{D}} \times {\text{D }} + \gamma_{\text{S}} \times {\text{S }} + \gamma_{{{\text{S}} \times {\text{D}}}} \times {\text{SD}}) \, + \, ({\text{e }} + \eta_{\text{D}} \times {\text{D }} + \eta_{\text{S}} \times {\text{S }} + \hfill \\ \eta_{{{\text{S}} \times {\text{FS}}}} \times {\text{SD}}) \hfill \\ \end{gathered} $$

Nested models were compared by χ2 tests, in which the χ2 statistic is computed by subtracting the value of −2LL (log-likelihood) for the full model from that for a reduced model (χ2 = −2LL1 − (−2LL0)). This statistic is χ2 distributed with degrees of freedom (df) equal to the difference in the number of parameters estimated in the two models (Δdf = df 1 − df 0). In addition to the χ2 test statistic, Akaike’s Information Criterion (AIC) was computed. The AIC is obtained as: 2k − 2ln(L), where k is the number of parameters in the model, and L is the maximized value of the likelihood function for the estimated model. Thus, with a smaller number of parameters (k) AIC will be lower, and indicates a better model fit. Since all nested models were compared to the full model, the order in which the regression weights were tested did not influence the results. We calculated the genetic and environmental variance components from the parameter estimates for the path coefficients from the best fitting models, such that, for instance, the unstandardized genetic variance equals a2 in the group that did not experience divorce. The unstandardized genetic variance in the group of children who experienced divorce is equal to (a + αD)2. The standardized variance components were computed by dividing the unstandardized variance components by the total variances. The standardized genetic variance components are referred to as the heritabilities. We defined p < 0.05 as statistically significant in all analyses.

Results

Tests of mean effects and twin correlations

In Table 1, means of internalizing and externalizing problems at ages 3 years (pre-divorce) and 12 years (post-divorce) are presented conditional on sex and divorce status. Mean effects of divorce status were significant at age 12 (INT12 βD = 1.39; EXT12 βD = 1.65), showing that children from divorced families had significantly higher mean scores than children from non-divorced families. Mean effects of sex were significant for all outcomes (INT3 βS = −0.17; INT12 βS = 0.30; EXT3 βS = −1.41; EXT12 βS = −0.83). Boys had higher mean scores for INT3, EXT3, and EXT12, and girls had higher mean scores for INT12. Variances (see Table 2, column Total variance) were larger in children from divorced families (INT3: βD = 3.23, χ2(1) = 21.93, p < 0.05; INT12: βD = 3.97, χ2(1) = 10.06, p < 0.05; EXT12: βD = 4.14, χ2(1) = 4.34, p < 0.05), although for EXT3 this effect was only significant for girls (βS×D = 14.77, χ2(1) = 8.83, p < 0.05). For EXT12, the variance was lower in girls than in boys (βS = −7.11, χ2(1) = 179.25, p < 0.05). Twin correlations for INT3, INT12, EXT3, and EXT12 are presented in Table 3. The pattern of MZ and DZ correlations revealed genetic, shared environmental, and non-shared environmental influences on internalizing and externalizing problems at ages 3 and 12. There was no indication for genetic influences on INT12 for the divorced group as the MZ correlation (r = 0.60) was not significantly greater than the DZ correlation (r = 0.63) (χ² (1) = 0.29, p = 0.59). As the correlational structure of the data was not identical for children from married and divorced families, it was justified to start genetic modelling with divorce status as a moderator on the genetic and environmental path coefficients.

Table 1 Maximum likelihood estimates of means of INT3, INT12, EXT3, and EXT12, conditional on sex and divorce status
Full size table
Table 2 Unstandardized and standardized estimates of A, C, and E for internalizing and externalizing problems at ages 3 and 12, conditional on sex and divorce status
Full size table
Table 3 Twin correlations for INT3, INT12, EXT3, and EXT12, conditional on divorce status
Full size table

Model fitting results

Table 4 presents the model fitting results for internalizing and externalizing problems at ages 3 and 12 years. The parameter estimates of the best fitting models are presented in Table 5. For INT3, the best fitting model showed divorce status significantly modified the non-shared environmental variance (ηD = 0.56). For INT12, the best fitting model included the effects of divorce status on all three variance components (αD = −2.01, γD = 2.34, ηD = 0.72). Also, the genetic and non-shared environmental variance were modified by sex × divorce status (αS×D = 2.27; ηS×D = −0.80). For EXT3, sex modified the amount of shared environmental variance (γS = −1.26), as did sex × divorce status (γS×D = 2.06). For EXT12, the best fitting model showed that sex significantly modified the genetic and non-shared environmental variance (αS = −1.65, ηS = −0.59). The non-shared environmental variance was also modified by divorce status (ηD = 1.05) and by sex × divorce status (ηS×D = −0.72).

Table 4 Model fitting results
Full size table
Table 5 Parameter estimates of the best models with their 95% confidence intervals
Full size table

Differences in estimates of ACE

From the parameter estimates of the best fitting models, the standardized and unstandardized contributions of genetic and environmental influences were calculated, which are presented in Table 2. For INT3, the non-shared environmental variance was larger in children from divorced families than in children from non-divorced families, resulting in a lower heritability for children from divorced families (0.49 vs. 0.57). For INT12, the amount of shared environmental variance was larger in boys and girls from divorced families, resulting in lower heritabilities than in children from non-divorced families. Boys from divorced families also showed less unstandardized genetic variance and more non-shared environmental variance. For EXT3, the amount of shared environmental variance was higher for girls from divorced families than for girls from non-divorced families, resulting in a lower heritability for girls from divorced families (0.31 vs. 0.41). For EXT12, boys, and to a lesser extent girls, from divorced families showed more non-shared environmental variance than boys and girls from non-divorced families. Also, boys showed more genetic variance than girls. This resulted in lower heritabilities for children from divorced families.

Overall, genetic factors explained at least 25% of the total variances, except for INT12, for which the heritability was estimated near zero for boys from divorced families. With regard to INT12 and EXT3 (girls only), shared environmental influences were relatively more important for children from divorced families. Non-shared environmental influences were more important for children from divorced families with respect to INT3 and EXT12.

To summarize, parental divorce was associated with higher levels of psychopathology and larger variances of internalizing and externalizing problems. The differences in variance were environmental in origin, and resulted in lower heritabilities for children from divorced families and for children whose parents would later divorce.

Discussion

The present study contributes to the field of G × E interaction research by investigating whether the genetic architecture of internalizing and externalizing problems differed between children living in intact versus divorced families. Children’s behavioral problems were assessed pre-divorce (at age 3 years) as well as post-divorce (at age 12 years). Although age-specific genetic and environmental influences [6] may have contributed to a somewhat different pattern of interaction effects at ages 3 and 12 years, we showed that at both ages heritabilities were slightly lower for children from divorced families than for children from intact families. These results support the hypothesis that genetic influences are less important in the context of parental divorce [10, 32, 35], which is consistent with findings from several other studies [11, 13, 18, 37]. However, as the absolute amounts of genetic variance were similar between children from divorced and intact families (except for boys’ internalizing problems at age 12), genetic effects were not suppressed by the experience of parental divorce. Heritabilities were lower for children from divorced families because for these children there was more environmental variance. This finding emphasizes the importance of considering both the relative and the absolute contributions of genes and environment (i.e., the standardized and unstandardized estimates of genetic and environmental variances) when studying G × E. A sole focus on standardized estimates (i.e., heritabilities) may provide a limited picture of what is going on, and can even lead to misinterpretation [13].

One of the unique aspects of this study is the assessment of pre-divorce as well as post-divorce problems. Although we did not find evidence for higher mean levels of pre-divorce internalizing problems at age 3 in this study as well as in our earlier report [35], our present results show a larger non-shared environmental variance for pre-divorce internalizing problems. Thus, pre-divorce environmental factors increase heterogeneity in children’s internalizing problems without significantly increasing mean levels of internalizing problems.

Our findings are in agreement with previous studies that showed increasing environmental variances of internalizing problems with age [6, 9], which may likely be explained by mounting environmental influences as a result of children’s transition to school, during which they develop social relations with other children, and must cope with several new demands. Our study shows that shared environmental influences become even more important for children of divorced families, which may be an expression of the impact of divorce on family dynamics, such as family relationships. It may also partly reflect children’s understanding of the family situation or twins’ shared access to possible support in relationships outside the family. More research is needed to define which specific environmental factors are of interest in explaining the effects of parental divorce on variation in children’s problem behaviors.

The significant moderation of the shared and non-shared environment by divorce status indicates a greater heterogeneity in internalizing and externalizing problems among children with divorced parents. This underlines the importance of focusing on etiologies of individual differences in the effects of parental divorce instead of focusing on general mean effects.

The ways in which children cope with stressors associated with parental divorce, such as spending less time with the non-custodial parent, parents beginning new romantic relationships, and the intensity, frequency, and style of parental conflict [11, 21], may be important mediators and moderators of the relation between parental divorce and behavioral problems [12]. One could hypothesize that exposure to parental conflict could enhance shared environmental effects, while coping strategies would be more likely to enhance genetic effects.

There seems to be a stronger G × E effect for internalizing problems than for externalizing problems, since divorce moderates genetic effects on internalizing problems, but not those on externalizing problems. One possible explanation for this difference is that additive genetic factors are the main source of stability in externalizing problems over time, whereas most of the stability in internalizing problems can be explained by shared environmental factors [6]. The high genetic contribution to stability in externalizing problems results from the fact that a subset of genes is active at multiple ages throughout childhood. Also, heritability estimates of externalizing problems are generally more similar at different ages during childhood than heritability estimates of internalizing problems [6].

G × E approaches to human behavior are still at the initial stages of inquiry [27]. Research to the role of the environment in genetic expression may be complicated if the population is heterogeneous in the type of G × E interactions [35]. Some interactions may apply only to subgroups of the population, that differ for instance in their social experiences or that have specific genetic variants. For researchers interested in identifying specific genes associated with externalizing or internalizing problems, our results suggest that greater efficiency in detecting associations with specific genes could be obtained by limiting samples to children who have not experienced parental divorce, because then the relative influence of genetics is the largest.

Besides some powerful aspects of this study, such as the large number of participating families, there are some important limitations. Firstly, the observed rate of parental divorce in our twin sample (i.e., 10%) is not representative for the general population. In a large prospective study of Dutch adolescents, about 20% had experienced parental divorce before the age of 12 years [30]. Non-response could be related to divorce status. Non-responders are not necessarily unwilling to return the questionnaire, but as parental divorce often involves change of address, surveys could have been sent to wrong addresses when the families moved without informing the NTR. Secondly, divorced parents may over-report children’s problems or they may be more sensitive to them. Over-reporting could be due to parental characteristics such as anxiety or depression [14], to parents’ concern about their children’s mental health, or to their own distress associated with the divorce [26]. Part of the increased environmental variance found in children from divorced families could therefore be due to over-reporting. However, earlier studies to the effects of rater bias and unreliability showed that these effects distort the estimates of the shared and non-shared environmental factors to a small degree only. For both the internalizing and the externalizing scale, rater bias accounted for at most 13% of the variance, and measurement error and unreliability accounted for less than 11% of the variance [4, 5]. Thirdly, as younger children may be differently affected by parental divorce than older children and adolescents [17], our results cannot be generalized to children of other ages. The heritability estimates obtained in our moderation analyses control for common variance between the moderators and child behavioral problems. Therefore, our estimates may not be directly comparable to studies that only estimate heritability.

In conclusion, environmental influences become more important in explaining variation in children’s problem behaviors in the context of parental divorce. Replication of these G × E effects is necessary. It would be interesting to extend our study with adult data of these children, in order to investigate if the experience of parental divorce in childhood also affects variability in problem behaviors in adulthood.