Childhood parenting and adolescent internalizing and externalizing symptoms: Moderation by multilocus hypothalamic–pituitary–adrenal axis-related genetic variation

Abstract

Genetic variants that regulate hypothalamic–pituitary–adrenal (HPA) axis function have been demonstrated to moderate the association between parenting and mental health. However, extant research has focused primarily on (i) effects of individual genes or (ii) maternal as opposed to paternal parenting. Using a multilocus genetic profile score (MGPS) approach, the current study is the first to examine the moderation effect of multilocus HPA-axis related genetic variants on the association of both maternal and paternal parenting with adolescent internalizing and externalizing symptoms. In a sample of 772 Chinese Han adolescents (M _age = 16.48 ± 1.40 years; 50.1% girls), a theory-driven MGPS was calculated using six polymorphisms within HPA-axis related genes (CRHR1, NR3C1, NR3C2, FKBP5, COMT, and HT1RA). Results showed that the MGPS interacted with both maternal and paternal parenting in the association with adolescent internalizing symptoms, but not externalizing symptoms. Consistent with the differential susceptibility model, adolescents with high versus low MGPS exhibited not only more internalizing symptoms when exposed to low quality of parenting but also less internalizing symptoms when exposed to high quality of parenting. The current findings highlight the potential value of using a multilocus approach to understanding gene-by-environment interaction (G×E) effects underlying mental health. Within such G×E effects, not only maternal but also paternal parenting should be addressed.

Introduction

It has been well documented that internalizing (e.g., depression, anxiety) and externalizing (e.g., conduct problems, hyperactivity) symptoms in adolescence have long-term negative consequences for social relationships, academic performance, and physical health (Achenbach et al., 2016; Ehrenreich et al., 2016; Thapar et al., 2012). Previous research has found evidence for both genetic and environmental influences on adolescent internalizing and externalizing symptoms (Bakermans-Kranenburg & Van IJzendoorn, 2015). The heritability estimates of internalizing and externalizing symptoms vary from 25% to 51% and from 35% to 77%, respectively (e.g., Gjone & Stevenson, 1997; Scourfield et al., 2004). Both internalizing and externalizing symptoms are influenced by many genes of small effect (Plomin, 2013). For example, the 5-HTTLPR and the MAOA-uVNTR have been repeatedly linked to internalizing and externalizing problems (Levinson, 2006; Veroude et al., 2016). Furthermore, there is evidence suggesting that genetic variants that relate to the dysfunction of the hypothalamic–pituitary–adrenal (HPA) axis increase risk for internalizing and externalizing symptoms (Heim et al., 2009; Le Francois et al., 2008; Mulder et al., 2017; Young et al., 2005; Zhou et al., 2017). Similarly, environmental factors such as parenting have been repeatedly found to be associated with both internalizing and externalizing symptoms (Cao et al., 2018; Ong et al., 2018; Pinquart, 2017a, 2017b). Furthermore, it is increasingly recognized that HPA-axis related genes and parenting rarely act alone but interact with each other to explain individual differences in internalizing and externalizing symptoms (e.g., Bevilacqua et al., 2012; Sulik et al., 2015; Zhang et al., 2016; Zhou et al., 2017; Wang et al., 2018). However, extant research has focused primarily on interaction effects of (i) individual HPA-axis related genes or (ii) maternal parenting as opposed to paternal parenting. Adopting a multilocus genetic profile score (MGPS) approach (Belsky & Beaver, 2011; Pagliaccio et al., 2014; Vrshek-Schallhorn et al., 2015) in a sample of Chinese Han adolescents, this study aimed to examine the interaction effect of multilocus HPA-axis related genetic variants on the associations of both maternal and paternal parenting with adolescent internalizing and externalizing symptoms.

According to the attachment theory, childhood parenting is an influential environmental factor in the development of both internalizing and externalizing symptoms in childhood and adolescence (Cao et al., 2018; Fearon et al., 2016; Ong et al., 2018; Pinquart, 2017a, 2017b). Unresponsive and harsh parenting could foster a relatively stable insecure internal working model (IWM) of children, with a representation of the self as unlovable or incompetent, which might increase risk for adolescent internalizing symptoms (Fearon et al., 2016). Furthermore, this insecure IWM might also impede awareness of others’ mental states and reduce the empathy toward others, leading to more adolescent externalizing symptoms (Fearon et al., 2016). A myriad of longitudinal research has supported that adolescents who have experienced unsupportive or hostile parenting of their mothers or fathers in childhood are more likely to display elevated internalizing and externalizing symptoms (Cao et al., 2018; Ong et al., 2018; Pinquart, 2017a, 2017b). However, not all adolescents who have experienced poor parenting show maladjustment. This might depend on plasticity factors that reside in the individuals, especially candidate genes (i.e., gene-by-environment (G×E) interaction), although many of these G×E studies are underpowered (Bakermans-Kranenburg & Van IJzendoorn, 2015; Belsky & Pluess, 2009).

The HPA-axis is one of the main stress-responsive physiological systems (Guerry & Hastings, 2011). Hyperactivity of the HPA-axis has been consistently found to trigger internalizing symptoms, whereas hypoactivity has been related to heightened externalizing symptoms (Buitelaar, 2013; Guerry & Hastings, 2011; Platje et al., 2013). Several genetic variants that modulate HPA-axis functioning (see Table 1) have been examined for possible associations with adolescent internalizing and externalizing symptoms, yielding mixed findings (e.g., Bevilacqua et al., 2012; Heim et al., 2009; Le Francois et al., 2008; Wang et al., 2018; Young et al., 2005; Zhang et al., 2016; Zhou et al., 2017). For instance, several studies demonstrated that individuals with the T allele of FKBP5 rs1360780 are more susceptible to developing depressive symptoms than non-carriers (Zannas & Binder, 2014). However, a recent meta-analysis failed to support this association (Hernández-Díaza et al., 2019).

