Introduction
Internalizing and externalizing behaviors are common behavior problems among young children that are tightly linked to an increased risk of worse interpersonal relationships, academic underachievement, school dropout, and later psychopathology, especially in adolescence (Basten et al., Reference Basten, Tiemeier, Althoff, van de Schoot, Jaddoe, Hofman, Hudziak, Verhulst and van der Ende2016; McGrath & Elgar, Reference McGrath, Elgar and Wright2015). The crucial role of caregiver parenting in the etiology of child internalizing and externalizing behaviors has been well established (see Pinquart, Reference Pinquart2017a, Reference Pinquart2017b, for meta-analyses), even though not all individuals experiencing poor parenting show behavior problems (Belsky & Pluess, Reference Belsky and Pluess2009; Boyce & Ellis, Reference Boyce and Ellis2005). Recently, epigenetic modifications, especially DNA methylation (DNAm) of a stress-related gene – the nuclear receptor subfamily 3, group C, member 1 (NR3C1) – have been suggested to help explain individual differences in the parenting–behavior problems link (i.e., Epigenome × Environment interaction, Epi × E; Cecil et al., Reference Cecil, Neumann and Walton2023). However, the direct evidence for this Epi × E interaction, especially in the context of Chinese culture, is still lacking (only one study on German adolescents; Radtke et al., Reference Radtke, Schauer, Gunter, Ruf-Leuschner, Sill, Meyer and Elbert2015). Moreover, a series of important issues, including the relevant dynamic effects, function models, and potential moderators (parental and child sex and genotypes), remain unresolved. Therefore, using a one-and-a-half-year longitudinal design and multimodal data (epigenetic, genetic, and multi-informant behavioral data), the current study aimed to systematically examine the interactive effect of NR3C1 methylation × caregiver parenting on child internalizing and externalizing behaviors in Chinese preschoolers.
Over the past decade, a rapidly growing body of literature has underscored the significant role of epigenetic modifications in the pathology of child behavior problems (Cecil et al., Reference Cecil, Neumann and Walton2023). Epigenetic modifications refer to alterations that influence gene expression without altering the DNA sequence. As the most stable epigenetic marker, DNAm involves the addition of a methyl group to individual cytosines in the context of cytosine–phosphate–guanine (CpG) nucleotide sites (Cecil et al., Reference Cecil, Neumann and Walton2023). When this addition occurs in gene promoters, it is often associated with transcriptional gene silencing or reduced expression, leading to physiological, psychological, and behavioral consequences (Cecil et al., Reference Cecil, Neumann and Walton2023). In epigenetic research on mental health, DNAm of the NR3C1 gene has gained considerable attention due to its pivotal regulatory effect on the function of the hypothalamic–pituitary–adrenal (HPA) axis; higher and lower HPA-axis reactivity have been found to trigger child internalizing and externalizing behaviors, respectively (Koss & Gunnar, Reference Koss and Gunnar2018).
To be specific, the NR3C1 gene, also known as the glucocorticoid receptor gene, could modulate the negative feedback of the HPA-axis by encoding glucocorticoid receptors (McGowan et al., Reference McGowan, Sasaki, D’Alessio, Dymov, Labonté, Szyf, Turecki and Meaney2009). Although its methylation levels differ across various regions and are usually low at the promoter region, methylation occurring in the exon 1F region of the promoter (NR3C1-1F), which contains the canonical nerve growth factor-inducible protein A (NGFI-A) transcription factor binding sites (see Supplementary Figure S1), has been found to have a significant effect on child development (see Watkeys et al., Reference Watkeys, Kremerskothen, Quide, Fullerton and Green2018, for a review). In rats and humans, these alterations in NR3C1-1F could downregulate gene expression by physically impeding NGFI-A from binding to the DNA sequence (McGowan et al., Reference McGowan, Sasaki, D’Alessio, Dymov, Labonté, Szyf, Turecki and Meaney2009). In a pioneering study (Weaver et al., Reference Weaver, Cervoni, Champagne, D’Alessio, Sharma, Seckl, Dymov, Szyf and Meaney2004), greater methylation of NR3C1 exon 17 (homologous to human exon 1F) in the hippocampal region was found in rats exposed to poor maternal care, which resulted in reduced gene expression, decreased density of hippocampal glucocorticoid receptors, and worse physiological (e.g., higher HPA-axis reactivity) and behavioral outcomes (e.g., depression-like behaviors). Inspired by this groundbreaking work, subsequent human studies have examined the associations between NR3C1 methylation in various tissues (e.g., postmortem brain tissue, blood, buccal cells, and saliva) and child behavior problems, especially internalizing behaviors, but have yielded conflicting results (see Watkeys et al., Reference Watkeys, Kremerskothen, Quide, Fullerton and Green2018, for a review). Nevertheless, growing research has begun to focus on this emerging field and target NR3C1 methylation in an attempt to unravel the antecedents of child behavior problems.
Beyond the main effects, NR3C1 methylation has been suggested as a promising factor that may explain individual differences in the parenting–behavior problems link (Bosmans et al., Reference Bosmans, Young and Hankin2018; Radtke et al., Reference Radtke, Schauer, Gunter, Ruf-Leuschner, Sill, Meyer and Elbert2015). According to attachment theory (Fearon et al., Reference Fearon, Groh, Bakermans-Kranenburg, van Ijzendoorn, Roisman and D.2016), caregiver parenting, especially in the early years of life, plays an essential role in both child internalizing and externalizing behaviors by affecting the formation of the internal working model, a cognitive framework for how the child understands the world, the self, and others. Substantial cross-sectional and longitudinal evidence supports the links between caregiver parenting and child behavior problems, while the concurrent association was found to be stronger than the longitudinal one (see de Roo et al., Reference de Roo, Veenstra and Kretschmer2022; Pinquart, Reference Pinquart2017a, Reference Pinquart2017b, for reviews). However, not all children experiencing poor parenting develop behavior problems. Theories such as the adaptive calibration model (ACM) of stress responsivity and biological sensitivity to the context (BSC) propose that the stress reactivity relevant to the HPA-axis may be a key moderator in these individual differences. The ACM emphasizes that beyond coordinating the organism’s allostatic response to physical and psychosocial challenges, another important function of the HPA axis is to encode and filter information coming from social and physical environments, thereby regulating the susceptibility of organisms and contributing to changes in physiological and behavioral outcomes (Del Giudice et al., Reference Del Giudice, Ellis and Shirtcliff2011). Similarly, the BSC theory (Boyce & Ellis, Reference Boyce and Ellis2005; Ellis & Boyce, Reference Ellis and Boyce2011), a theoretical foundation of ACM, suggests that the factors contributing to higher HPA-axis reactivity (e.g., hormones, genes, and epigenetic modifications) predispose individuals to being more susceptible to both negative and positive environmental influences. A great deal of evidence from basal cortisol, cortisol reactivity, or genetic variation relevant to the HPA-axis function has supported this “Stress reactivity × Environment” interaction hypothesis (e.g., Cao & Rijlaarsdam, Reference Cao and Rijlaarsdam2023; Kuhlman et al., Reference Kuhlman, Geiss, Vargas and Lopez-Duran2018). However, as a pivotal regulator of the HPA-axis function, NR3C1 methylation has been investigated as a biological mediator of stress-induced behavior problems (Chen & Cao, Reference Chen and Cao2024; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016), but has rarely been tested as a moderator for this interaction hypothesis.
To the best of our knowledge, there has been only one study investigating the “Stress reactivity × Environment” interaction involving NR3C1 methylation on child mental health. In a small sample of 46 German adolescents (aged 11–21 years old), Radtke et al. (Reference Radtke, Schauer, Gunter, Ruf-Leuschner, Sill, Meyer and Elbert2015) provided preliminary evidence for the moderation effect of NR3C1-1F methylation on the association between retrospectively reported childhood maltreatment and adolescent concurrent behavior problems comprising internalizing and externalizing symptoms. They revealed that, consistent with the hypothesis of ACM and BSC, adolescents with both higher methylation levels of NR3C1-1F (especially at canonical NGFI-A binding sites relevant to higher HPA-axis reactivity) and who were exposed to childhood maltreatment showed the most severe behavior problems. While not limited to mental health, another prospective longitudinal study conducted among 487 American children aged 7–16 years with a one-and-a-half-year interval reported a very similar, significant three-way interaction on child attachment (Bosmans et al., Reference Bosmans, Young and Hankin2018). That is, in a high-stress context (e.g., peer victimization) experienced during the interval, only children with higher NR3C1-1F methylation levels covering canonical NGFI-A binding sites and low maternal support reported higher levels of later anxious attachment. It should be noted that the prospective design adopted in this study outweighs a retrospective one used in the study of Radtke et al. (Reference Radtke, Schauer, Gunter, Ruf-Leuschner, Sill, Meyer and Elbert2015), as it greatly reduces recall biases and measurement errors for environmental exposures caused by retrospective assessment. Beyond this, it should also be kept in mind that the longitudinal research could provide an avenue for testing a potential, but rarely investigated, time-specific effect of the Epi × E interaction. This is crucial given the possible dynamic nature of epigenetic effects. For instance, a recent novel finding from an EWAS-based meta-analysis (Neumann et al., Reference Neumann, Sammallahti, Cosin-Tomas, Reese, Suderman, Alemany, Almqvist, Andrusaityte, Arshad, Bakermans-Kranenburg, Beilin, Breton, Bustamante, Czamara, Dabelea, Eng, Eskenazi, Fuemmeler, Gilliland and Cecil2025) revealed that, across child developmental outcomes (general psychopathology, ADHD, sleep, BMI, etc.), effect sizes of the concurrent effect of DNAm were larger compared to the longitudinal effect (“childhood DNAm → childhood outcomes” vs. “birth DNAm → childhood outcomes”). Particularly, a similar effect was found among CpG sites linked to corticosteroid hormone secretion involving ADHD (Neumann et al., Reference Neumann, Sammallahti, Cosin-Tomas, Reese, Suderman, Alemany, Almqvist, Andrusaityte, Arshad, Bakermans-Kranenburg, Beilin, Breton, Bustamante, Czamara, Dabelea, Eng, Eskenazi, Fuemmeler, Gilliland and Cecil2025). Unfortunately, no studies to date have used a (prospective) longitudinal design to test the extent to which the Epi × E effect regarding NR3C1 methylation varied across development (e.g., concurrent vs. longitudinal), especially during the preschool period, which is the first severity and prevalence peak (sensitive period) of behavior problems (Kavčič et al., Reference Kavčič, Podlesek and Zupančič2012).
On top of the potential dynamic nature, another important but unresolved issue that remains in extant research is how the Epi × E interaction works on child behavior problems, especially in terms of differential susceptibility or diathesis-stress. The traditional diathesis-stress model proposes a “risk” factor or “diathesis” conception (e.g., higher HPA-axis reactivity or related characteristics), which predisposes individuals to developmental dysfunction when exposed to adversity (Monroe & Simons, Reference Monroe and Simons1991). In contrast, the evolutionary-based differential susceptibility model reconstructs the “risk” nature of these factors into “plasticity” and argues that individuals with these factors could not only suffer from adverse environments but also benefit from supportive environments (Bakermans-Kranenburg & van IJzendoorn, Reference Bakermans-Kranenburg and van IJzendoorn2015). In fact, the aforementioned ACM and BSC theory support the differential susceptibility notion more strongly (Boyce & Ellis, Reference Boyce and Ellis2005; Del Giudice et al., Reference Del Giudice, Ellis and Shirtcliff2011; Ellis & Boyce, Reference Ellis and Boyce2011). Recently, the updated differential susceptibility 2.0 model further emphasizes a novel and very interesting perspective that the “plasticity” factors may vary across environmental exposure types (Belsky et al., Reference Belsky, Zhang and Sayler2022). For instance, the plasticity factor, such as a higher polygenic score, was only found to be correlated with greater susceptibility to the effects of child care but not the family environment (Zhang et al., Reference Zhang, Widaman and Belsky2023). However, whether the NR3C1 methylation interacts with caregiver parenting in a manner consistent with differential susceptibility or diathesis-stress has not been elucidated, although this would guide how we apply this biological factor to the prevention of and intervention for child behavior problems.
Besides, when considering the NR3C1 methylation × caregiver parenting interaction on child behavior problems, several moderators, especially parental sex, should not be neglected, especially in the context of Chinese culture. Although attachment theory emphasizes the equivalent importance of mothers and fathers in child development, empirical studies support that parental roles should account for ethnic or cultural customs and values (Cao & Rijlaarsdam, Reference Cao and Rijlaarsdam2023; Chuang & Su, Reference Chuang and Su2009). For instance, under the influence of Confucianism, Chinese fathers are often assigned a more central position in the family and primarily responsible for disciplining children, while Chinese mothers are mainly responsible for daily care and emotional companionship, especially in early life (Chuang & Su, Reference Chuang and Su2009). This cultural dynamic may drive a more pronounced Epi × E interaction effect regarding paternal parenting on child externalizing behaviors or maternal parenting on child internalizing behaviors. In addition, the moderating roles of child sex and NR3C1 genotypes could be explored. As early as preschool, boys are often at greater risk of externalizing behaviors and have higher glucocorticoid receptor expression, while girls show a higher – albeit not very robust – likelihood of developing internalizing behaviors and have higher plasma cortisol levels (Kokras et al., Reference Kokras, Hodes, Bangasser and Dalla2019; Wade et al., Reference Wade, Cairney and Pevalin2002). However, no significant sex differences in NR3C1 methylation alterations were preliminarily reported (Chen & Cao, Reference Chen and Cao2024). Besides, it remains unclear to what extent the NR3C1 BclI (rs41423247) and Tth111I (rs10052957) polymorphisms – two common single nucleotide polymorphisms (SNPs) that play an important role in HPA-axis reactivity and glucocorticoid sensitivity (Derijk, Reference Derijk2009; Spijker & Van Rossum, Reference Spijker and Van Rossum2012) – individually or additively moderate the Epi × E interaction. This rarely investigated allele-specific effect integrates the effects of epigenetic, genetic, and environmental factors together.
Current study
Although NR3C1 methylation has been suggested to help explain the individual differences in environmental susceptibility, the direct evidence, especially regarding the potential dynamic effects, function models, and moderating roles of parental and child sex and NR3C1 genotypes in the context of Chinese culture, is still lacking. Based on a one-and-a-half-year longitudinal design and multimodal data (epigenetic, genetic, and multi-informant behavioral data), the present study aimed to systematically investigate the interactive effect of NR3C1 methylation × caregiver parenting on child internalizing and externalizing behaviors in Chinese preschoolers. As shown in Figure 1, the following research questions were investigated: (i) whether there was a dynamic moderation effect involving NR3C1 methylation that varied between concurrent and longitudinal effects; (ii) how the Epi × E interaction functioned, especially in a manner consistent with differential susceptibility or diathesis-stress; and (iii) whether parental sex, child sex, and NR3C1 BclI and Tth111I polymorphisms (individually or additively) moderated the Epi × E interaction, particularly in Chinese families; it should be noted that of these three potential moderators, parental sex was the primary focus. Based on the novel findings regarding the dynamic effect of DNAm, the NR3C1 methylation might drive a more pronounced Epi × E interaction in the concurrent, but not longitudinal, effect. In light of relevant theories (e.g., ACM and BSC), the Epi × E interaction may function in a manner consistent with differential susceptibility such that children with higher versus lower NR3C1 methylation levels might not only suffer from low-quality parenting but also benefit from high-quality parenting, exhibiting a “for better and for worse” interaction nature. This effect may vary between parental sex in accordance with culture norms of Chinese families. However, we did not formulate hypotheses about the moderation effects of child sex or NR3C1 BclI and Tth111I polymorphisms due to the limited current evidence.
The conceptual model of hypotheses in the current study. Note. In the current study, a series of models were established to investigate the main and interaction effects of NR3C1 methylation and caregiver parenting on child internalizing and externalizing behaviors at T1 and T2, respectively. Through these models, three research questions were the primary focus: (i) whether there was a dynamic moderation involving NR3C1 methylation that varied between concurrent (Epi × E_T1 → outcomes_T1) and longitudinal effects (Epi × E_T1 → outcomes_T2); (ii) to which extent the Epi × E effect functioned in a manner consistent with differential susceptibility or diathesis-stress; and (iii) whether parental sex, child sex, and NR3C1 genotypes moderated the Epi × E effect in Chinese families; of these three moderators, parental sex was the primary one; the effect of NR3C1 genotypes included the individual and additive effects of BclI and Tth111I polymorphisms. Except for the model where the moderation effect of child sex was examined, child sex and age were entered as covariates; corresponding behavior problems at T1 were controlled for in the models at T2.
Method
Participants and procedure
Participants were drawn from a project investigating the epigenetic mechanisms underlying preschoolers’ mental health in Shandong Province, China. A few stress-related candidate genes were analyzed from libraries in this project, including NR3C1, NR3C2 (nuclear receptor subfamily 3 group C member 2), FKBP5 (FK506-binding protein 5), OXTR (oxytocin receptor), and BDNF (brain-derived neurotrophic factor). At baseline, a total of 224 families with children aged 36–72 months (M = 47.33 ± 9.60 months; 42.5% girls) were recruited from three kindergartens in Jinan city in September 2021. All children were Han Chinese and had no current diagnosis of any major physical, neurological, or pervasive developmental disorder. Table 1 shows the detailed characteristics of the sample at baseline.
Characteristics of the current sample (N = 224)
Note.
a 1,000¥ = 138USD.
All percentages are valid percentages.
Families were recruited through electronic brochures that were distributed in online class management groups. The entire data collection was conducted via one-on-one home visits. At baseline, both mothers and fathers were asked to complete questionnaires on rearing practices and child behavior problems, and children’s buccal cells were obtained by research assistants using a cheek swab kit (Rellagene, Shanghai, China) for DNA extraction. Children refrained from eating or drinking for 30 minutes prior to buccal cell collection. One and a half years later (T2), parents whose children were still in kindergarten completed the questionnaire assessing child behavior problems again (N = 113, M age = 63.60 ± 7.68 months; 45.7% girls). Small gifts (e.g., toys) were given to children for their participation. This study was approved by the University’s ethics committee. Informed consent was obtained from all parents and school principals prior to data collection.
Measures
Caregiver parenting
At T1, mothers and fathers evaluated their own and partner’s parenting quality using the widely used Chinese version of the Egna Minnen Beträffande Uppfostran for parents (EMBU-P; Castro et al., Reference Castro, de Pablo, Gómez, Arrindell and Toro1997). Three commonly investigated dimensions were assessed: warmth (7 items, e.g., “I/my partner have/has shown with words and gestures that I/my partner liked the child”), rejection (7 items, e.g., “I/my partner have punished the child even for small offenses”), and overprotection (9 items, e.g., “When my child has come back home, I/my partner always let him/her account for what he/she had been doing”). Each item was rated on a Likert scale ranging from 1 (not true at all) to 4 (almost always true). In this study, Cronbach’s α for the three dimensions ranged from .77 to .82. After rejection and overprotection items were reverse-coded, average scores for each dimension were first calculated with a maximum of 25% missing data allowed. Then, self-report and partner-report parenting scores for both mothers and fathers were computed by averaging standardized scores of three dimensions (r ≥ .14, p < .05; see Supplementary Table S1). Confirmatory factor analysis provided validity for this dimension combination (χ2/df ≤ 1.31, p < .05, RMSEA ≤ .04, SRMR ≤ .06, CFI ≥ .95, TLI ≥ .91). Finally, maternal parenting (M = −0.01 ± 0.87, range from −3.17 to 2.36) and paternal parenting (M = 0.01 ± 0.86, range from −2.39 to 2.91) were both derived by averaging their standardized self-report and partner-report scores (r ≥ .31, p < .001; see Supplementary Figure S3). High scores indicated high-quality maternal or paternal parenting. Notably, the mean scores and variances for maternal and paternal parenting scores of the current sample were higher than those of clinical (e.g., major depressive disorder) and some of the community child samples at a small-to-large effect size level (e.g., Cheng et al., Reference Cheng, Liu, Zheng, Jin, Wang, Liu and Wu2023; Kang et al., Reference Kang, Li, Zhao, Cui, Qin, Fang, Chen and Liu2024; p < .05; 0.32 ≤ Cohen’s d ≤ 1.56).
DNA methylation and genotype
Genomic DNA isolated from buccal swab samples was used for the measurement of DNAm and genotyping. The samples were assayed for quality by determining the concentration and purity using a Nanodrop 2,000 (NanoDrop Technologies, Wilmington, DE, USA). According to the manufacturer’s instructions, each DNA sample was subjected to bisulfite conversion using the EZ DNA Methylation-Gold Kit (ZYMO, CA, USA), which converted unmethylated cytosine bases to uracil. After polymerase chain reaction amplification (HotStart Taq polymerase, TaKaRa, Dalian, China) and library construction, samples were high-throughput sequenced (paired-end sequencing, 2 × 150 bp) using Illumina HiSeq/NovaSeq (Illumina, CA, USA). The sequencing primers and genomic locations of the assays were determined in accordance with Parade et al. (Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016). For each CpG site of the 1F region, the β value was used to represent its methylation levels, defined as the proportion of methylated signals among all signals (see Supplementary Table S1). The eight canonical NGFI-A binding sites in the region (see Supplementary Figure S1 and S2) were the primary focus of the current study (hereinafter referred to as NR3C1 methylation; McGowan et al., Reference McGowan, Sasaki, D’Alessio, Dymov, Labonté, Szyf, Turecki and Meaney2009), and their average methylation was used for primary analyses to avoid the correction for multiple testing (Chen & Cao, Reference Chen and Cao2024; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016)Footnote 1 .
The BclI (rs41423247) and Tth111I (rs10052957) polymorphisms were genotyped using a Sequenom (San Diego, CA, USA) chip-based MALDI-TOF mass spectrometry platform according to standard techniques (Cao & Rijlaarsdam, Reference Cao and Rijlaarsdam2023). Genotype frequencies of BclI (GG, CG, or CC) were 63.1%, 32.4%, or 4.5%, and those of Tth111I (GG, AG, or AA) were 83.3%, 16.2%, or 0.5%, both of which were in accordance with Hardy–Weinberg equilibrium (p > .05). The minor allele frequencies were both greater than 5% and comparable to frequencies identified in Asian ancestry samples. The haplotype analysis showed no linkage disequilibrium between BclI and Tth111I (D’ = .18, r 2 = .21; Carlson et al., Reference Carlson, Eberle, Rieder, Yi, Kruglyak and Nickerson2004). Given that the C allele of BclI and the A allele of Tth111I might be susceptible alleles associated with higher HPA-axis reactivity (Derijk, Reference Derijk2009; Spijker & Van Rossum, Reference Spijker and Van Rossum2012), these two polymorphisms were both coded in a dominant way (BclI: 0 = GG homozygote vs. 1 = at least one C allele; Tth111I: 0 = GG homozygote vs. 1 = at least one A allele), as in other studies (e.g., Chen & Cao, Reference Chen and Cao2024). Moreover, according to previous studies (e.g., Chen & Cao, Reference Chen and Cao2024; Davies et al., Reference Davies, Cicchetti and Hentges2015), a genetic composite score (GCS) was calculated by combining these two dominant codes and used to investigate the additive moderation effects of these two SNPs.
Internalizing and externalizing behaviors
The Chinese version of the Child Behavior Checklist for ages 1.5–5 (CBCL/1.5–5; Achenbach & Rescorla, Reference Achenbach and Rescorla2000) was used for both mothers and fathers to measure their children’s internalizing and externalizing behaviors at T1 and T2, respectively. The checklist has a total of 60 items, including two subscales: internalizing (36 items; emotionally reactive, anxious-depressed, somatic complaints, and withdrawn subscales; e.g., “Unhappy, sad, or depressed”) and externalizing (24 items; attention problems and aggressive subscales; e.g., “Can’t sit still”). Each item was rated from 0 (not true) to 2 (very true or often true). The Cronbach’s α for the two subscales ranged from .79 to .93, indicating acceptable to excellent reliability. Given significant intercorrelations (T1: r ≥ .20, p < .01; T2: r ≥ .41, p < .001), average scores of mother and father reports were calculated to obtain the final scores of internalizing behaviors (T1: 7.86 ± 5.23; T2: 3.24 ± 3.81; range from 0 to 72) and externalizing behaviors (T1: 10.22 ± 5.92; T2: 4.43 ± 4.86; range from 0 to 48) at two waves. Higher scores reflected higher levels of child behavior problems.
Missing data
In the current sample, 5 (2.2%) children had missing data on paternal parenting at T1. Mainly due to transition to primary school, 111 (49.6%) children did not have data on behavior problems at T2. Little’s missing completely at random test indicated that the missingness of data was not missing completely at random, χ2(47) = 75.29, p = .005. Nonrespondents at T2 were significantly older than respondents (t = 6.76, p < .001), but did not differ in any other study variables. Child age (the predictor of missingness) was thus controlled for in all subsequent models to fit missing at random assumptions. As such, missing data were handled using the gold standard approach of the full information maximum likelihood (Chen & Cao, Reference Chen and Cao2024), which allowed all 224 participants to be retained in the primary analyses (Little et al., Reference Little, Jorgensen, Lang and Moore2014).
Statistical analyses
Primary analyses
Three main steps were conducted in our primary analyses using Mplus 8.3. First, descriptive statistics and correlation analyses were conducted for study variables and covariates. Second, a series of nested models were established to investigate the main and interaction effects of NR3C1 methylation and caregiver parenting on child internalizing and externalizing behaviors at T1 and T2. Notably, while testing these Epi × E interactions, we incorporated the moderation effect of parental sex simultaneously. Specifically, (i) the freely estimated models, where maternal parenting, paternal parenting, and their respective interaction with NR3C1 methylation were all included, were compared with (ii) the constrained models, where regression weights of parenting and methylation × parenting were constrained to be equal between mothers and fathers, via the Wald chi-square test. If the omnibus Wald chi-square test was significant, a post hoc test was conducted to reveal which effect (main or interaction effect) significantly differed between mothers and fathers; otherwise, a nonsignificant statistic indicated no moderation effect of parental sex. Thus, via these model comparisons, (iii) the final models were established. Child sex and age were entered as covariates (Cao & Rijlaarsdam, Reference Cao and Rijlaarsdam2023); corresponding behavior problems at T1 were controlled for in the models at T2. Continuous independent variables were mean-centered prior to analysis. To reduce the effect of data bias on model estimation caused by nonnormal distributions, especially in methylation levels, robust maximum likelihood estimation with robust standard errors (MLR) was used. Finally, simple slope analyses were performed to probe the significant interactions. More importantly, re-parameterized regressions were conducted to test to what extent the interactions would be consistent with hypotheses of differential susceptibility or diathesis-stress (see details on the Supplementary “Re-parameterized regression model” section; Widaman et al., Reference Widaman, Helm, Castro-Schilo, Pluess, Stallings and Belsky2012). Notably, if equal effects for mothers and fathers were found in the above nested models, average caregiver parenting was used in the subsequent simple slope analysis and re-parameterized regression model (Stocker et al., Reference Stocker, Masarik, Widaman, Reeb, Boardman, Smolen, Neppl and Conger2017).
Secondary analyses
The potential moderating roles of child sex and NR3C1 genotypes in the Epi × E interaction were explored in two steps of secondary analyses. First, the Wald chi-square test was conducted to explore the respective moderation effects of child sex, BclI and Tth111I polymorphisms by constraining regression weights of the interaction model to be equal between dummy-coded group moderators (i.e., sex: 0 = boys, 1 = girls; BclI: 0 = GG homozygote, 1 = at least one C allele; Tth111I: 0 = GG homozygote, 1 = at least one A allele). When testing the moderation effect of child sex, it was not controlled as a covariate. Second, the additive moderation effects of these two SNPs were tested by constructing the three-way interaction term of GCS × methylation × parenting (ΔR 2).
In both primary nested models and secondary moderation analyses, Bonferroni correction was conducted to control for Type-I error in testing the interaction effects of methylation × parenting on child behavioral problems; the corrected threshold for the p-value was .0125 (i.e., .05/4; internalizing and externalizing behaviors at T1 and T2, respectively).
Sensitivity analyses
First, the primary Epi × E findings were tested again (i) among children with data available for both waves (N = 113); (ii) when child behavior problems reported by fathers and mothers were analyzed as the dependent variables, respectively; (iii) after further controlling for parental depressive symptoms at T1; and (iv) when internalizing and externalizing behaviors at the same wave were controlled for each other. Alternatively, the model was also run on child behavior problems at T2 without controlling for corresponding behavior problems at T1. Besides the focused canonical NGFI-A binding sites in our primary analyses, the roles of average methylation of the entire exon 1F region (i.e., 46 CpG sites) and the five noncanonical NGFI-A binding sites were explored. The roles of individual parenting dimensions were tested. Additionally, given preliminary evidence for the biological embedding mechanism of NR3C1 methylation underlying stress-induced psychopathology (Chen & Cao, Reference Chen and Cao2024; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016), we also tested the mediation effect of NR3C1 methylation in the parenting–behavior problems link. Finally, a machine learning-based k−fold cross validation analysis (Kuhn et al., Reference Kuhn, Jed, Steve, Andre, Chris, Allan and Tyler2023; Melhem et al., Reference Melhem, Porta, Oquendo, Zelazny, Keilp, Iyengar, Burke, Birmaher, Stanley, Mann and Brent2019) was further used to test the robustness of the current findings (see details in the Supplementary “Method” section).
Results
Preliminary analyses
Descriptive statistics and bivariate correlations among all study variables are shown in Table 2. Internalizing and externalizing behaviors were correlated with each other within each wave (r ≥ .75, p < .001), and both declined dramatically from T1 to T2 (8.79 ≤ t ≤ 9.46, p < .001; large effect: Cohen’s d = 0.83 and 0.89), consistent with previous findings regarding the developmental trajectory of behavior problems during the preschool period (Kavčič et al., Reference Kavčič, Podlesek and Zupančič2012). NR3C1 methylation was correlated with maternal parenting (r = .13, p < .05), but not with paternal parenting or any child behavior problems (−.09 ≤ r ≤ .11, p > .05). The dummy-coded BclI and Tth111I polymorphisms were not individually or additively related to any study variables (−.11 ≤ r ≤ .10, p > .05). Maternal and paternal parenting (intercorrelation: r = .63, p < .001) were negatively correlated with internalizing and externalizing behaviors at T1 (r ≥ .30, p < .001), but not at T2 (−.09 ≤ r ≤ .12, p > .05). Boys were more likely than girls to exhibit more externalizing behaviors at T1 (r = −.15, p < .05).
Descriptive statistics and correlations among study variables
Note. N T1 = 219; N T2 = 113; mainly due to the primary school entrance, 111 (49.6%) children did not have data on behavior problems at T2.
Sex: 0 = boys, 1 = girls; Age = age of children at T1 (months); MPQ = maternal parenting quality; PPQ = paternal parenting quality; APQ = average parenting quality; methylation = average methylation of canonical NGFI-A binding sites; BclI: 0 = GG homozygote, 1 = at least one C allele; Tth111I: 0 = GG homozygote, 1 = at least one A allele; GCS = the genetic composite score of BclI and Tth111I polymorphisms, ranging from 0 to 2; INT = internalizing behaviors; EXT = externalizing behaviors.
*p < .05, **p < .01, ***p < .001.
NR3C1 methylation × caregiver parenting interaction on internalizing behaviors
After model comparisons (1.36 ≤ Wald statistic ≤ 8.32, p > .0125), both constrained models were chosen as final models for internalizing behaviors at T1 and T2, where the main effect of parenting and the interaction effect of methylation × parenting did not significantly differ between parental sex.
Models regarding internalizing behaviors at T1
As shown in Table 3, after controlling for child sex and age, caregiver parenting (β = −.20, p < .001), rather than NR3C1 methylation (β = −.04, p = .596), was negatively associated with child internalizing behaviors at T1. The interaction of NR3C1 methylation × caregiver parenting was significant (β = .09, p = .002, ΔR 2 = .022). Simple slope analysis showed that (see Figure 2a), only among children with lower but not higher NR3C1 methylation levels, caregiver parenting negatively predicted child internalizing behaviors at T1 (methylation at M − 2SD (value = −0.72): B (SE) = −4.34 (0.82), t = −5.32, p < .001 vs. methylation at M + 2SD (value = 0.72): B (SE) = −0.51 (0.71), t = −0.72, p = .472). Exploratory individual-site analyses revealed that this interaction was mainly manifested at sites of CpG 20 and 21 (see Supplementary Table S2).
Interactions between NR3C1 methylation and caregiver parenting on child behavior problems at T1. Note. Canonical NGFI-A methylation × caregiver parenting quality on (a) child internalizing behaviors and (b) externalizing behaviors at T1. Simple slopes for canonical NGFI-A methylation were plotted at M ± 2SD (−0.72 vs. 0.72). Gray shaded area represents the 95% CI of the cross-over point C of the interactions on the maternal/paternal parenting quality axis. The 95% CI of C on internalizing and externalizing behaviors ranged from –0.16 to 1.22 and from –0.17 to 0.71, respectively. The caregiver parenting quality axis is the average score of maternal and paternal parenting quality, ranging from −2.43 to 2.34.
Regression models predicting child behavior problems from interactions between NR3C1 methylation and caregiver parenting
Note. 224 participants were retained in all models given that missing data were handled using the gold standard approach of full information maximum likelihood.
Sex: 0 = boys, 1 = girls; Age = age of children at T1 (months); Methylation = average methylation of canonical NGFI-A binding sites.
MPQ = maternal parenting quality; PPQ = paternal parenting quality; INT = internalizing behaviors; EXT = externalizing behaviors.
a MPQ/PPQ means the constrained, equal effect regarding mothers and fathers.
b The corrected threshold of p-value for the interaction effects was .0125.
Re-parameterized regression models (see Supplementary Table S3) demonstrated that this Epi × E interaction aligned more with differential susceptibility than diathesis-stress. The point estimate and 95% CI of the crossover point (C = 0.57, 95% CI [−0.09, 1.22]) both fell within the range of average caregiver parenting quality (range from −2.43 to 2.34). The model comparison supported that differential susceptibility outweighed diathesis-stress (ΔR 2 = .018, p = .035). Thus, lower-level NR3C1 methylation might be a “plasticity” rather than “risk” factor: children with lower versus higher NR3C1 methylation levels not only suffered from low-quality parenting, exhibiting more internalizing behaviors, but also benefited from high-quality parenting, exhibiting fewer internalizing behaviors (see Figure 2a).
Models regarding internalizing behaviors at T2
As shown in Table 3, only child internalizing behaviors at T1 significantly predicted internalizing behaviors at T2 (β = .34, p = .008), whereas no other effects regarding NR3C1 methylation, caregiver parenting, or their interaction were significant (−.05 ≤ β ≤ .06, p > .0125). The findings at both waves suggest a time-specific (concurrent but not longitudinal) Epi × E interaction on child internalizing behaviors during the preschool period.
NR3C1 methylation × caregiver parenting interaction on externalizing behaviors
Similarly, both constrained models were chosen as final models for externalizing behaviors at T1 and at T2, where the main effect of parenting and the interaction effect of methylation × parenting did not significantly differ between parental sex (2.01 ≤ Wald statistic ≤ 6.60, p > .0125).
Models regarding externalizing behaviors at T1
As shown in Table 3, after controlling for child sex and age, caregiver parenting (β = −.24, p < .001) rather than NR3C1 methylation (β = .02, p = .800) was negatively associated with child externalizing behaviors at T1. The interaction of NR3C1 methylation × caregiver parenting was significant (β = .12, p < .001, ΔR 2 = .046). As shown in Figure 2b, among those with lower but not higher NR3C1 methylation levels, caregiver parenting was negatively associated with child externalizing behaviors at T1 (methylation at M − 2SD (value = −0.72): B (SE) = −6.51 (0.82), t = −7.96, p < .001 vs. methylation at M + 2SD (value = 0.72): B (SE) = −0.20 (0.63), t = −0.32, p = .749). This interaction was mainly manifested at sites of CpG 17, 20, and 21 (see Supplementary Table S2).
As suggested by the re-parameterized regression models (see Supplementary Table S3), this interaction was also in line with the hypothesis of differential susceptibility but not diathesis-stress (Figure 2b). Specifically, the point estimate and 95% CI of the crossover point of the interaction (C = 0.27, 95% CI [−0.17, 0.71]) were within the range of observed values of average caregiver parenting quality. Meanwhile, the differential susceptibility model significantly fitted the data better than the diathesis-stress model (ΔR 2 = .043, p < .001). The finding supports the plastic effect of lower-level NR3C1 methylation toward caregiver parenting with regard to child externalizing behaviors.
Models regarding externalizing behaviors at T2
As shown in Table 3, child externalizing behaviors at T2 were significantly predicted by corresponding behaviors at T1 (β = .35, p = .005). No other effects regarding NR3C1 methylation, caregiver parenting, or their interaction fell into significance (−.06 ≤ β ≤ .08, p > .0125). Again, the findings at both waves suggest a time-specific (concurrent but not longitudinal) Epi × E interaction on child externalizing behaviors during the preschool period.
Secondary analyses: moderation effects of child sex and NR3C1 genotypes
For either behavior problems at T1 or T2, the Epi × E interaction did not vary between child sex (0.82 ≤ Wald statistic ≤ 5.91, p > .0125; see Table 4); meanwhile, NR3C1 BclI and Tth111I polymorphisms did not individually or additively related to the NR3C1 methylation (−.11 ≤ r ≤ .02, p > .05) or moderate the Epi × E interaction (0.31 ≤ Wald statistic ≤ 4.45, .001 ≤ ΔR 2 ≤ .003, p > .0125).
Regression models testing the moderation effects of child NR3C1 BclI and Tth111I polymorphisms
Note. 224 participants were retained in all models given that missing data were handled using the gold standard approach of full information maximum likelihood.
Moderators: child sex: 0 = boys, 1 = girls; BclI: 0 = GG homozygote, 1 = at least one C allele; Tth111I: 0 = GG homozygote, 1 = at least one A allele.
GCS = the genetic composite score of BclI and Tth111I polymorphisms, ranging from 0 to 2; INT = internalizing behaviors; EXT = externalizing behaviors.
a When testing the moderation effect of child sex, only child age was included as a covariate.
b When testing the moderation effects of NR3C1 genotypes, child sex and age were included as covariates.
c The corrected threshold of p-value for the interaction effects was .0125.
Sensitivity analyses
First, our primary findings of the concurrent, but not longitudinal, Epi × E effect remained robust (i) among children with data available for both waves (N = 113; .18 ≤ β ≤ .21, p < .001, .060 ≤ ΔR 2 ≤ .080; Table S4) and (ii) after controlling for parental depressive symptoms at T1 (.09 ≤ β ≤ .12, p < .01, .022 ≤ ΔR 2 ≤ .040; Table S5). While using indicators of mother and father reports, respectively, the Epi × E findings did not vary between mothers vs. fathers or boys vs. girls (0.01 ≤ Wald statistic ≤ 9.02, p > .0125); the Epi × E effect was significant on child externalizing (.08 ≤ β ≤ .09, p < .05, .006 ≤ ΔR 2 ≤ .016), but not internalizing behaviors (.04 ≤ β ≤ .06, p > .05, .000 < ΔR 2 ≤ .009; Table S6), at a p < .05 level. Again, such a more pronounced effect on externalizing behaviors was also reported when child behavior problems were controlled for each other (β = .07, p < .001, ΔR 2 = .013; Table S7). Meanwhile, the Epi × E interaction remained nonsignificant on child behavior problems at T2, even when not controlling for their corresponding behaviors at T1 (.04 ≤ β ≤ .05, p > .0125, .007 ≤ ΔR 2 ≤ .009; Table S8).
Third, average methylation levels of the entire exon 1F region, but not the noncanonical NGFI-A binding sites, still drove a significant but smaller interaction effect on child externalizing behaviors at T1 (β = .07, p = .004; .004 ≤ ΔR 2 ≤ .018; Tables S9 and S10). Fourth, the Epi × E findings were replicated in most of the individual parenting dimensions (Epi × E involving rejection and overprotection → internalizing; Epi × E involving warmth, rejection, and overprotection → externalizing; .017 ≤ ΔR 2 ≤ .040; Table S11). Fifth, we failed to capture the mediation effect of NR3C1 methylation between caregiver parenting and child behavior problems at both waves (all 95% CIs of indirect effects included 0; Table S12). Finally, the k−fold cross-validation test showed that child internalizing and externalizing behaviors at T1 (.35 ≤ r ≤ .45; p < .001) but not T2 (.19 ≤ r ≤ .21; p > .0125) predicted by the Epi × E interaction models in the training sets were positively correlated with the actual child internalizing and externalizing behaviors scores in the test sets. That is, the current findings of the time-specific (concurrent but not longitudinal) Epi × E interaction on child behavioral problems were validated internally and had an acceptable generalization ability.
Discussion
Although NR3C1 methylation has been suggested to help explain the individual differences in environmental susceptibility, the direct evidence, especially regarding the dynamic effects, function models, and moderation of parental and child sex and NR3C1 genotypes, is still lacking. Based on a one-and-a-half-year longitudinal design and multimodal data (epigenetic, genetic, and multi-informant behavioral data), the current study was the first to systematically investigate whether and how NR3C1 methylation interacted with caregiver parenting to affect the development of internalizing and externalizing behaviors among Han Chinese preschoolers. Three main findings are highlighted here. First, NR3C1 methylation significantly moderated the concurrent, rather than longitudinal, relationship between caregiver parenting and child internalizing and externalizing behaviors, suggesting a dynamic Epi × E interaction nature. Two sites, CpG 20 and CpG 21, contributed more to this effect. Second, this interaction functioned in a manner consistent with differential susceptibility. Children with lower versus higher NR3C1 methylation levels exhibited not only more behavior problems when exposed to low-quality parenting but also fewer behavior problems when exposed to high-quality parenting. Third, such a time-specific Epi × E effect did not vary between parental sex, child sex, and NR3C1 BclI and/or Tth111I polymorphisms.
To date, the vast majority of previous research has primarily adopted a cross-sectional design and focused more on the concurrent effect of epigenetics on child behavior problems (see Watkeys et al., Reference Watkeys, Kremerskothen, Quide, Fullerton and Green2018, for a review), whereas little is known about the extent to which DNAm functions in a dynamic way across development, especially during some sensitive periods. One recent study preliminarily showed that DNAm involving corticosteroid hormone secretion exerted a larger concurrent main effect than the longitudinal effect on ADHD (“childhood DNAm → childhood outcome” vs. “birth DNAm → childhood outcome”; Neumann et al., Reference Neumann, Sammallahti, Cosin-Tomas, Reese, Suderman, Alemany, Almqvist, Andrusaityte, Arshad, Bakermans-Kranenburg, Beilin, Breton, Bustamante, Czamara, Dabelea, Eng, Eskenazi, Fuemmeler, Gilliland and Cecil2025). Based on the multimodal data from a two-wave longitudinal design, the current study extended this interesting time-specific epigenetic main effect into the Epi × E domain and provided the first evidence that NR3C1 methylation only interacted with caregiver parenting on the concurrent, but not subsequent (e.g., one and a half years later) or prospective changes in child internalizing and externalizing behaviors in the preschool period.
Although the mechanism is far from clear, the explanation of this intriguing finding may need to take into account (i) the stages of development and (ii) the dynamics of DNAm itself. For instance, the preschool period is a crucial period for the development of behavior problems, such that child internalizing and externalizing behaviors show their first severity and prevalence peak at around the age of 3 and then decline subsequently until adolescence (Kavčič et al., Reference Kavčič, Podlesek and Zupančič2012). The first four years are also pivotal for the growth of a key stress regulation area – the hippocampus, which is rich in glucocorticoid receptors and greatly affected by NR3C1 methylation (Eichenbaum, Reference Eichenbaum2013). Thus, the dynamic Epi × E effect may be a mapping of “the closer you get to the sensitive period, the greater the effect will be.” Alternatively, regardless of the specific period, the shrinking longitudinal Epi × E effect may be induced by (potentially large) dynamic changes in DNAm across development. Recall that preliminary evidence reports that NR3C1 methylation shows weak stability, with only a small to medium correlation between intervals from 3 months to 2 years (range: −.008 to .33; e.g., Di Sante et al., Reference Di Sante, Ismaylova, Nemoda, Gouin, Yu, Caldwell, Vitaro, Szyf, Tremblay and Booij2018). This possibility warrants repeated measures of DNAm data, whereas, unfortunately, we were unable to validate it. Besides, the confounding effect of some unmeasured variables should not be ignored; for instance, some other environmental factors (e.g., peer victimization; Watkeys et al., Reference Watkeys, Kremerskothen, Quide, Fullerton and Green2018) may also drive the changes in NR3C1 methylation. Taken together, the current finding preliminarily highlights a very interesting time-specific Epi × E effect involving NR3C1 methylation; we call for more longitudinal studies with multiple repeated measures of DNAm, especially at ideally more sensitive developmental periods, to better capture its dynamic nature.
This study further contributed to the field by revealing, for the first time, that the lower-than-average levels of NR3C1 methylation that are tightly associated with lower HPA-axis reactivity may be a plasticity factor, exerting the Epi × E interaction in a manner consistent with differential susceptibility rather than diathesis-stress. This finding was only partially consistent with the ACM and BSC theories, as both highlight that higher, but not lower, HPA-axis reactivity would show greater plasticity in response to the environment (Del Giudice et al., Reference Del Giudice, Ellis and Shirtcliff2011; Ellis & Boyce, Reference Ellis and Boyce2011). Despite the fact that testing this interaction hypothesis from epigenetics is still in its infancy, both “stress hypo-reactivity × Environment” and “stress hyper-reactivity × Environment” findings have been reported in research regarding HPA-axis reactivityFootnote 2 (e.g., basal cortisol or cortisol reactivity; Kuhlman et al., Reference Kuhlman, Geiss, Vargas and Lopez-Duran2018; Kushner et al., Reference Kushner, Barrios, Smith and Dougherty2016). The lower HPA-axis reactivity has been proposed to be a long-term maladaptive outcome of chronic and repeated stress exposure, whereas the opposite, higher HPA-axis reactivity, is more likely induced by acute stress (Tarullo & Gunnar, Reference Tarullo and Gunnar2006). Furthermore, these stress hypo-reactivity alterations, including DNAm, may have pleiotropy effects: (i) poor coping and recovery difficulties in response to subsequent stress (i.e., “for worse” effect; Kushner et al., Reference Kushner, Barrios, Smith and Dougherty2016), but (ii) self-protection of the brain (the hippocampus region in particular) against the toxic effects of hypercortisolism, especially under supportive environments (i.e., “for better” effect; Kuhlman et al., Reference Kuhlman, Geiss, Vargas and Lopez-Duran2018). Alternatively, as proposed by the differential susceptibility 2.0 model (Belsky et al., Reference Belsky, Zhang and Sayler2022), lower-than-average levels of NR3C1 methylation are likely a “specific” plasticity factor toward mild, rather than severe (e.g., maltreatment; Radtke et al., Reference Radtke, Schauer, Gunter, Ruf-Leuschner, Sill, Meyer and Elbert2015), rearing environments. Notably, although few Epi × E studies have tested this crossover interaction, there is increasing evidence regarding characteristics of polygenic scores, difficulty temperament, or sensory sensitivity that supports this evolutionary-based view of differential susceptibility, rather than the traditional diathesis-stress model (Assary et al., Reference Assary, Krebs and Eley2023; Bakermans-Kranenburg & van IJzendoorn, Reference Bakermans-Kranenburg and van IJzendoorn2015). However, one caveat is that the “plasticity” factors may vary across environmental exposure or person characteristic types (Belsky et al., Reference Belsky, Zhang and Sayler2022), as well as across the observed environment range based on the samples (Xu et al., Reference Xu, Yang and Cao2025). More work is encouraged to deepen the understanding of the Epi × E pattern regarding NR3C1 methylation and parenting in this nascent field.
Contrary to our hypothesis, the current study reported that neither the main effect of parenting nor the interaction effect of methylation × parenting varied between parental sex, regardless of who reported the child behavior problems. Relevant to this finding, the role differences between Chinese mothers and fathers in the family have been found to become dramatically weaker in recent years, especially with the significant increase in women’s educational levels and the influence of multiculturalism (e.g., gender equality) in China (Li & Lamb, Reference Li, Lamb and Roopnarine2016). Moreover, the Epi × E effect did not vary between child sex and was not individually or additively moderated by NR3C1 BclI and/or Tth111I polymorphisms. Notably, despite not being found in our study, there is evidence supporting an allele-specific effect in epigenetic moderation regarding the FKBP5 gene (another important stress-related gene; Mulder et al., Reference Mulder, Rijlaarsdam, Luijk, Verhulst, Felix, Tiemeier, Bakermans-Kranenburg and Van Ijzendoorn2017). The roles of these potential moderators – parental and child sex, and NR3C1 genotypes – should be further investigated, especially in other groups beyond the current Han Chinese sample. Consistent with previous research (see Watkeys et al., Reference Watkeys, Kremerskothen, Quide, Fullerton and Green2018, for a review), no allele-specific effect of BclI and Tth111I polymorphisms on the NR3C1 methylation was reported. Future studies with larger sample sizes and more comprehensive genotyping – examining not only polymorphisms located directly at the CpG sites themselves, but also those in upstream or downstream regulatory regions (Gaunt et al., Reference Gaunt, Shihab, Hemani, Min, Woodward, Lyttleton, Zheng, Duggirala, McArdle, Ho, Ring, Evans, Davey Smith and Relton2016) – are warranted to validate and extend the current finding. Longitudinal methylation data with repeated sampling would provide greater leverage in disentangling environmentally induced epigenetic changes from those driven by stable genetic variation.
Besides, the current findings highlight the effect of canonical, rather than noncanonical, NGFI-A binding sites, whereas when considering the effect of the entire 1F region, the Epi × E effect sizes were reduced in magnitude. Methylation alterations of CpG 20 and CpG 21, in particular, contributed more to this dynamic Epi × E effect. As noted, the inhibitory effect of NR3C1-1F methylation on glucocorticoid receptor expression is primarily driven by canonical NGFI-A binding sites (McGowan et al., Reference McGowan, Sasaki, D’Alessio, Dymov, Labonté, Szyf, Turecki and Meaney2009). Previous research has emphasized the effect of methylation at these sites on psychopathology. For instance, among 39 CpG sites (from CpG 9 to CpG 47) of NR3C1-1F, Bakusic et al. (Reference Bakusic, Vrieze, Ghosh, Bekaert, Claes and Godderis2020) found that only methylation of CpG 20 significantly predicted a reduction in depressive symptoms among adults over an 8-week interval. Methylation alterations of both CpG 20 and CpG 21 have been reported to be more likely affected by stress exposures (e.g., maternal anxiety; Hompes et al., Reference Hompes, Izzi, Gellens, Morreels, Fieuws, Pexsters, Schops, Dom, Van Bree, Freson, Verhaeghe, Spitz, Demyttenaere, Glover, Van den Bergh, Allegaert and Claes2013). Such potential differential effects across individual CpG sites call for a more targeted investigation of their molecular mechanisms. Given that parental psychopathology is often significantly correlated with unfavorable rearing practices, child epigenetic modifications, and behavior problems (e.g., Aaron et al., Reference Aaron, Kaplan and Black2024; Hill et al., Reference Hill, Pickles, Wright, Quinn, Murgatroyd and Sharp2019), we checked the Epi × E findings after controlling for parental depressive symptoms. The results remained robust, highlighting the unique effect driven by caregiver parenting. Nevertheless, due to the complex nature of child behavior problems, it is necessary to consider more factors involving parental psychopathology or psychological stress when exploring the biological effect of epigenetics on the parenting–behavior problems link. In addition, a more robust interaction effect was revealed on externalizing, rather than internalizing, behaviors when these problems were controlled for each other or adopted indicators of mother and father reports, respectively. This might be due to a more stable relationship between HPA-axis function and externalizing behaviors during preschool (Koss & Gunnar, Reference Koss and Gunnar2018). Despite the mediation effect of NR3C1 methylation from parenting to child behavior problems was not captured, this possibility has been proposed by the biological embedding model (Miller et al., Reference Miller, Chen and Parker2011) and reported in some studies (e.g., Chen & Cao, Reference Chen and Cao2024; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016). Prospective longitudinal research across multiple developmental phases (e.g., preschool, adolescence) is beneficial to further ascertain whether the roles (mediator and/or moderator) of NR3C1 methylation differ across different periods or have conditional effects.
This study is the first to examine the interaction of NR3C1 methylation × caregiver parenting on child internalizing and externalizing behaviors in the context of Chinese culture, which has several strengths. First, the use of a one-and-a-half-year longitudinal design allowed us to test the time-specific moderation of NR3C1 methylation, facilitating the understanding of the dynamic Epi × E effect in the preschool period. Second, this study highlighted how the lower levels of NR3C1 methylation functioned in a manner consistent with differential susceptibility with standardized statistical tests for the first time. Third, by adopting multimodal data, the rarely investigated moderating roles of parental and child sex and NR3C1 genotypes were explored, particularly in the Han Chinese sample. Fourth, multi-informant measures for both caregiver parenting and child behavior problems reduced reporting bias to some extent. Finally, a series of rigorous analyses were conducted to expand the findings and ensure their robustness; for instance, the role of each individual CpG site and the alternative mediation effect of NR3C1 methylation were tested.
However, there are some limitations that need to be considered cautiously. First, the lack of repeated data, especially NR3C1 methylation, across multiple waves and/or across multiple developmental phases (e.g., preschool, adolescence) limits our understanding of the time-specific nature of the Epi × E interaction. For instance, it remains unclear to what extent the current findings can be generalized to other developmental stages or are specific to the preschool period. The sensitive window of temporal dynamics of DNAm itself remains poorly understood (Boyce & Kobor, Reference Boyce and Kobor2015). Moreover, although we tested the Epi × E effect on child behavioral problems, causal inference cannot be established without longitudinal repeated data, as the bidirectional association may exist between epigenetic modifications and behavioral problems (Klengel et al., Reference Klengel, Pape, Binder and Mehta2014). Besides, incorporating downstream indicators of DNAm, such as glucocorticoid receptor mRNA expression, blood glucocorticoid levels, inflammatory biomarkers, and neural stress reactivity, is encouraged (e.g., McGowan et al., Reference McGowan, Sasaki, D’Alessio, Dymov, Labonté, Szyf, Turecki and Meaney2009; Na et al., Reference Na, Chang, Won, Han, Choi, Tae, Yoon, Kim, Joe, Jung, Lee, Ham and DeAngelis2014; van der Knaap et al., Reference van der Knaap, Oldehinkel, Verhulst, van Oort and Riese2015). For instance, whether the Epi × E effect functions through altered blood glucocorticoid and IL-6 and TNF-α levels or via stress regulation in relevant brain areas (e.g., amygdala, prefrontal cortex, or hippocampus) needs to be tested to clarify the involved pathway (Na et al., Reference Na, Chang, Won, Han, Choi, Tae, Yoon, Kim, Joe, Jung, Lee, Ham and DeAngelis2014; van der Knaap et al., Reference van der Knaap, Oldehinkel, Verhulst, van Oort and Riese2015; Weaver et al., Reference Weaver, Cervoni, Champagne, D’Alessio, Sharma, Seckl, Dymov, Szyf and Meaney2004). Second, this study mainly focused on the NR3C1-1F region (NGFI-A binding sites in particular) – one of the most investigated regions in extant research regarding the methylation mechanism of behavior problems (Cecil et al., Reference Cecil, Neumann and Walton2023). Nevertheless, the effects of other regions (e.g., 1H and 1D; Palma-Gudiel et al., Reference Palma-Gudiel, Cordova-Palomera, Leza and Fananas2015) and other important candidate genes (e.g., FKBP5; Tyrka et al., Reference Tyrka, Ridout and Parade2016) should also be considered; more importantly, how the epigenetic modifications of these candidate genes work together (e.g., interactively or additively) to elucidate the complex, polygenic nature of the biological background of behavior problems needs to be explored (Huls & Czamara, Reference Huls and Czamara2020). Third, given its significant correlation with methylation levels in brain tissues (Braun et al., Reference Braun, Han, Hing, Nagahama, Gaul, Heinzman, Grossbach, Close, Dlouhy, Howard, Kawasaki, Potash and Shinozaki2019), noninvasive buccal swab samples were collected to assess DNAm, as performed in the vast majority of previous child research (Chen & Cao, Reference Chen and Cao2024). However, it is important to note that the proportion of epithelial cells in buccal samples shows great heterogeneity. For instance, although Theda et al. (Reference Theda, Hwang, Czajko, Loke, Leong and Craig2018) reported relatively high epithelial cell proportions (ranging from 75 to 95%) in 20 swab samples, Merrill et al. (Reference Merrill, Konwar, Fatima, Dever, MacIsaac, Letourneau, Giesbrecht, Dewey, England-Mason, Lewis, Wang, Teh, Meaney, Gonzalez, Noll, De Weerth, Bush, O’Donnell, Stewart and Kobor2025), based on 4,626 children’s samples, observed a broader range of 50 to 100%. Therefore, the current findings need to be considered with caution, and future studies would benefit from correcting this confounding effect; for instance, the bioinformatics deconvolution method can be adopted. Fourth, although both maternal and paternal reports of parenting and child behavior problems were utilized, the potential for reporting bias caused by shared method variance remains and may limit the interpretability of the findings. Future studies could incorporate more informants or observational assessments to provide more objective and ecologically valid measurements of these variables. Fifth, considering the overall high level of education of the present sample (e.g., 98.2% of mothers with college degrees), the extent to which the results can be generalized to other ethnic/cultural groups or high-risk children with low socioeconomic status needs to be verified. Finally, and importantly, a significant limitation is that the current study was not preregistered, which greatly limits transparency in analytic decision-making and replicability (Nosek et al., Reference Nosek, Ebersole, DeHaven and Mellor2018). We strongly recommend future research to adopt preregistered designs to minimize analytic biases and allow for more robust replication of the current findings.
Based on a one-and-a-half-year longitudinal design and multimodal data, the current study provided preliminary evidence that the NR3C1 methylation moderated the concurrent, but not longitudinal, associations between caregiver parenting and child internalizing and externalizing behaviors in Chinese families. Consistent with the differential susceptibility hypothesis, children with lower versus higher NR3C1 methylation not only suffered from low-quality parenting but also benefited from high-quality parenting. This time-specific Epi × E effect did not vary between parental and child sex and was not moderated by NR3C1 BclI and Tth111I polymorphisms individually or additively. The findings highlight the dynamic nature of Epi × E interaction in the preschool period and advance our knowledge of identifying NR3C1 methylation as a plasticity factor that works in a “for better and for worse” manner. These findings may help establish precise preventive interventions from a perspective of “for whom” and “when”. For instance, by targeting children with low levels of NR3C1 methylation, especially during the critical preschool years, preventive interventions (e.g., mindfulness-based parenting interventions) might achieve maximum efficiency with minimal effort to reduce child behavior problems.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0954579425100850.
Data availability statement
The data and code are available from the corresponding author upon reasonable request.
Acknowledgments
We are extremely grateful to all the families who participated in the current study, as well as the research assistants who contributed to data collection.
Funding statement
This study was supported by the National Natural Science Foundation of China (32571247, 31800936), Natural Science Foundation of Shandong Province (ZR2023MC118), and Future Plan of Young Scholars in Shandong University (21330089964217) to Cong Cao.
Competing interests
All authors declared that they have no conflicts of interest.
Ethics approval
Informed consent was obtained from all involved parents and school principals. The present study was approved by the Ethics Committee of the School of Nursing and Rehabilitation, Shandong University. All procedures performed in studies involving human participants were in accordance with the ethical standards of the ethics committee and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards.
Pre-registration statement
The present research was not pre-registered. It was conducted under the framework of a project funded by the National Natural Science Foundation of China and approved by the Ethics Committee of the School of Nursing and Rehabilitation, Shandong University.
Open access
