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., Reference Achenbach, Ivanova, Rescorla, Turner and Althoff2016; Ehrenreich et al., Reference Ehrenreich, Beron and Underwood2016; Thapar et al., Reference Thapar, Collishaw, Pine and Thapar2012). Previous research has found evidence for both genetic and environmental influences on adolescent internalizing and externalizing symptoms (Bakermans-Kranenburg & Van IJzendoorn, Reference Bakermans-Kranenburg and Van Ijzendoorn2015). The heritability estimates of internalizing and externalizing symptoms vary from 25% to 51% and from 35% to 77%, respectively (e.g., Gjone & Stevenson, Reference Gjone and Stevenson1997; Scourfield et al., Reference Scourfield, van den Bree, Martin and McGuffin2004). Both internalizing and externalizing symptoms are influenced by many genes of small effect (Plomin, Reference Plomin2013). For example, the 5-HTTLPR and the MAOA-uVNTR have been repeatedly linked to internalizing and externalizing problems (Levinson, Reference Levinson2006; Veroude et al., Reference Veroude, Zhang-James, Fernàndez-Castillo, Bakker, Cormand and Faraone2016). 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., Reference Heim, Bradley, Mletzko, Deveau, Musselman, Nemeroff, Ressler and Binder2009; Le Francois et al., Reference Le Francois, Czesak, Steubl and Albert2008; Mulder et al., Reference Mulder, Rijlaarsdam, Luijk, Verhulst, Felix, Timeier, Bakermans-Kranenburg and van IJzendoorn2017; Young et al., Reference Young, Abelson and Cameron2005; Zhou et al., Reference Zhou, Fan, Peng, Li, Long, Zhou and Liang2017). Similarly, environmental factors such as parenting have been repeatedly found to be associated with both internalizing and externalizing symptoms (Cao et al., Reference Cao, Rijlaarsdam, Voort, Ji, Zhang and Bakermans-Kranenburg2018; Ong et al., Reference Ong, Eilander, Saw, Xie and Broekman2018; Pinquart, Reference Pinquart2017a, Reference Pinquart2017b). 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., Reference Bevilacqua, Carli, Sarchiapone, George, Goldman, Roy and Enoch2012; Sulik et al., Reference Sulik, Eisenberg, Spinrad, Lemery-Chalfant, Swann, Silva, Reiser, Stover and Verrelli2015; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016; Zhou et al., Reference Zhou, Fan, Peng, Li, Long, Zhou and Liang2017; Wang et al., Reference Wang, Shelton and Dwivedi2018). 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, Reference Belsky and Beaver2011; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg and Craske2015) 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., Reference Cao, Rijlaarsdam, Voort, Ji, Zhang and Bakermans-Kranenburg2018; Fearon et al., Reference Fearon, Groh, Bakermans-Kranenburg, van IJzendoorn, Roisman and Cichetti2016; Ong et al., Reference Ong, Eilander, Saw, Xie and Broekman2018; Pinquart, Reference Pinquart2017a, Reference Pinquart2017b). 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., Reference Fearon, Groh, Bakermans-Kranenburg, van IJzendoorn, Roisman and Cichetti2016). 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., Reference Fearon, Groh, Bakermans-Kranenburg, van IJzendoorn, Roisman and Cichetti2016). 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., Reference Cao, Rijlaarsdam, Voort, Ji, Zhang and Bakermans-Kranenburg2018; Ong et al., Reference Ong, Eilander, Saw, Xie and Broekman2018; Pinquart, Reference Pinquart2017a, Reference Pinquart2017b). 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, Reference Bakermans-Kranenburg and Van Ijzendoorn2015; Belsky & Pluess, Reference Belsky and Pluess2009).
The HPA-axis is one of the main stress-responsive physiological systems (Guerry & Hastings, Reference Guerry and Hastings2011). Hyperactivity of the HPA-axis has been consistently found to trigger internalizing symptoms, whereas hypoactivity has been related to heightened externalizing symptoms (Buitelaar, Reference Buitelaar2013; Guerry & Hastings, Reference Guerry and Hastings2011; Platje et al., Reference Platje, Jansen, Raine, Branje, Doreleijers, De Vries-Bouw, Popma, van Lier, Koot, Meeus and Vermeiren2013). 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., Reference Bevilacqua, Carli, Sarchiapone, George, Goldman, Roy and Enoch2012; Heim et al., Reference Heim, Bradley, Mletzko, Deveau, Musselman, Nemeroff, Ressler and Binder2009; Le Francois et al., Reference Le Francois, Czesak, Steubl and Albert2008; Wang et al., Reference Wang, Shelton and Dwivedi2018; Young et al., Reference Young, Abelson and Cameron2005; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016; Zhou et al., Reference Zhou, Fan, Peng, Li, Long, Zhou and Liang2017). 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, Reference Zannas and Binder2014). However, a recent meta-analysis failed to support this association (Hernández-Díaza et al., Reference Hernández-Díaza, González-Castro, Tovilla-Zárate, Juárez-Rojop, López-Narváez, Pérez-Hernández, Rodríguez-Pérez and Genis-Mendoza2019).
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, Reference Ellis and Boyce2008). 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, Reference Ellis and Boyce2008). The BSC theory is consistent with the G×E hypothesis of differential susceptibility but contrasts with the diathesis-stress model (Bakermans-Kranenburg & Van IJzendoorn, Reference Bakermans-Kranenburg and Van Ijzendoorn2015; Belsky & Pluess, Reference Belsky and Pluess2009). 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., Reference Bevilacqua, Carli, Sarchiapone, George, Goldman, Roy and Enoch2012; Mulder et al., Reference Mulder, Rijlaarsdam, Luijk, Verhulst, Felix, Timeier, Bakermans-Kranenburg and van IJzendoorn2017; Sulik et al., Reference Sulik, Eisenberg, Spinrad, Lemery-Chalfant, Swann, Silva, Reiser, Stover and Verrelli2015; Wang et al., Reference Wang, Shelton and Dwivedi2018; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016; Zhou et al., Reference Zhou, Fan, Peng, Li, Long, Zhou and Liang2017). 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., Reference Mulder, Rijlaarsdam, Luijk, Verhulst, Felix, Timeier, Bakermans-Kranenburg and van IJzendoorn2017; Sulik et al., Reference Sulik, Eisenberg, Spinrad, Lemery-Chalfant, Swann, Silva, Reiser, Stover and Verrelli2015; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016).
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., Reference Bevilacqua, Carli, Sarchiapone, George, Goldman, Roy and Enoch2012; Mulder et al., Reference Mulder, Rijlaarsdam, Luijk, Verhulst, Felix, Timeier, Bakermans-Kranenburg and van IJzendoorn2017; Wang et al., Reference Wang, Shelton and Dwivedi2018), others observed that this was true for CC homozygotes (de Castro-Catalaa et al., Reference de Castro-Catalaa, Peña, Kwapil, Papiol, Sheinbaum, Cristóbal-Narváez, Ballespí, Barrantes-Vidal and Rosa2017) or observed no interaction at all (Bryushkova et al., Reference Bryushkova, Zai, Chen, Pappa, Mileva, Tiemerier, Bakermans-Kranenburg, Kennedy and Beitchman2016). 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., Reference Border, Johnson, Evans, Smolen, Berley, Sullivan and Keller2019; Cao et al., Reference Cao, Cao and Chen2019; Duncan & Keller, Reference Duncan and Keller2011; Zhang & Belsky, Reference Zhang and Belsky2020). 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., Reference Sullivan, Daly and OʼDonovan2012). 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, Reference Belsky and Beaver2011; Keers et al., Reference Keers, Coleman, Lester, Roberts, Breen, Thastum, Bögels, Schneider, Heiervang, Meiser-Stedman, Nauta, Creswell, Thirlwall, Rapee, Hudson, Lewis, Plomin and Eley2016; for a review, see also Zhang & Belsky, Reference Zhang and Belsky2020). The MGPS, for example, has shown greater predictive validity and power than individual genes examined in isolation (Belsky & Beaver, Reference Belsky and Beaver2011; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014; Starr & Huang, Reference Starr and Huang2019; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg and Craske2015). 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., Reference Mckenna, Hammen and Brennan2021).
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., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017; Feurer et al., Reference Feurer, Mcgeary, Knopik, Brick and Gibb2017; McKenna et al., Reference Mckenna, Hammen and Brennan2021; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014; Starr & Huang, Reference Starr and Huang2019). 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., Reference van IJzendoorn, Belsky and Bakermans-Kranenburg2012) 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., Reference Chen-Bouck, Patterson and Chen2019). 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, Reference Chuang and Su2009).
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, Reference Belsky and Widaman2018; Stocker et al., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017; Widaman et al., Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012), 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., Reference Gerlsma, Arrindell, Van der Veen and Emmelkamp1991; Yue et al., Reference Yue, Li, Jin and Ding1993), 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: χ 2(693) = 2056.79, p < .001, RMSEA = .05, SRMR = .06, CFI = .89, TLI = .89; for paternal parenting: χ 2(693) = 1960.34, p < .001, RMSEA = .05, SRMR = .06, CFI = .90, TLI = .89). As such, consistent with previous studies (Cao et al., Reference Cao, Rijlaarsdam, Voort, Ji, Zhang and Bakermans-Kranenburg2018; Stocker et al., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017), 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, Reference Tabachnick and Fidell2012).
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., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017; Feurer et al., Reference Feurer, Mcgeary, Knopik, Brick and Gibb2017; McKenna et al., Reference Mckenna, Hammen and Brennan2021; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014; Starr & Huang, Reference Starr and Huang2019), 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., Reference Feurer, Mcgeary, Knopik, Brick and Gibb2017; Owens et al., Reference Owens, Harrison, Burkhouse, McGeary, Knopik, Palmer and Gibb2016). 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., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017; Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg and Craske2015). 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., Reference Cao, Rijlaarsdam, Voort, Ji, Zhang and Bakermans-Kranenburg2018; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016). 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%, χ 2 (1) = 30.86, p < .001), less Met alleles of rs4680 (28.2% vs. 44.8%; χ 2 (1) = 5.56, p < .05), but more G alleles of rs6295 (72.3% vs. 51.6%; χ 2 (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 (χ 2 (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, Reference Belsky and Beaver2011; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014). We permitted up to one missing genotype (≤ 20%) per person (Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg and Craske2015; 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, Reference Belsky and Beaver2011). 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, Reference Goodman2001), 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., Reference Kieling, Hutz, Genro, Polanczyk, Anselmi, Camey, Hallal, Barros, Victora, Menezes and Rohde2013; Quach et al., Reference Quach, Nguyen, Williams and Sciberras2018), 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 (Reference Keller2014), 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., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017; 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, Reference Li and Lee2012; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016). 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., Reference Mckenna, Hammen and Brennan2021; Zhang et al., Reference Zhang, Cao, Wang, Ji and Cao2016). 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, Reference Belsky and Widaman2018; Stocker et al., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017), 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., Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012). 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, Reference Storey2002); 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., (Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017), 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., Reference Feurer, Mcgeary, Knopik, Brick and Gibb2017; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014, Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2015). 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 2 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., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017; Feurer et al., Reference Feurer, Mcgeary, Knopik, Brick and Gibb2017; Mckenna et al., Reference Mckenna, Hammen and Brennan2021; Pagliaccio et al., Reference Pagliaccio, Luby, Bogdan, Agrawal, Gaffrey, Belden, Botteron, Harms and Barch2014; Starr & Huang, Reference Starr and Huang2019). 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 2 ≤ .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 2, the R 2 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 2 = .009). This interaction survived after the internal replication analysis (subsample 1: b (SE) = −0.30 (0.13), p = .018, q = .036, R 2 = .012; subsample 2: b (SE) = −0.25 (0.10), p = .016, q = .032, R 2 = .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 2 = .016). This interaction was internally replicated (subsample 1: b (SE) = −0.39 (0.12), p = .001, q = .004, R 2 = .025; subsample 2: b (SE) = −0.31 (0.11), p = .005, q = .020, R 2 = .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 2 = .009, p = .004; for paternal parenting: ΔR 2 = .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, Reference Belsky and Pluess2009). 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 0+B 1(X MPQ-C)+B 2((X MPQ-C) • X MGPS)+B 3 X sex +B 4 X age+B 5 X PPQ +B 6 X sex• X MPQ; b Y INT = B 0+B 1(X PPQ-C)+B 2((X PPQ-C) • X MGPS)+B 3 X sex+B 4 X age+B 5 X MPQ+B 6 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 2 = .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 2 = .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 2 ≤ .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., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017; 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 2 ≤ .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, Reference Pinquart2017a, Reference Pinquart2017b). Our findings further support these meta-analyses (Pinquart, Reference Pinquart2017a, Reference Pinquart2017b) 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, Reference Pinquart2017a).
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 2 = .009) and fathers (R 2 = .016; Cohen, Reference Cohen1988). 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., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017). 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., Reference Chen-Bouck, Patterson and Chen2019). 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, Reference Ho and Lamb1987), whereas Chinese mothers are more involved in daily activities, such as feeding, dressing, and caring (Ho, Reference Ho and Lamb1987). As such, the opinion or influence of fathers might be more explicitly weighted by their children, especially by those with genetic plasticity (Khaleque & Rohner, Reference Khaleque and Rohner2012).
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, Reference Ellis and Boyce2008). 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., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017) but also more sustained attention towards positive emotion cues (Owens et al., Reference Owens, Harrison, Burkhouse, McGeary, Knopik, Palmer and Gibb2016). 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., (Reference Klengel, Mehta, Anacker, Rex-Haffner, Pruessner, Pariante, Pace, Mercer, Mayberg, Bradley, Nemeroff, Holsboer, Heim, Ressler, Rein and Binder2013) 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., Reference Klengel, Mehta, Anacker, Rex-Haffner, Pruessner, Pariante, Pace, Mercer, Mayberg, Bradley, Nemeroff, Holsboer, Heim, Ressler, Rein and Binder2013). 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., Reference Alink, Ijzendoorn, Bakermans-Kranenburg, Mesman, Juffer and Koot2010; Buitelaar, Reference Buitelaar2013; Guerry & Hastings, Reference Guerry and Hastings2011). 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., Reference Alink, Ijzendoorn, Bakermans-Kranenburg, Mesman, Juffer and Koot2010) or among adolescents with persistent externalizing symptoms (Platje et al., Reference Platje, Jansen, Raine, Branje, Doreleijers, De Vries-Bouw, Popma, van Lier, Koot, Meeus and Vermeiren2013).
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., Reference Surtees, Wainwright, Willis-Owen, Luben, Day and Flint2006) 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, Reference Bakermans-Kranenburg and Van Ijzendoorn2015). 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., Reference Demo, Small and Savin-Williams1987; Ong et al., Reference Ong, Eilander, Saw, Xie and Broekman2018). 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., Reference Gerlsma, Arrindell, Van der Veen and Emmelkamp1991; Yue et al., Reference Yue, Li, Jin and Ding1993). 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., Reference Di Iorio, Carey, Michalski, Corral-Frias, Conley, Hariri and Bogdan2017). 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., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017). 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., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg and Craske2015), 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., Reference Feurer, Mcgeary, Knopik, Brick and Gibb2017), and that research might be advanced by moving beyond core SNPs to subsequently incorporate additional SNPs (Vrshek-Schallhorn et al., Reference Vrshek-Schallhorn, Stroud, Mineka, Zinbarg and Craske2015). 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0954579421001620.
Acknowledgments
We are extremely grateful to all the families who took part in this study. We gratefully acknowledge the research assistants who contributed to the data generation.
Funding statement
This work was supported by the National Natural Science Foundation of China (31800936), Natural Science Foundation of Shandong Province (ZR2019BC111), Social Science Planning Fund of Shandong Province (18DJYJ01), and Future Plan of Young Scholars in Shandong University (21330089964217) to C.C. The work of J.R. was supported by the Netherlands Organization for Scientific Research (NWO ZonMw VENI, grant no. 91618147).
Conflicts of interest
The authors declare that they have no conflict of interest.
Ethical standards
Informed consent was obtained from all involved adolescents, parents, and school principals. All procedures followed were in accordance with the ethical standards of the ethics committee on human experimentation of Shandong University and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.