Table 1. The coding scheme of HPA-axis-related SNPs and their associations with HPA-axis function and related biological processes

Note: ^aA SNP genotype was coded as 1 if it was previously associated with HPA-axis hyperactivity and related biological processes; a SNP genotype was coded as 0 if it was previously associated with HPA-axis hypoactivity and related biological processes. A multilocus genetic profile score (MGPS) was created by summing the scores across these 6 SNPS, where higher genetic profile scores indicate a larger number of alleles that have previously been associated with heightened HPA-axis stress reactivity.

Beyond main effects, genetic variants that modulate HPA-axis function might operate as differential susceptibility factors moderating associations of environmental factors with adolescent internalzing and externalizing symptoms (Ellis & Boyce, 2008). The theory of biological sensitivity to context (BSC theory) proposes that, as a result of temperamental, physiological and genetic makeup, heightened stress reactivity may reflect increased BSC. This increased biological sensitivity has the potential to elicit negative health effects under conditions of adversity versus positive health effects under conditions of support (Ellis & Boyce, 2008). The BSC theory is consistent with the G×E hypothesis of differential susceptibility but contrasts with the diathesis-stress model (Bakermans-Kranenburg & Van IJzendoorn, 2015; Belsky & Pluess, 2009). Whereas the diathesis-stress model views genetic predisposition as a “risk” factor and emphasizes the negative effects of “risk” alleles in contextual adversity, the differential susceptibility model argues that a genetic predisposition should be considered as a “plasticity” factor, which may predispose carriers not only to suffer from adverse environments but also to benefit from supportive environments. Several recent studies have supported this BSC theory by showing an interaction between HPA-axis related genes and parenting on both adolescent internalizing and externalizing symptoms (e.g., Bevilacqua et al., 2012; Mulder et al., 2017; Sulik et al., 2015; Wang et al., 2018; Zhang et al., 2016; Zhou et al., 2017). Specifically, adolescents with alleles of HPA-axis related genes linked to heightened stress reactivity (e.g., high basal cortisol and cortisol reactivity) showed not only more internalizing or externalizing symptoms when exposed to poor parenting but also less symptoms when exposed to supportive parenting. However, it is worth noting that the majority of these G×E studies have focused primarily on maternal parenting, whereas relatively little attention has been given to paternal parenting (e.g., Mulder et al., 2017; Sulik et al., 2015; Zhang et al., 2016).

Furthermore, despite promising findings relating to the interaction of HPA-axis related genes and parenting in the association with mental health, conflicting findings remain in the literature. For instance, whereas some studies on FKBP5 rs1360780 have demonstrated that the T allele moderated the association of poor parenting with internalizing and externalizing symptoms (Bevilacqua et al., 2012; Mulder et al., 2017; Wang et al., 2018), others observed that this was true for CC homozygotes (de Castro-Catalaa et al., 2017) or observed no interaction at all (Bryushkova et al., 2016). One possible explanation for these conflicting findings is that the vast majority of research has examined G×E interactions using individual genes, which could limit the power to find robust G×E effects (Border et al., 2019; Cao et al., 2019; Duncan & Keller, 2011; Zhang & Belsky, 2020). Indeed, there is a growing consensus that mental health has significant polygenic additive underpinnings, with individual genes exerting very small effects that are difficult to capture (Sullivan et al., 2012). As a result, various polygenetic approaches have been developed and used in recent studies, such as: (i) polygenic scores derived from genome-wide association studies and (ii) MGPSs derived theoretically from combining the effects of multiple candidate genes based on the “biological plausibility” (e.g., Belsky & Beaver, 2011; Keers et al., 2016; for a review, see also Zhang & Belsky, 2020). The MGPS, for example, has shown greater predictive validity and power than individual genes examined in isolation (Belsky & Beaver, 2011; Pagliaccio et al., 2014; Starr & Huang, 2019; Vrshek-Schallhorn et al., 2015). For instance, in a sample of 381 Caucasian mother–child dyads, a MGPS examining the multilocus influence of three HPA-axis related genes (i.e., CRHR1, NR3C1, and FKBP5) was found to moderate the relationship between maternal prenatal stress and offspring depression in early adulthood. However, none of the individual HPA-axis related genes yielded a significant G×E interaction when examined in isolation from each other (McKenna et al., 2021).

The extant literature using the HPA-axis related MGPS has primarily focused on interactive effects with various stress exposures (e.g., chronic and acute stress, interpersonal stress) on depression or the structure and function of related brain regions (e.g., amygdala; Di Iorio et al., 2017; Feurer et al., 2017; McKenna et al., 2021; Pagliaccio et al., 2014; Starr & Huang, 2019). To the best of our knowledge, no study to date has utilized this MGPS approach to examine interaction effects between HPA-axis related genes and parenting underlying both internalizing and externalizing symptoms. Furthermore, previous studies on the HPA-axis related MGPS often included Caucasian samples, whereas fewer studies were conducted in Chinese Han samples. The inclusion of a Chinese Han sample is particularly important in light of evidence suggesting that (i) different ethnicities may differ in allele frequencies, linkage disequilibrium and allele functions (van IJzendoorn et al., 2012) and (ii) different ethnical or cultural customs and values may shape different child-rearing practices and distinct parental roles in the family (Chen-Bouck et al., 2019). For instance, due to the influences of Confucianism, Chinese fathers are often viewed as yi jia zhi zhu (i.e., “master of the family”) and assigned in a more powerful position in the family than mothers (Chuang & Su, 2009).

To address these gaps in the extant literature, the current study aimed to investigate the extent to which the multilocus HPA-axis related genetic variants moderated the associations of maternal and paternal parenting with adolescent internalizing and externalizing symptoms in a sample of Chinese Han adolescents. To this end, we established a theory-driven MGPS by counting the number of alleles showing previously established associations with heightened stress reactivity in six HPA-axis related genes: CRHR1 rs110402, NR3C1 rs41423247, NR3C2 rs5522, FKBP5 rs1360780, COMT rs4680, and HTR1A rs6295 (see Table 1). Higher scores on the MGPS indicated a larger number of alleles that have previously been associated with heightened stress reactivity of the HPA-axis. Perceived parenting as well as internalizing and externalizing symptoms were based on adolescent self-reports. We hypothesized that the HPA-axis related MGPS would interact with maternal and paternal parenting in the association with adolescent internalizing and externalizing symptoms. Based on re-parameterized regression models (Belsky & Widaman, 2018; Stocker et al., 2017; Widaman et al., 2012), we expected that these interactions would be consistent with the differential susceptibility model, where for increasing numbers of alleles that previously associated with heightened stress reactivity, the sensitivity toward maternal and paternal parenting might gradually increase. Moreover, in light of the traditions of Confucianism, we hypothesized that the interactive effect would be particularly evident for fathers.

Methods

Participants and procedure

Participants were drawn from an ongoing study investigating genetics and adolescent psychological adjustment in Shandong Province, China. A total of 789 adolescents (M _age = 16.46 ± 1.39 years; 50.1% girls) were recruited from 27 classes of three technical secondary schools located in the cities of Jinan, Rizhao, and Tai’an. All participants were of Chinese Han ethnicity and had no current diagnosis of major physical, neurological, or pervasive developmental disorder. Given the focus on the role of parenting in childhood, only adolescents who were living together with their parents during childhood were included in this study (N = 772; M _age = 16.48 ± 1.40 years; 50.1% girls). The median monthly family income was between CNY 4,000 (about $ 570) and CNY 5,000 (about $ 710). The majority of parents graduated from high school (mothers: 20.9%; fathers: 31.0%), junior high school (mothers: 47.5%; fathers: 47.5%), or primary school (mothers: 23.2%; fathers: 9.6%), whereas a small proportion of mothers (8.4%) and fathers (11.9%) received an education beyond high school.

Adolescents were recruited by electronic brochures that were distributed in online management groups of their classes. A total of 45.6% adolescents were willing to participate in this study. About one week after recruitment, participants were instructed by research assistants to (i) complete questionnaires on parenting and internalizing and externalizing symptoms in their classroom and (ii) provide buccal cells using the cheek swabs kit (Rellagene, Shanghai, China) for DNA extraction. Adolescents refrained from eating, smoking, drinking, or chewing gum in 30 minutes prior to buccal cells collection. Adolescents received small gifts such as notebooks and pencils for their participation. This study was approved by the Ethics Committee of Shandong University. Informed consent was obtained from adolescents, their parents, and headmasters prior to data collection.

Measures

Parenting

Adolescents retrospectively assessed their perceptions of the rearing behaviors of their mothers and fathers in childhood by the widely used Chinese version of the Egna Minnen av Barndoms Uppfostran for Adolescents (EMBU-A, i.e., “My Memories of Upbringing for Adolescents”; Gerlsma et al., 1991; Yue et al., 1993), using a Likert scale from 1 (not true at all) to 4 (almost always true). Three commonly investigated dimensions were assessed: warmth (19 items, e.g., “My mother/father comforted me when I was sad”), punishment (12 items, e.g., “My mother/father gave me more corporal punishment than I had deserved”) and rejection (8 items, e.g., “My mother/father often told me that she/he did not approve of my behavior at home”). The Cronbach’s αs for these dimensions ranged from .81 to .92. After punishment and rejection were reverse-coded, we calculated the mean scores of each dimension allowing a maximum of 25% missing data. The three dimensions significantly correlated with each other for both mothers and fathers (r ≥ .46, p < .001). Confirmatory factor analysis constraining six correlations between the errors in six pairs of items provided the validity of combining these three dimensions into one score for each parent (for maternal parenting: χ ²(693) = 2056.79, p < .001, RMSEA = .05, SRMR = .06, CFI = .89, TLI = .89; for paternal parenting: χ ²(693) = 1960.34, p < .001, RMSEA = .05, SRMR = .06, CFI = .90, TLI = .89). As such, consistent with previous studies (Cao et al., 2018; Stocker et al., 2017), the three dimensions were standardized and averaged into one parenting score for each parent, with higher scores reflecting higher levels of positive parenting. Data inspection showed 10 and 6 outliers (z-score > 3.29) for maternal and paternal parenting, respectively, which were winsorized (i.e., transformed to match the next highest value; Tabachnick & Fidell, 2012).

HPA-axis related multilocus genetic profile score

In order to extend existing research in the field of HPA-axis related MGPSs (Di Iorio et al., 2017; Feurer et al., 2017; McKenna et al., 2021; Pagliaccio et al., 2014; Starr & Huang, 2019), we included into our MGPS not only core SNPs that directly regulate HPA-axis function (CRHR1, NR3C1, NR3C2, FKBP5) but also peripheral SNPs (COMT, HTR1A) that indirectly regulate HPA-axis function (Feurer et al., 2017; Owens et al., 2016). Given the important roles of COMT rs4680 and HTR1A rs6296 in dopamine and serotonin signaling, respectively, these SNPs have previously been used in MGPSs of dopaminergic or serotonergic risk (Stocker et al., 2017; Vrshek-Schallhorn et al., 2015). As such, these peripheral SNPs might be less specific to an MGPS of HPA-axis related risk. Nevertheless, HTR1A rs6295 was selected as a peripheral SNP based on (i) the interrelations between rs6295, serotonin signaling, and HPA-axis function and (ii) the association between rs6295 and a variety of HPA-axis related outcomes, including neural, biological, and behavioral functioning (see Table 1). Similarly, COMT rs4680 was selected based on its robust association with HPA-axis stress responses (see Table 1). To aid comparisons across studies, we additionally examined multiple MGPSs including only core (vs. peripheral) SNPs (see sensitivity analyses).

Buccal cells of adolescents were used to extract genomic DNA and genotyping using a Sequenom chip-based MALDI-TOF mass spectrometry platform (San Diego, CA, USA) according to standard techniques (Cao et al., 2018; Zhang et al., 2016). No SNPs showed significant deviations from the Hardy−Weinberg equilibrium and the minor allele frequencies for all SNPs were greater than .05 (see Supplementary Table S1). Compared to Caucasian samples (National Genomics Data Center Members and Partners, 2020), our Chinese Han sample had less G alleles of rs110402 (12.3% vs. 48.3%, χ ² (1) = 30.86, p < .001), less Met alleles of rs4680 (28.2% vs. 44.8%; χ ² (1) = 5.56, p < .05), but more G alleles of rs6295 (72.3% vs. 51.6%; χ ² (1) = 12.00, p < .001). For the other three SNPs, no significant gene frequency differences were found between our Chinese Han sample and Caucasian samples (χ ² (1) ≤ 3.03, p ≥ .082, see Supplementary Table S1). The MGPS was aggregated based on the previously established effects of individual SNPs on the HPA-axis regulation, where a higher MGPS indicated a larger number of alleles previously associated with heightened stress reactivity of the HPA-axis (see Table 1; Belsky & Beaver, 2011; Pagliaccio et al., 2014). We permitted up to one missing genotype (≤ 20%) per person (Vrshek-Schallhorn et al., 2015; Table S1). The HPA-axis related MGPS ranged from 0 to 6. As only two participants scored 6, we aggregated 6 and 5 (Belsky & Beaver, 2011). The frequencies for the scores within this MGPS were as follows: 0 (N = 48, 6.2%), 1 (N = 139, 18.0%), 2 (N = 226, 29.3%), 3 (N = 226, 29.3%), 4 (N = 108, 14.0%), and 5–6 (N = 25, 3.2%).

Internalizing and externalizing symptoms

Adolescent internalizing and externalizing symptoms in the past 6 months were assessed with the Chinese version of the widely used Strengths and Difficulties Questionnaire (Goodman, 2001), completed by the adolescents. Each item was rated on a 3-point Likert scale, ranging from 0 = not true to 2 = certainly true. Consistent with previous studies in children and adolescents (Kieling et al., 2013; Quach et al., 2018), internalizing symptoms were assessed using the emotional problems subscale (5 items, e.g., “often unhappy, down-hearted or tearful”), and externalizing symptoms were assessed using the subscales of conduct problems (5 items, e.g., “often fights with other young people or bullies them”) and hyperactivity/inattention problems (5 items, e.g., “easily distracted, concentration wanders”). The respective items were summed to create internalizing (M _int = 2.93 ± 2.33, range = 0–10, Cronbach’s α = .75) and externalizing symptoms scores (M _ext = 6.54 ± 2.91, range = 0–20, Cronbach’s α = .61) allowing a maximum of 25% missing data. Higher scores on the internalizing and externalizing scales reflect higher levels of symptoms.

Covariates

Adolescent sex and age were included as covariates. Following Keller (2014), covariates were included as both main and two-way interaction effects with MGPS or parenting (i.e., sex×MGPS, sex×parenting, age×MGPS, age×parenting), to minimize the possibility that significant findings are incorrectly attributed to interactive effects of G×E when the true underlying effect is the result of interactions with covariates.

Statistical analysis

Main analyses

Our main analyses proceeded in six main steps. First, because the interpretation of an additive MGPS for multiple genes can be confounded if genes have differential effects, we conducted an equal gene model to ensure that there was an equal effect across different SNPs in both main and interaction effects (Stocker et al., 2017; see Supplementary Information section “Equal gene model”). Second, bivariate correlations between all main variables and covariates were tested. Third, linear regression analyses were used to examine both main and interaction effects of MGPS and parenting on adolescent internalizing and externalizing symptoms. Although the interaction effects of maternal and paternal parenting were modeled separately, we simultaneously included their main effects (Li & Lee, 2012; Zhang et al., 2016). Nonsignificant interaction effects for covariates were removed from the model by backward elimination. To reduce multicollinearity between product terms, continuous independent variables were mean-centered prior to analyses.

Fourth, we ran an internal replication analysis by randomly splitting our sample and rerunning our regression analyses (i.e., split-half validation; McKenna et al., 2021; Zhang et al., 2016). Fifth, a simple slope test was used to probe the interaction effects that internally replicated. Sixth, in line with previous research (e.g., Belsky & Widaman, 2018; Stocker et al., 2017), we tested the hypothesis that these internally replicated interactions would be consistent with differential susceptibility by conducting re-parameterized regression models with a “linear×linear interaction” (Widaman et al., 2012). For further details, see Supplementary Information section “Re-parameterized regression model”.

We used the false discovery rate (q < .05) to control for multiple testing (Storey, 2002); 2 models (MGPS×maternal parenting, MGPS×paternal parenting) × 2 outcomes (internalizing symptoms, externalizing symptoms) = 4. In addition, a series of follow-up analyses was conducted to examine the role of individual SNPs contained in the MGPS.

Sensitivity analyses

A series of sensitivity analyses was conducted to test the robostness of our main results. Of particular note, using our existing, sequenced data, we were able to conduct a sensitivity analysis using a close approximation of the MGPS previously established by Di Iorio et al., (2017), including CRHR1 rs110402, FKBP5 rs1360780, and NR3C2 rs5522/rs4635799 haplotype. Our approximation MGPS only differed in that NR3C2 rs5522 but not the NR3C2 rs5522/rs4635799 haplotype was included. Moreover, we retested the G×E interactions using the MGPS based on core genes that directly (vs. indirectly) regulate the HPA-axis function (for details, see Supplementary Materials). There is little evidence supporting three-way interactions between sex, HPA-axis related MGPSs, and environmental factors in predicting mental health symptoms (Feurer et al., 2017; Pagliaccio et al., 2014, 2015). Nevertheless, we examined the three-way interaction of sex×MGPS×parenting in a sensitivity analysis. Finally, interaction effects of maternal and paternal parenting with MGPS were also examined simultaneously.

Power analysis

In previous candidate-gene studies, the R ² of interactive effects between the HPA-axis related MGPS and environmental factors on child and adolescent mental health symptoms and related neuro-endophenotypes ranged from .02 to .08 (Di Iorio et al., 2017; Feurer et al., 2017; Mckenna et al., 2021; Pagliaccio et al., 2014; Starr & Huang, 2019). For our linear regression analysis examining the interactive effect of the HPA-axis related MGPS by parenting (α = .05), the statistical power ranged from 97.5% to 100.0% for the full sample (N = 772) and from 79.1% to 100.0% for the two split-half subsamples (N1 = 386, N2 = 386; G*Power 3.1.9.2).

Results

Preliminary analyses

As shown in Table 2, the constrained, equal gene effect model did not significantly differ from the freely estimated gene effect model (ΔR ² ≤ .008, p ≥ .71). This finding suggests an equal effect across each of the six SNPs regarding their main effects and interactive effects, thereby supporting the additivity across the effects of six SNPs and providing the statistical evidence for the rational for creating a MGPS used in this study. In addition, none of the interactions among any two SNPs (i.e., G1×G2 and G1×G2×parenting) was significant, suggesting that the SNPs were not epistatic to each other. Taken together, these tests of the individual SNPs support our MGPS approach.

Table 2. Model comparisons testing equality of effects of individual SNPs on adolescent internalizing and externalizing symptoms

Note: Adj. R ², the R ² estimate corrected for model complexity or the number of predictors; G×MPQ = interaction between gene and maternal parenting quality; G×PPQ = interaction between gene and paternal parenting quality.

As can be seen in Table 3, the MGPS was unrelated to all study variables. Maternal and paternal parenting quality (intercorrelation: r = .76, p < .001) both negatively correlated with internalizing (r _maternal = −.31, p < .001; r _paternal = −.29, p < .001) and externalizing (r _maternal = −.33, p < .001; r _paternal = −.31, p < .001) symptoms (intercorrelation: r = .54, p < .001). Girls were more likely than boys to experience higher quality of paternal parenting and higher levels of internalizing symptoms. Adolescent age was positively correlated with the quality of maternal and paternal parenting and was negatively correlated with externalizing symptoms.

Table 3. Bivariate correlations between the multilocus genetic profile score, parenting and adolescent internalizing and externalizing symptoms

Note. ^aSex: boys = 0, girls = 1; ^b50.1% girls; MGPS = multilocus genetic profile score; MPQ = maternal parenting quality; PPQ = paternal parenting quality;

***p < .001, **p < .01, *p < .05.

G×E interactions and internalizing symptoms

MGPS×maternal parenting interaction

As shown in Table 4, the main effect of MGPS on adolescent internalizing symptoms was not significant. Maternal parenting negatively associated with adolescent internalizing symptoms. As expected, the MGPS significantly interacted with maternal parenting in its association with adolescent internalizing symptoms (b (SE) = −0.26 (0.08), p = .003, q = .006, R ² = .009). This interaction survived after the internal replication analysis (subsample 1: b (SE) = −0.30 (0.13), p = .018, q = .036, R ² = .012; subsample 2: b (SE) = −0.25 (0.10), p = .016, q = .032, R ² = .010; see Supplementary Table S2).

Table 4. Regression models predicting adolescent internalizing and externalizing symptoms from interactions between the multilocus genetic profile score and parenting

Note. MGPS = multilocus genetic profile score; MPQ = maternal parenting quality; PPQ = paternal parenting quality; ^aNonsignificant interaction effects for covariates (e.g., sex×MGPS) were removed from the models by backward elimination; ^bThe false discovery rate (FDR, q < .05) was used to control for multiple testing; 2 models (MGPS×MPQ and MGPS×PPQ models) × 2 outcomes (internalizing and extertnalizing symptoms) = 4.

A simple slope test showed that for increasing numbers of alleles that previously associated with heightened stress reactivity, the sensitivity toward maternal parenting gradually increased (Figure 1a). Specifically, the slopes of maternal parenting on adolescent internalizing symptoms for MGPS of 0, 1, and 2 were non-significant, MGPS = 0, b (SE) = 0.17 (0.27), p = .529; MGPS = 1, b (SE) = −0.06 (0.22), p = .787; MGPS = 2, b (SE) = −0.30 (0.19), p = .115. The slopes for MGPS of 3, 4, and 5–6 were all significant, ranging from b (SE) = −0.54 (0.19), p = .005 for index value of 3 to b (SE) = −1.01 (0.27), p < .001 for an index value of 5–6.

Figure 1. Interaction of the multilocus genetic profile score with (a) maternal parenting and (b) paternal parenting on adolescent internalizing symptoms.

MGPS×paternal parenting interaction

As shown in Table 4, paternal parenting negatively associated with adolescent internalizing symptoms. The MGPS significantly interacted with paternal parenting in relation to adolescent internalizing symptoms (b (SE) = −0.32 (0.08), p < .001, q < .001, R ² = .016). This interaction was internally replicated (subsample 1: b (SE) = −0.39 (0.12), p = .001, q = .004, R ² = .025; subsample 2: b (SE) = −0.31 (0.11), p = .005, q = .020, R ² = .013; see Supplementary Table S2). A simple slope test showed that for increasing numbers of alleles that previously associated with heightened stress reactivity, the sensitivity toward paternal parenting gradually increased (Figure 1b). The slopes of paternal parenting on internalizing symptoms for MGPS of 0 and 1 were nonsignificant, MGPS = 0, b (SE) = 0.22 (0.27), p = .418; MGPS = 1, b (SE) = −0.10 (0.21), p = .631. The slopes for MGPS from 2 to 5–6 were all significant, ranging from b (SE) = −0.42 (0.18), p = .020 for index value of 2 to b (SE) = −1.38 (0.26), p < .001 for an index value of 5–6.

Main and interaction effects of covariates

As shown in Table 4, girls were more likely than boys to exhibit more internalizing symptoms. Adolescent age was unrelated to adolescent internalizing symptoms. Additionally, adolescent sex significantly moderated the associations of both maternal and paternal parenting with internalizing symptoms, such that maternal and paternal parenting negatively associated with girls’ but not boys’ internalizing symptoms (for girls: p < .001; for boys: p ≥ .129). No other significant interactions between covariates (i.e., sex and age) with MGPS or maternal, paternal parenting were observed for internalizing symptoms.

G×E hypothesis of differential susceptibility versus diathesis stress

As shown in Table 5, for both the significant MGPS×maternal parenting and MGPS×paternal parenting interactions, the re-parameterized regression model showed that the differential susceptibility models fitted the data better than the diathesis-stress models (for maternal parenting: ΔR ² = .009, p = .004; for paternal parenting: ΔR ² = .016, p < .001). Moreover, the point estimates of the crossover points of the interactions on the maternal and paternal parenting axes were −0.25 (95%CI [−0.82 to 0.33]) and −0.17 (95%CI [−0.59, 0.25]), respectively, which were both within the range of observed values (ranges for maternal and paternal parenting were from −2.58 to 1.42 and from −2.62 to 1.34, respectively). Taken together, the results suggest that the MGPS×maternal parenting and MGPS×paternal parenting interactions were consistent with the differential susceptibility model, but not the diathesis stress model. That means, the alleles that related to heightened stress reactivity should be “plasticity” rather than “risk” in nature (Belsky & Pluess, 2009). Individuals with more versus less “plasticity alleles” (i.e., high vs. low MGPS) not only suffered from low quality of parenting but also benefited from high quality of parenting, thus exhibited more and less internalizing symptoms, respectively, functioning as a “for better and for worse” pattern.

Table 5. Reparameterized regression model for differential susceptibility and diathesis-stress hypotheses on internalizing symptoms

Note. MGPS = multilocus genetic profile score; MPQ = maternal parenting quality; PPQ = paternal parenting quality; ^a Y _INT = B ₀+B ₁(X _MPQ-C)+B ₂((X _MPQ-C) • X _MGPS)+B ₃ X _sex +B ₄ X _age+B ₅ X _PPQ +B ₆ X _sex• X _MPQ; ^b Y _INT = B ₀+B ₁(X _PPQ-C)+B ₂((X _PPQ-C) • X _MGPS)+B ₃ X _sex+B ₄ X _age+B ₅ X _MPQ+B ₆ X _sex• X _PPQ; where Y _INT = internalizing symptoms, X _MPQ = maternal parenting quality, X _PPQ = paternal parenting quality, X _MGPS = multilocus genetic profile score, X _age = adolescent age; X _sex = adolescent sex, C = crossover point.

G×E interactions and externalizing symptoms

MGPS×maternal parenting interaction

As shown in Table 4, the main effect of MGPS on adolescent externalizing symptoms was not significant. Maternal parenting negatively associated with adolescent externalizing symptoms. However, the MGPS did not significantly interact with maternal parenting in the association with adolescent externalizing symptoms (b (SE) = −0.12 (0.10), p = .251, q = .251, R ² = .002).

MGPS×paternal parenting interaction

As shown in Table 4, paternal parenting negatively associated with externalizing symptoms. However, the MGPS did not significantly interact with paternal parenting in relation to externalizing symptoms (b (SE) = −0.21 (0.10), p = .046, q = .061, R ² = .005).

Main and interaction effects of covariates

As shown in Table 4, neither the main effects of adolescent sex and age nor their interactions with MGPS and maternal/paternal parenting were significant for adolescent externalizing symptoms.

Follow-up analyses: Tests of individual SNPs

First, when examining the role of individual SNPs, none of the polymorphisms yielded a significant interaction with maternal and paternal parenting on internalizing symptoms after multiple testing correction (R ² ≤ .006, p ≥ .013, q ≥ .312; Table 6). Furthermore, excluding one of the six individual SNPs from the MGPS at a time did not change the significant MGPS×parenting interactions on internalizing symptoms (see Supplementary Table S3). These results suggest that, the MGPS provided greater predictive validity than individual SNPs examined in isolation, and the MGPS results were not driven by any individual SNP.

Table 6. Results for regression models predicting adolescent internalizing and externalizing symptoms from interactions between individual SNPs and parenting

Note. G×MPQ = interaction between gene and maternal parenting quality; G×PPQ = interaction between gene and paternal parenting quality; ^aThe false discovery rate (FDR, q < .05) was used to control for multiple testing; 6 SNPs × 2 models (MGPS×MPQ and MGPS×PPQ models) × 2 outcomes (internalizing and externalizing symptoms) = 24.

Sensitivity analyses

As part of our sensitivity analyses, we highlight here four findings. First, the interactions of MGPS on the associations of both maternal and paternal parenting with adolescent internalizing symptoms remained significant, and their effect sizes remained largely unchanged in sensitivity analyses (i) controlling for externalizing symptoms (Supplementary Table S4), (ii) without correcting manternal and paternal parenting for each other (Supplementary Table S5), (iii) considering the individual parenting dimensions (i.e., maternal/paternal warmth and harsh discipline; Supplementary Table S6), (iv) using the MGPS without collapsing groups with 5 or 6 gene variants into one single group (range = 0–6; Supplementary Table S7), (v) using the linear 0/1/2 coding of the MGPS (Supplementary Table S8), (vi) using a close approximation of a HPA-axis-related MGPS established in previous research (Di Iorio et al., 2017; Supplementary Table S9), and (vii) using the MGPS based on core genes of CRHR1, NR3C1, NR3C2, and FKBP5 that directly regulate the HPA-axis function (for details see Supplementary Information section “Results” and Supplementary Table S10). Second, we observed a nonsignificant interaction effect of MGPS and maternal/paternal parenting for the two subscales of externalizing symptoms (i.e., conduct disorder and hyperactivity; Supplementary Table S11). Third, for both internalizing and externalizing symptoms, the three-way interactions of adolescent sex by MGPS by parenting were nonsignificant (R ² ≤ .001; p ≥ .439, q ≥ .919; Supplementary Table S12). Fourth, there was no significant three-way interaction of MGPS, maternal parenting, and paternal parenting (Supplementary Table S13). However, when interaction effects of maternal and paternal parenting with MGPS were examined simultaneously on internalizing symptoms, only the MGPS by paternal parenting interaction remained significant (Supplementary Table S14).

Discussion

Using a recently developed MGPS approach, this study is the first to examine the moderation effect of multilocus HPA-axis related genetic variants in the association of both maternal and paternal parenting with adolescent internalizing and externalizing symptoms in a sample of Chinese Han adolescents. We highlight here three main findings. First, both maternal and paternal parenting, but not the HPA-axis related MGPS, associated with adolescent internalizing and externalizing symptoms. Second, the HPA-axis related MGPS showed an interaction with both maternal and paternal parenting on adolescent internalizing symptoms. No interaction was observed for externalizing problems. Third, consistent with the differential susceptibility model but not the diathesis-stress model, adolescents with high versus low MGPS exhibited not only more internalizing symptoms when exposed to low quality of parenting but also less internalizing symptoms when exposed to high quality of parenting.

The finding that both maternal and paternal parenting quality negatively associated with internalizing and externalizing symptoms is consistent with a recent series of meta-analyses demonstrating that the effects of parental warmth and control on child internalizing and externalizing symptoms did not vary between mothers and fathers (Pinquart, 2017a, 2017b). Our findings further support these meta-analyses (Pinquart, 2017a, 2017b) in that adolescent sex significantly moderated the associations of maternal and paternal parenting with internalizing symptoms, but not externalizing symptoms. Maternal and paternal parenting negatively associated with internalizing symptoms in girls but not in boys. This might be due to a higher prevalence of internalizing symptoms in girls versus boys, with restricted variance in boys potentially leading to weaker associations with parenting (Pinquart, 2017a).

Our study provided preliminary evidence of an interaction between HPA-axis related MGPS and parenting on adolescent internalizing symptoms, with small effect sizes comparable in magnitude for mothers (R ² = .009) and fathers (R ² = .016; Cohen, 1988). Moreover, there was no evidence for sex differences in these G×E interactions. This finding extends previous research showing that the interaction between a dopaminergic and serotonergic-related MGPS (5-HTTLPR, DRD2, DRD4, and COMT) and parenting on adolescent mental health symptoms (i.e., depression, anxiety and hostility) did not differ significantly across Caucasian mothers and fathers (Stocker et al., 2017). However, when the MGPS×maternal parenting and the MGPS×paternal parenting interactions were included simultaneously in our regression model, only the MGPS×paternal parenting interaction remained significant. This unique contribution of paternal parenting might be explained by cultural factors. Collectivism is a central social pattern in China, emphasizing the interdependence among individuals, as well as family harmony and reputation (Chen-Bouck et al., 2019). Based on the traditions of Confucianism, Chinese fathers are often viewed as the master of the family, who enjoy higher interpersonal power and prestige. Chinese fathers seem to show more active involvement in educating and disciplining children (Ho, 1987), whereas Chinese mothers are more involved in daily activities, such as feeding, dressing, and caring (Ho, 1987). As such, the opinion or influence of fathers might be more explicitly weighted by their children, especially by those with genetic plasticity (Khaleque & Rohner, 2012).

The finding that adolescents with high MGPS (i.e., more alleles that linked to heightened stress reactivity) exhibited not only more internalizing symptoms when exposed to low quality of parenting, but also less internalizing symptoms when exposed to high quality of parenting, was in line with the BSC theory. The BSC theory posits that heightened stress reactivity, as a result of temperamental, physiological, and genetic makeup, could elevate individuals’ biological sensitivity to both supportive and adverse contexts (Ellis & Boyce, 2008). Consistently, previous research found that individuals carrying alleles linked to heightened stress reactivity exhibited not only more hyperactivity of the right amygdala towards negative emotion cues (Di Iorio et al., 2017) but also more sustained attention towards positive emotion cues (Owens et al., 2016). Although the mechanism underlying this G×E interaction remains unclear, epigenetic modifications of HPA-axis related genes in response to the early rearing environment might be one possible explanation. For instance, Klengel et al., (2013) found significantly decreased DNA methylation at intron 7 of the FKBP5 gene within T allele carriers who also suffered from childhood trauma. This epigenetic modification in turn might affect the function of immune cells and brain areas associated with stress regulation (e.g., hippocampus), both of which could increase the risk for internalizing symptoms (Klengel et al., 2013). More epigenetic and neuropsychological research is needed to fully understand how HPA-axis related genetic characteristics are linked to sensitivity to the environment.

Our finding of no significant interaction between parenting and HPA-axis related MGPS on adolescent externalizing symptoms supports previous research suggesting that the link with HPA-axis functioning is weaker for externalizing than for internalizing symptoms (Alink et al., 2010; Buitelaar, 2013; Guerry & Hastings, 2011). However, a promising avenue for future research is the assessment of HPA-axis related MGPS×parenting interactions on externalizing symptoms using clinical samples (Alink et al., 2010) or among adolescents with persistent externalizing symptoms (Platje et al., 2013).

The current study extended previous research in five key ways. First, a theory-driven HPA-axis related MGPS was established, which increased the G×E power by aggregating genetic variants. Using this MGPS, we also tested the two competitive G×E hypotheses (differential susceptibility vs. diathesis-stress). Second, unlike previous research focusing primarily on maternal parenting, this study investigated both maternal and paternal parenting in G×E interactions. Meanwhile, both internalizing and externalizing symptoms were investigated. Third, the inclusion of a sample of Chinese Han adolescents is a strength of this study, as previous research in the field mainly focused on Caucasian samples. Fourth, to ensure the robustness of our findings, we (a) controlled for covariates using both main and interaction effects, (b) internally replicated our findings by randomly splitting our sample and rerunning the regression analyses, and (c) showed that the MGPS results were not fully accounted for by individual genes. Finally, we also demonstrated an equal effect across each of the six SNPs regarding their main effects and interactive effects and found no interactions among any two SNPs (i.e., G1×G2 and G1×G2×parenting), thereby providing further support for the MGPS approach.

Despite these strengths, several limitations should be considered. First, this study used adolescent perceived, retrospective ratings to assess childhood parenting, which may introduce recall bias. Another limitation of this and previous G×E studies in the field (e.g., Surtees et al., 2006) is that measures of perceived parenting and mental health symptoms were both based on adolescent self-reports. For example, adolescents with more internalizing symptoms may be more likely to report poor parenting, rendering the study susceptible to measurement error and reduced power to find significant G×E interactions (Bakermans-Kranenburg & Van Ijzendoorn, 2015). However, it is also possible that the development of adolescent mental health is influenced more by perceived parenting than by parenting that is rated or observed by others (Demo et al., 1987; Ong et al., 2018). Furthermore, the measure we used to assess adolescent perceived parenting, the EMBU-A, has been widely used and validated in Chinese Han adolescents and various other adolescent samples throughout the world (Gerlsma et al., 1991; Yue et al., 1993). Nevertheless, research with multiple reporters is needed to verify our findings. Second, while the MGPS approach increased G×E power relative to an individual-gene approach, it is also faced with challenges. For instance, the MGPS assumes additivity within a given system but ignores potential epistatic interactions (Di Iorio et al., 2017). It is worth noting that the present study showed no low-order epistatic effects between genes and furthermore supported the additivity across genes using the recently developed equal gene model (Stocker et al., 2017). Our MGPS included not only core SNPs that directly regulate HPA-axis function but also peripheral SNPs that indirectly regulate HPA-axis function, with COMT rs4680 and HTR1A rs6295 being more central to the dopamine and serotonin systems, respectively. Although it might be argued that an HPA-axis related MGPS including only core SNPs is more beneficial for understanding etiology (Vrshek-Schallhorn et al., 2015), our sensitivity analysis suggested that an MGPS including both core and peripheral SNPs shows a slightly, albeit not significantly, larger interaction effect than an MGPS including only core SNPs. It has been previously noted that a greater coverage of HPA-axis related genes will likely yield stronger associations (Feurer et al., 2017), and that research might be advanced by moving beyond core SNPs to subsequently incorporate additional SNPs (Vrshek-Schallhorn et al., 2015). Future research is needed to test the replicability of our findings in independent samples and other ethnicities.

In summary, HPA-axis related genes showed an additive polygenic plasticity towards both maternal and paternal parenting in association with adolescent internalizing symptoms, but not externalizing symptoms. The polygenic plasticity of HPA-axis related genes worked in a “for better and for worse” manner that was consistent with the differential susceptibility hypothesis. The current findings highlight the potential value of using a multilocus approach to understanding G×E effects underlying mental health. Meanwhile, this study suggests that, within such G×E effects, not only maternal but also paternal parenting should be addressed.