Each year in the United States four million reports of abuse involving over seven million children are made to child protective services (US Department of Health and Human Services, Administration on Children, Youth and Families, Children's Bureau, 2017). Exposure to early environmental adversity (e.g., abuse, trauma, or contextual stressors) confers substantial risk for the development of psychopathology and lifelong risks of chronic disorders of health and well-being (see Cicchetti & Toth, Reference Cicchetti, Toth and Cicchetti2016; Cohen, Janicki-Deverts, & Miller, Reference Cohen, Janicki-Deverts and Miller2007; Norman et al., Reference Norman, Byambaa, De, Butchart, Scott and Vos2012; Shonkoff, Boyce, & McEwen, Reference Shonkoff, Boyce and McEwen2009, for reviews). Despite the prevalence of early environmental adversity and clear scientific consensus linking this adversity to detrimental health outcomes, the mechanisms underlying this connection are not well understood (McCrory & Viding, Reference McCrory and Viding2015).
One process by which maltreatment and other early adversities may influence the development of later psychopathology is through modification of the physiologic stress response system, and in particular, the hypothalamus–pituitary–adrenal (HPA) axis (Heim & Binder, Reference Heim and Binder2012; Tyrka, Burgers, Philip, Price, & Carpenter, Reference Tyrka, Burgers, Philip, Price and Carpenter2013). In response to stressful stimuli, the HPA axis is activated, and glucocorticoids are released, exerting cellular responses by binding at the intracellular glucocorticoid receptor (GR). GRs are distributed throughout the body and brain, where they regulate basal physiologic function and promote adaptive responses to acute stressors (de Kloet, Joels, & Holsboer, Reference de Kloet, Joëls and Holsboer2005; Kadmiel & Cidlowski, Reference Kadmiel and Cidlowski2013). Further, activation of the GR through cortisol binding at the hypothalamus and pituitary engages a negative feedback loop that inhibits further release of cortisol and prevents damaging effects of extreme or chronic activation (Herman, McKlveen, Solomon, Carvalho-Netto, & Myers, Reference Herman, McKlveen, Solomon, Carvalho-Netto and Myers2012; Laryea, Muglia, Arnett, & Muglia, Reference Laryea, Muglia, Arnett and Muglia2015). However, studies of rodents, children, and adults demonstrate that severe (e.g., child maltreatment) or chronic (e.g., contextual stressors associated with poverty) early life stress produces long-term alterations in glucocorticoid signaling that contribute to “wear and tear” on this system across the life span (McEwen et al., Reference McEwen, Bowles, Gray, Hill, Hunter, Karatsoreos and Nasca2015; Tyrka, Ridout, & Parade, Reference Tyrka, Ridout and Parade2016). This stress-induced dysregulation of the HPA axis has been implicated in the pathogenesis of stress-related psychiatric disorders (Barden, Reference Barden2004; Braquehais, Picouto, Casas, & Sher, Reference Braquehais, Picouto, Casas and Sher2012; Doom & Gunnar, Reference Doom and Gunnar2013).
Substantial research has been aimed at determining whether early life stress is associated with epigenetic changes to the promoter region of the GR gene, nuclear receptor subfamily 3, group C, member 1 (NR3C1), that may alter GR expression and HPA axis homeostasis and responses to stress (see Turecki & Meaney, Reference Turecki and Meaney2016; Tyrka, Ridout, et al., Reference Tyrka, Ridout and Parade2016, for reviews). Epigenetic modifications to the genome allow for altered gene expression but do not change the DNA sequence and thus permit elaboration of the genome beyond what is determined by DNA base coding. The most highly studied and best characterized epigenetic modification, DNA methylation (DNAm), usually involves a direct covalent, chemical modification of a cytosine base lying sequentially adjacent to a guanine base (thus a CpG dinucleotide; CG/GC). DNAm is associated with the silencing of gene transcription that appears to be mediated in one of two ways (Bird, Reference Bird2002). First, wide swaths of DNA can be methylated, and the shear density of methylation precludes transcription factor binding to DNA sites, thus silencing gene expression. Second, when methylation levels are low, as in the case of NR3C1, small changes result in redistributing the transcriptional landscape, affecting translational isoform production, and orchestrating the final proteomic landscape that results in gene silencing (see Leenen, Muller, & Turner, Reference Leenen, Muller and Turner2016; Meaney, Reference Meaney2010, for reviews).
Maltreatment and other early adversities have been linked in a majority of studies to increased DNAm within the NR3C1 gene, resulting in reduced gene expression and decreased hippocampal GR (e.g., McGowan et al., Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009; Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004), and increased cortisol reactivity (e.g., Oberlander et al., Reference Oberlander, Weinberg, Papsdorf, Grunau, Misri and Devlin2008; Stroud et al., Reference Stroud, Papandonatos, Rodriguez, McCallum, Salisbury, Phipps and Marsit2014). Further, methylation of NR3C1 in rodents has been shown to be responsive to many forms of early life stress, including prenatal stress (Lillycrop et al., Reference Lillycrop, Slater-Jefferies, Hanson, Godfrey, Jackson and Burdge2007; Szyf, Reference Szyf2013) and postnatal stress (Kundakovic, Lim, Gudsnuk, & Champagne, Reference Kundakovic, Lim, Gudsnuk and Champagne2013; Witzmann, Turner, Meriaux, Meijer, & Muller, Reference Witzmann, Turner, Meriaux, Meijer and Muller2012). For example, low levels of maternal care in rodents causes greater methylation of the promoter region of the hippocampal GR gene, which interferes with binding of the transcription factor nerve growth factor inducible protein A (NGFI-A), resulting in reduced NR3C1 gene expression (Kosten & Nielsen, Reference Kosten and Nielsen2014; Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004). These biological effects of early environmental stress have also been found in humans with studies finding prenatal (e.g., Oberlander et al., Reference Oberlander, Weinberg, Papsdorf, Grunau, Misri and Devlin2008; Ostlund et al., Reference Ostlund, Conradt, Crowell, Tyrka, Marsit and Lester2016) and early childhood (e.g., Romens, McDonald, Svaren, & Pollak, Reference Romens, McDonald, Svaren and Pollak2015; Tyrka et al., Reference Tyrka, Parade, Eslinger, Marsit, Lesseur, Armstrong and Seifer2015) adversity linked with increased NR3C1 methylation evidenced in early development (e.g., Conradt et al., Reference Conradt, Hawes, Guerin, Armstrong, Marsit, Tronick and Lester2016; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016) and long-lasting into adulthood (e.g., Tyrka, Price, Marsit, Walters, & Carpenter, Reference Tyrka, Price, Marsit, Walters and Carpenter2012). In addition, these associations have been found across a number of human cell and tissue types including DNA from postmortem hippocampal tissue (e.g., Labonte et al., Reference Labonte, Yerko, Gross, Mechawar, Meaney, Szyf and Turecki2012; McGowan et al., Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009), placenta (e.g., Conradt et al., Reference Conradt, Lester, Appleton, Armstrong and Marsit2013; Kertes et al., Reference Kertes, Kamin, Hughes, Rodney, Bhatt and Mulligan2016), umbilical cord blood (e.g., Hompes et al., Reference Hompes, Izzi, Gellens, Morreels, Fieuws, Pexsters and Claes2013; Mulligan, D'Errico, Stees, & Hughes, Reference Mulligan, D'Errico, Stees and Hughes2012), peripheral blood (e.g., Romens et al., Reference Romens, McDonald, Svaren and Pollak2015; van der Knaap et al., Reference van der Knaap, Riese, Hudziak, Verbiest, Verhulst, Oldehinkel and van Oort2014), and saliva (e.g., Melas et al., Reference Melas, Wei, Wong, Sjoholm, Aberg, Mill and Lavbratt2013; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016). In sum, the association between early adversity and methylation of NR3C1 has been consistently observed across species and tissues (Vinkers et al., Reference Vinkers, Kalafateli, Rutten, Kas, Kaminsky, Turner and Boks2015) with replication of findings across multiple laboratories (Palma-Gudiel, Cordova-Palmera, Leza, & Fananas, Reference Palma-Gudiel, Cordova-Palomera, Leza and Fananas2015; Turecki & Meaney, Reference Turecki and Meaney2016).
Emerging evidence also suggests that methylation of NR3C1 is associated with behavior and psychopathology (e.g., Conradt et al., Reference Conradt, Lester, Appleton, Armstrong and Marsit2013; Kosten, Huang, & Nielsen, Reference Kosten, Huang and Nielsen2014; Na et al., Reference Na, Chang, Won, Han, Choi, Tae and Ham2014; Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004) as well as potentially functioning as a mediator for the link between early adversity and increased internalizing problems (Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016). However, findings as to the direction of NR3C1 methylation effects on behavior and psychopathology have been mixed with some studies finding positive associations with internalizing problems (e.g., Dadds, Moul, Hawes, Mendoza Diaz, & Brennan, Reference Dadds, Moul, Hawes, Mendoza Diaz and Brennan2015; van der Knaap, van Oort, Verhulst, Oldehinkel, & Riese, Reference van der Knaap, van Oort, Verhulst, Oldehinkel and Riese2015), while others have found negative associations with internalizing (Tyrka, Parade, et al., Reference Tyrka, Parade, Welch, Ridout, Price, Marsit and Carpenter2016) and posttraumatic stress (Labonte, Azoulay, Yerko, Turecki, & Brunet, Reference Labonte, Azoulay, Yerko, Turecki and Brunet2014; Vukojevic et al., Reference Vukojevic, Kolassa, Fastenrath, Gschwind, Spalek, Milnik and de Quervaine2014; Yehuda et al., Reference Yehuda, Flory, Bierer, Henn-Haase, Lehrner, Desarnaud and Meaney2015). One explanation for differential effects may be a developmental progression from dysregulation in early childhood to psychopathology in adulthood. Specifically, early acute stress-related hypermethylation results in emotional and behavioral dysregulation in early childhood (Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016) but then chronic environmental adversity throughout development may result in hypomethylation (Tyrka, Parade, et al., Reference Tyrka, Parade, Welch, Ridout, Price, Marsit and Carpenter2016) as an adaptation to repeated stress exposures that, for some, results in sustained psychopathology in adulthood (Labonte et al., Reference Labonte, Azoulay, Yerko, Turecki and Brunet2014; Vukojevic et al., Reference Vukojevic, Kolassa, Fastenrath, Gschwind, Spalek, Milnik and de Quervaine2014).
There has been substantial progress in behavioral epigenetics over the last decade (Lester, Conradt, & Marsit, Reference Lester, Conradt and Marsit2016). A strong emphasis of behavioral epigenetics research has been on understanding how epigenetic processes play a vital, potentially explanatory, role in the biological embedding of early environmental adversity and the genesis of adaptive and maladaptive development (Boyce & Kobor, Reference Boyce and Kobor2015). A major limitation, however, is that a basic understanding of how methylation of glucocorticoid signaling genes change over time has yet to be achieved, with little understanding of the role of risk factors that impact these developmental trajectories (Tyrka, Ridout, et al., Reference Tyrka, Ridout and Parade2016; Vinkers et al., Reference Vinkers, Kalafateli, Rutten, Kas, Kaminsky, Turner and Boks2015). Prior work examining the effect of early stress on methylation NR3C1 has been exclusively cross-sectional, greatly constraining our understanding of epigenetic processes over time. In addition, the extant cross-sectional data does not allow for examination of how maladaptive patterns of stress-related methylation changes are altered or maintained over time in association with adversity and the development of psychopathology. Such questions are crucial to understanding pathophysiology and informing the next generation of prevention and intervention efforts.
The primary aim of the current study was to begin to address the lack of longitudinal designs in environmental epigenetic research. We utilized a prospective longitudinal design with repeated assessments (two waves approximately 6 months apart) of methylation of saliva NR3C1 to understand how child maltreatment and other environmental adversity (e.g., stressor associated with living in poverty or other traumatic events) contributes to within-child change in NR3C1 methylation over time. Specifically, using a sample of preschoolers, approximately half of whom have documented maltreatment within 6 months of the baseline assessment, we examined the effect of maltreatment and other environmental adversity on short-term change in methylation to begin to characterize within-child epigenetic changes during this sensitive developmental period. We assessed DNA methylation of the well-characterized exon 1F, and extended this research to exon 1D, which also has support from cross-sectional studies for a role in these processes (Hompes et al., Reference Hompes, Izzi, Gellens, Morreels, Fieuws, Pexsters and Claes2013; Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016; van der Knaap et al., Reference van der Knaap, Riese, Hudziak, Verbiest, Verhulst, Oldehinkel and van Oort2014; Weder et al., Reference Weder, Zhang, Jensen, Yang, Simen, Jackowski and O'Loughlin2014).
Method
Participants
Two hundred sixty families residing in the Northeast enrolled in this study. One child from each family was included in the study. Children ranged in age from 3 to 5 years (M = 4.2, SD = 0.74), and 52.3% were female. The sample was racially and ethnically diverse (45.6% Hispanic, 27.7% White non-Hispanic, 16.3% Black, 21.9% biracial, and 2.7% other races). Nearly all children were living in poverty (91%). Most caregivers (95%) were biological mothers, 60% of caregivers had less than or equal to a high school diploma, 53.3% were single parents, 55% of caregivers were unemployed, and 86.5% of the families qualified for public assistance. Approximately half of the children (53.3%) had substantiated cases of moderate to severe child maltreatment within the past 6 months as described below. Baseline results from a subsample (n = 184; 70.7%) of the present sample were reported previously elsewhere (Tyrka et al., Reference Tyrka, Parade, Eslinger, Marsit, Lesseur, Armstrong and Seifer2015).
Procedure
Families with a maltreated child were identified from the local child welfare agency and an emergency maltreatment assessment service via record review. Families of children with no indicated case of maltreatment within the past 6 months were recruited at a pediatric medical clinic during a well-child visit as well as at childcare centers. Based on review of available medical records and parent report, children with a chronic illness, medication use, obesity, and failure-to-thrive were excluded. Those with acute illness or medication use were included no less than 2 weeks following resolution of illness and medication use. Families completed a baseline set of assessments at the time of initial study enrollment and a follow-up set of assessments 6 months following enrollment (M = 6.42 months, SD = 0.67 months). At each wave of assessment, families completed two home visits and questionnaires between the visits. The current report focuses on data from the first home visit during the baseline assessment during which caregivers completed interviews on child stress exposure and a baseline saliva sample for DNA isolation was collected from the children, as well as the first home visit during the follow-up assessment during which caregivers completed an interview on service utilization and a follow-up saliva sample for DNA isolation was collected from children.
Measures
Child maltreatment status
All families consented to examination of child welfare records to determine maltreatment status. Trained research staff coded the records using the System for Coding Subtype and Severity of Maltreatment in Child Protective Records (Barnett, Manly, & Cicchetti, Reference Barnett, Manly, Cicchetti, Cicchetti and Toth1993). Five maltreatment subtypes and severity scores ranging from 1 (least severe) to 5 (most severe) were derived. Children with a case of moderate to severe levels of maltreatment (score of 3–5) within the last 6 months were considered as part of the maltreated group (n = 134). Within the maltreated group, 41% of the sample had substantiated cases of multiple subtypes of maltreatment, 5 children had substantiated cases of physical abuse, 32 sexual abuse, 16 physical neglect/failure to provide, 36 physical neglect/lack of supervision, and 85 emotional maltreatment. The comparison group included children who had never had a substantiated case of maltreatment.
Additional adversity markers
Contextual stress interview
Caregivers completed a semistructured interview developed in our laboratory to assess contextual stressors experienced in the child's lifetime. Categories were death of a caregiver, separation from a caregiver, housing instability, inadequate food or clothing, and other stressful events, which included witnessing neighborhood violence or parental arrest. Interviews were conducted and scored by trained clinical social workers and PhD-level psychologists. The project coordinator reviewed each interview to ensure compliance with the scoring protocol. Each domain was scored positive if at least one episode occurred, and domains were summed to determine the number of contextual stressor categories the child experienced in his or her lifetime. Possible scores ranged from 0 (no stressors) to 5 (stressors in all five domains). In the current sample, the number of stressor categories ranged from 0 to 5 (M = 1.34, SD = 1.20).
Traumatic life events
The Diagnostic Infant and Preschool Assessment (Scheeringa & Haslett, Reference Scheeringa and Haslett2010) interview was conducted with caregivers to assess child experiences of traumatic life events. Interviews were conducted by trained clinical social workers and PhD-level psychologists, reviewed in a group supervision format, and scored based upon group consensus. Traumatic events in each domain were dichotomized (no trauma vs. ≥1 trauma), then summed to create a scale for number of types of traumas experienced in the child's lifetime. Physical and sexual abuse were not included because they were assessed as maltreatment (above). Possible scores range from 0 to 9 and in the present sample a mean of 1.01 and SD of 1.05 were observed.
Other adversity composite
Confirmatory factor analysis using Mplus was conducted to estimate a single factor latent adversity variable based on model fit and theoretical interpretability with the following indicators: the number of lifetime contextual stressors, the number of traumatic life events, single-parent status (coded as 0 or 1), and the child's lifetime experiences of homelessness (8%). Single-parent status and experiences of homelessness were collected as part of the demographic questionnaire. A robust weighted least squares estimator using a diagonal weight matrix was used because it provides the best option for modeling categorical or ordered data. Model fit was excellent, χ2 (2, N = 260) = 1.44, p > .15, root mean square error of approximation = 0.00, 95% confidence interval (CI) [0.00 –0.11], comparative fit index = 1.0, weighted root mean square residual = 0.245, and standardized factor loadings were as follows: homelessness (.71), single-parent status (.44), lifetime stressors (.89), and the number of traumatic life events (.62). The factor score of the latent adversity variable was saved and used as observed in primary analyses.
NR3C1 methylation
Saliva samples were obtained using the Oragene DISCOVER kits (OGR-575) for Assisted Collections (DNA Genotek, Kanata, Ontario, Canada) at baseline and follow-up visits, and DNA was isolated following the manufacturer's instructions. Sodium bisulfite modification was performed with 500 ng of DNA using the EZ DNA methylation Kit (Zymo Research, Irvine, CA). For DNA methylation detection, bisulfite pyrosequencing was employed in two locations within the NR3C1 region: promoter of exon 1D and promoter of exon 1F (three assays; Figure 1) at the two time points. CpG numbering was taken from Palma-Gudiel et al. (Reference Palma-Gudiel, Cordova-Palomera, Leza and Fananas2015) to improve comparability of results to past studies. PyroMark Assay Design software version 2.1.15 (Qiagen) was used to design the pyrosequencing assays. Amplification polymerase chain reactions and sequencing primers (Integrated DNA Technologies, Inc., Coralville, IA) and the genomic locations of the assays are provided in Table 1. The PyroMark polymerase chain reaction kit and forward and reverse primers were used to amplify specific regions of the NR3C1 promoter. Four forward pyrosequencing assays covering a total of 27 CpG loci were performed in triplicate using the PyroMark MD (Qiagen). The percentage DNA methylation at each CpG locus was quantified with the PyroMark CpG software, version 1.11 (Qiagen). All procedures were performed following manufacturer's protocols. The percentage of alleles that were methylated was used in statistical analyses.
(Color online) Glucocorticoid receptor gene NR3C1 region: promoter of exon 1D and promoter of exon 1F. CpG numbering taken from Palma-Gudiel et al. (Reference Palma-Gudiel, Cordova-Palomera, Leza and Fananas2015). Boxes around CpG site numbers represent nerve growth factor inducible protein A (NGFI-A) transcription factor binding sites according to McGowan et al. (Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009).
Pyrosequencing primers and assay sequences
a Primer sequences are given in the 5′ to 3′ direction.
Covariates
Child age in months, child sex, and the number of days between baseline and the follow-up assessment were all examined as predictors of NR3C1 methylation at baseline and change in methylation between baseline and the 6-month follow-up. If significant associations emerged, those variables were included as covariates in primary analyses.
Modeling ancestry differences using principal component analysis (PCA)
Allele frequency differences due to systematic ancestry differences could cause spurious associations. Thus, we used PCA to model ancestry differences in the current study using genome-wide single nucleotide polymorphism (SNP) markers. Genotyping was conducted using the Illumina Infinium PsychArray-24 beadchip (over 588,000 autosomal SNPs). Genotypes were called using GenomeStudio software V2011.1 and genotyping module version 1.8.4 (Illumina), and cleaned using standard quality control procedures. We first conducted linkage disequilibrium-based pruning to lessen overweight of principle components (PCs) of genetic variation by the contribution of correlated SNPs, and followed by PCA using PLINK (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender and Sham2007). The first two PCs obtained using the PLINK were used for controlling the potential population stratification (Price et al., Reference Price, Patterson, Plenge, Weinblatt, Shadick and Reich2006).
Data analytic plan
Prior to analyses, outliers, defined as values more than 3 SD from the mean, were Winsorized by setting them to the next highest value within 3 SD. Mplus 7.31 software (Muthén & Muthén, Reference Muthén and Muthén1998–2012) was used to conduct all analyses. To assess whether maltreatment and the other adversity composite predict baseline level and within-child changes over time in DNA methylation, we utilized latent change score (LCS) models. Of relevance to our main research questions, the flexibility of LCS allows for a simultaneous analysis of the predictors of individual differences in the initial level of and subsequent changes in key variables similar to a latent growth curve model but with two time points (McArdle, Reference McArdle2009). Primary LCS models examined maltreatment and the adversity composite simultaneously as predictors of baseline level and change in NR3C1 methylation. See Figure 2 for a depiction of the primary LCS model.
Primary latent change score models. *Freely estimated parameter. PC1-2, DNA-based principal components used to account for potential population stratification.
Missing data for methylation of single CpG sites was a maximum of 4.2% at baseline and 30.4% at 6 months with methylation data available for 86.5% of the sample at follow-up. Data were missing completely at random (Little's missing completely at random test, p > .15), and full information maximum likelihood estimation techniques were used for inclusion of all available data. To account for nonnormality, maximum likelihood estimation with robust standard errors was used. The following fit statistics were employed to evaluate model fit: χ2 (p > .05 excellent), comparative fit index (>0.90 acceptable, >0.95 excellent), root mean square error of approximation (<0.08 acceptable, <0.05 excellent), and standardized root mean square residual (<0.08 acceptable, <0.05 excellent; Hu & Bentler, Reference Hu and Bentler1999).
Results
Preliminary analyses
Covariates
In individual predictor LCS models child sex, child age, and the number of days between baseline and follow-up assessments were examined as predictors of baseline level and change in mean methylation of NR3C1 exon 1D and 1F. None of the three potential covariates were significantly related to baseline or change in methylation of 1D or 1F. Thus, these variables were not considered further in analyses, and only the top two genetic PCs were included as covariates to account for potential population stratification in primary LCS analyses.
CpG correlations
Correlations between CpG sites within 1D were mostly small to medium in magnitude, similar at baseline and follow-up, with no clear clustering within the 1D CpG sites. For correlations between CpG sites within 1F, some CpG sites were more closely associated, including modest correlations between immediately adjacent CpG sites from 27 to 32 (e.g., CpG 28 and 29 r = .35 whereas CpG 28 and 30 r = .16), and large correlations between putative NGFI-A binding region CpG sites 30–32 (average r = .66), and between CpG sites 40–47 (average r = .62). In other words, methylation at some CpG sites across the 1F promoter region fluctuate in a coordinated fashion whereas others are more loosely associated.
Primary analyses
We conducted primary analyses in two stages. First, LCS models were conducted with mean methylation across exons 1D and 1F. Second, based on the correlation results above and CFA models that did not support a single-factor structure for NR3C1 exon 1F (poor model fit for each model, complete results available upon request), LCS analyses were conducted with CpG bundles for 1F (i.e., average methylation across several sites) based on the literature and the above results: CpG sites 27–29 given higher than average interrelations (mean r = .42), CpGs 30–32 and 37–38 as they are known NGFI-A binding sites (McGowan et al., Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009), and 40–47 given higher than average interrelations (mean r = .59).
Mean methylation of NR3C1 exons 1D and 1F
An unconditional LCS model was estimated for each exon to characterized unconditional change over time. Across LCS models for exon 1D and 1F, the covariance of intercept (i.e., baseline) and change factors were significant and negative suggesting that preschoolers who had higher methylation values at baseline tended to decrease more rapidly across 6 months for exon 1D, β = –0.52, 95% CI [–0.63, –0.42], and 1F, β = –0.55, 95% CI [–0.65, –0.45]. In addition, mean rate of change and variances of intercept and change factors for exons 1D and 1F significantly differed from zero, all ps < .01, indicating potentially important individual variability in both starting point and change over time.
A summary of primary conditional LCS model results are presented in Table 2, including unstandardized path coefficients and 95% CI. Model fit for the LCS model for exon 1D and 1F was excellent. In regard to the mean 1D methylation model (Model 1, Table 2), maltreatment status, but not other adversity, positively predicted baseline level and negatively predicted change in methylation, such that maltreated preschoolers evidenced higher baseline levels of exon 1D methylation but then significant decreases in 1D methylation between baseline and follow-up assessments. When the LCS model was reversed to examine follow-up methylation levels as the intercept, maltreated preschoolers evidenced marginally (p < .10) lower methylation levels at follow-up. The results for predicting change in 1D were consistent whether maltreatment and other adversity were examined simultaneously or separately as predictors.
Summary of latent change score analyses for CpG bundles
Note: 1F NGFI-A, transcription factor binding sites as identified by McGowan et al. (Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009): 30–32 and 37–38. Bold cells represent p < .10. PC, principal component; RMSEA, root mean square error of approximation; CFI, comparative fit index; SRMR, standardized root mean square residual.
a DNA-based principal components used to account for potential population stratification. Predicting follow-up paths are from reverse latent change score models where follow-up functions as the intercept. The results from these models had nearly identical model fit.
†p < .10. *p < .05.
Results for the mean exon 1F methylation model (Model 2, Table 2) differed from those of 1D, such that other adversity, but not maltreatment status, positively predicted baseline methylation levels and neither predictor was significantly related to change in 1F methylation. Results for predicting change in 1F were consistent whether maltreatment and other adversity were examined simultaneously or separately as predictors.
NR3C1 exon 1F CpG bundles
LCS model results for exon 1F bundles are also presented in Table 2. Model fit for all models was excellent and comparable to model fit statistics reported above. For the CpG 27–29 bundle (Model 3, Table 2), other adversity was positively related to baseline levels of methylation, such that children who faced more adversity, including contextual stress and trauma exposure, evidenced higher levels of methylation in this region. Further, maltreatment status was negatively related to change in methylation for the 27–29 region, such that maltreated preschoolers evidenced decreases in methylation. Once the LCS model was reversed, maltreated preschoolers evidenced significantly lower methylation levels at follow-up, b = –0.20, 95% CI [–0.35, –0.06]. In other words, for this region, maltreated and comparison children did not significantly differ in methylation level at baseline, but maltreated children had significant decreases in methylation over time and ended with lower levels of methylation compared to comparison children. However, though maltreatment status was not a statistically significant predictor of baseline methylation for exon 1F 27–29, the effect size was nearly identical to that of 1D, where the difference was statistically significant. Finally, for the NGFI-A binding site bundle (Model 4, Table 2) and the CpG 40–47 bundle (Model 5, Table 2), neither maltreatment nor adversity significantly predicted baseline or change in methylation (though the effect of adversity on baseline methylation levels of 40–47 approached significance).
Post hoc analyses
Following primary LCS analyses, post hoc analyses were conducted using individual CpG sites to examine whether specific sites drove significant effects and because this study is the first to examine prospective short-term within-child change in NR3C1 methylation.
LCS analysis of individual CpG sites
Given the number of models (27) for LCS analyses with individual CpG sites, p values are not interpreted in favor of examining the pattern of results across CpG sites via effect sizes and 95% bootstrapped CI. A summary of LCS model results for individual CpG sites is presented in Table 3. Each row of Table 3 represents a different model based on an individual CpG site. The maltreatment and adversity columns present results of maltreatment or adversity predicting baseline levels and change in the methylation outcome for each model. Model fit for all models was excellent and comparable to model fit statistics reported above.
Summary of latent change score analyses for individual CpG sites
Note: Bold values represent 95% confidence intervals (CI) that either did not include zero or had one end approximately at zero but an effect size similar to other nearby CpG sites. Δ, latent change. Each row represents a separate model with both maltreatment and other adversity as simultaneous predictors.
Overall, maltreated children evidenced higher baseline levels of methylation across most exon 1D CpG sites but few 1F CpG sites. Further, higher levels of adversity were related to higher baseline levels of methylation across much of the 1F, but not 1D, promoter region, though most of this effect was outside of known NGFI-A binding sites (comparable effect sizes for CpG 27–29 and 40–47). With regard to prediction of change over time, maltreated children evidenced decreases in methylation relative to comparison children across most 1D CpG sites (CpG 6 as the exception), but for 1F this effect was only observed consistently at CpG sites 27–29 with the largest difference observed for 1F CpG 27, which was nearly twice as large as the next largest difference. Finally, adversity rarely predicted change in methylation with wide confidence intervals observed for nearly all CpG sites. Intriguingly, for CpG sites 5, 8, and 9 of exon 1D, maltreatment status predicted decreases in methylation over time whereas other adversity predicted increases in methylation over time.
Post hoc CpG site bundle
To better characterize the effect of maltreatment on baseline level and change in NRC31 methylation, we combined CpG sites that evidenced consistency across analyses (exon 1D CpGs 3–5, 7–10 and exon 1F CpGs 27–29). The results of this LCS model served as the basis of Figure 3, which displays mean methylation at baseline and follow-up across these regions. As is shown in Figure 3, across these sites, maltreated children evidence higher baseline levels of methylation, β = 0.15, 95% CI [0.01, 0.29], significant decreases in methylation over time, β = –0.23, 95% CI [–0.38, –0.08], and then at follow-up, lower levels of methylation, β = –0.19, 95% CI [0.35, –0.02], relative to nonmaltreated preschoolers (whose methylation levels remained stable) after accounting for the effects of other adversity and potential population stratification.
Depiction of glucocorticoid receptor gene NR3C1 change over time separately by maltreatment and comparison children. Average methylation across multiple CpG sites (exon 1D CpGs 3–5, 7–10 and exon 1F CpGs 27–29) is presented separately for maltreated and comparison children. *p < .05, ns, p > .05. The graph represents multiple analyses to better understand study results. Maltreated children, n = 134, comparison group, n = 126.
Discussion
Results from the current study suggest that child maltreatment is associated with higher initial levels of NR3C1 promoter methylation within 6 months of documented maltreatment as well as the rate of within-child change in NR3C1 methylation approximately 1 year after documented maltreatment. The primary aim of the current study was to begin to address the lack of longitudinal methylation designs in environmental epigenetic research by utilizing a short-term longitudinal model to examine predictors of change in methylation of saliva DNA NR3C1 in a sample of children. While others have examined change in NR3C1 DNAm after nonrandomized psychological intervention with small samples (Roberts et al., Reference Roberts, Keers, Lester, Coleman, Breen, Arendt and Havik2015; Yehuda et al., Reference Yehuda, Daskalakis, Desarnaud, Makotkine, Lehrner, Koch and Bierer2013), the current study is the first to examine early environmental adversity as a predictor of NR3C1 DNAm change. Further, though several studies have examined the influence of child maltreatment and other early environmental adversity on later NR3C1 DNAm in adulthood (e.g., Tyrka et al., Reference Tyrka, Price, Marsit, Walters and Carpenter2012; Tyrka, Parade, et al., Reference Tyrka, Parade, Welch, Ridout, Price, Marsit and Carpenter2016), the current study is the first to explicitly model NR3C1 promoter methylation change in early childhood proximal to a documented incidence of child maltreatment or current adversity inherent within environments of poverty. Establishing trajectories of short-term change following child maltreatment and temporal precedence before the onset of significant psychopathology is crucial to drawing inferences about the development of psychopathology (Jones, Moore, & Kobor, Reference Jones, Moore and Koborin press).
In the current study, baseline associations between child maltreatment status and other adversities with NR3C1 promoter methylation were largely consistent with those previously reported using a subsample of participants in the current study (Tyrka et al., Reference Tyrka, Parade, Eslinger, Marsit, Lesseur, Armstrong and Seifer2015). Specifically, preschoolers with documented maltreatment evidenced higher mean baseline methylation at NR3C1 region 1D, but not 1F, relative to nonmaltreated but demographically similar comparison children. In addition, other early life adversity (e.g., contextual stress associated with poverty and trauma exposure), but not maltreatment, was associated with higher baseline mean methylation at region 1F. Results from the current study extend these previous results by examining child maltreatment and other adversity as predictors of change in methylation at these same regions while accounting for PCs used to adjust for potential population stratification. Findings indicated that though maltreated children evidenced higher baseline levels of methylation, they also had significant decreases in methylation over time, and at the 6-month follow-up their methylation levels were lower than those of nonmaltreated preschoolers (whose methylation levels remained stable over time). These results were observed across most of the 1D CpG sites examined and one region of 1F (CpGs 27 and 29), but limited support was found at other 1F sites including the NGFI-A transcription factor binding sites identified by McGowan et al. (Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009).
Baseline associations found in the current study were consistent with a large body of research examining the influence of early environmental adversity on NR3C1 methylation that has been replicated across species and cell populations as diverse as peripheral blood cells, salivary DNA (which is primarily of leukocyte origin), and central nervous system-derived cells (Turecki & Meaney, Reference Turecki and Meaney2016; Tyrka, Ridout, et al., Reference Tyrka, Ridout and Parade2016; Vinkers et al., Reference Vinkers, Kalafateli, Rutten, Kas, Kaminsky, Turner and Boks2015). The baseline effect sizes in the current study were small to medium in magnitude, which is consistent with effect sizes found in other studies. One of the most pressing problems in environmental epigenetics is that we currently do not know whether these small effect sizes constitute biologically meaningful differences (Vinkers et al., Reference Vinkers, Kalafateli, Rutten, Kas, Kaminsky, Turner and Boks2015). Nevertheless, there is evidence that small effect sizes have significant downstream effects on gene expression (Breton et al., Reference Breton, Marsit, Faustman, Nadeau, Goodrich, Dolinoy and Yousefi2017), and, specifically for NR3C1, significant effects of methylation on postmortem human GR expression are found in the hippocampus (Labonte et al., Reference Labonte, Yerko, Gross, Mechawar, Meaney, Szyf and Turecki2012; McGowan et al., Reference McGowan, Sasaki, D'Alessio, Dymov, Labonte, Szyf and Menaey2009).
The current study is the first to test this association longitudinally, and thus, results are in need of replication and extension to longer time intervals. The observation of dynamic stress-related changes in DNAm during a sensitive developmental period highlights the importance of longitudinal designs for environmental epigenetic research. Despite the lack of extant literature for interpretation of this longitudinal association, results from the current study are consistent with previous findings from two longitudinal twin studies that examined intraindividual longitudinal change in DNAm. First, though not examining NR3C1, Wong et al. (Reference Wong, Caspi, Williams, Craig, Houts, Ambler and Mill2010) observed changes in DNAm between ages 5 and 10 among monozygotic twin pairs at the promoter/regulatory regions of the dopamine receptor D4 (DRD4), serotonin transporter (SERT), and monoamine oxidase A (MAOA) genes that were primarily attributable to environmental, and not heritable, influences. Second, Levesque et al. (Reference Lévesque, Casey, Szyf, Ismaylova, Ly, Verner and Booij2014) examined the stability of genome-wide DNA methylation patterns over 3–6 months in a subsample of 8 adolescent monozygotic twins and found that NR3C1 methylation was part of a highly variable statelike gene network that may be highly responsive to changes in the environment. The results of the current study support and extend those of these two previous studies by finding that child maltreatment predicts changes in NR3C1 promoter methylation in early childhood.
However, the directionality of stress-related DNAm change was unexpected, and DNAm differences between maltreated and comparison children reversed between assessments for some CpG sites in the current study (see Figure 3). Previously reported associations between early adversity and NR3C1 promoter methylation have not been universal in direction. Though the majority of published studies find higher levels of methylation of NR3C1 with early environmental stress (Turecki & Meaney, Reference Turecki and Meaney2016), some studies have found the opposite (Daskalakis & Yehuda, Reference Daskalakis and Yehuda2014; Turecki & Meaney, Reference Turecki and Meaney2016; Tyrka, Parade, et al., Reference Turecki and Meaney2016). One hypothesis for the dynamic stress-related DNAm changes depicted in Figure 3 is that early acute stress-related hypermethylation results in emotional and behavioral dysregulation in early childhood (Parade et al., Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016) but then chronic or severe environmental adversity throughout development results in hypomethylation (Tyrka, Parade, et al., Reference Turecki and Meaney2016) as an adaption to repeated stress exposures that, for some, results in sustained psychopathology (Labonte et al., Reference Labonte, Azoulay, Yerko, Turecki and Brunet2014; Vukojevic et al., Reference Vukojevic, Kolassa, Fastenrath, Gschwind, Spalek, Milnik and de Quervaine2014). In line with this developmental progression of biological embedding hypothesis, the DNAm difference at baseline in the current study may represent an acute stress response associated with concurrent emotional and behavioral dysregulation as reported by Parade et al. (Reference Parade, Ridout, Seifer, Armstrong, Marsit, McWilliams and Tyrka2016) using a subsample (n = 174) of the current studies participants. After the initial abuse is stopped or reduced after state agency involvement, this initial hypermethylation and dysregulation may transition, for some, to hypomethylation.
The developmental progression hypothesis offered above is supported by results from the current study finding reductions in methylation over time as well as lower levels at follow-up, approximately 1 year after documented abuse, for maltreated children relative to nonmaltreated comparison children. In addition, this hypothesis is supported by research with adult populations consistently finding that hypomethylation is associated with posttraumatic stress disorder (e.g., Labonte et al., Reference Labonte, Azoulay, Yerko, Turecki and Brunet2014; Vukojevic et al., Reference Vukojevic, Kolassa, Fastenrath, Gschwind, Spalek, Milnik and de Quervaine2014; Yehuda et al., Reference Yehuda, Flory, Bierer, Henn-Haase, Lehrner, Desarnaud and Meaney2015). Furthermore, this hypothesis is also supported by recent research finding that within-subject decreasing DNAm levels over time at several genes, though not NR3C1, were related to the increasing levels of posttraumatic stress disorder symptoms over time among military veterans (Rutten et al., Reference Rutten, Vermetten, Vinkers, Ursini, Daskalakis, Pishva and Kenis2017). If the developmental progression hypothesis is correct, we would expect the current trajectory of change in NR3C1 promoter methylation to continue into middle childhood and adolescence with accompanying onset of posttraumatic stress disorder or severe mood dysregulation. Thus, further research utilizing longitudinal designs with follow-ups across developmental stages are needed to support or refute this hypothesis.
An alternative hypothesis for explaining the observed reduction in NR3C1 promoter methylation in the current study is that this change is an intervention-induced return to baseline after state agency involvement. The higher levels of NR3C1 methylation and emotional dysregulation seen in maltreated preschoolers might be reversed with early intervention. According to this hypothesis, elevated levels of DNAm would have persisted without intervention, and the observed reduction in DNAm is because of intervention in early childhood before stabilization of the epigenome later in development. However, this hypothesis was not supported by post hoc analyses in the current study that found no association between service utilization (assess via an interview at the follow-up assessment) and change in NR3C1 methylation among maltreated preschoolers nor that service utilization moderated the association between maltreatment and change in methylation. Yet, the services interview in the current study was limited, without indictors of quality, frequency, or intensity of services nor effectiveness of intervention services on treatment targets. Therefore, future longitudinal research would benefit from designs that more effectively parse the contextual factors that might promote or inhibit methylation following severe adversity.
There are several limitations of the current study that should be noted. It is possible that the longitudinal results observed in the present study could be accounted for by unmeasured confounders. Baseline associations in the current study are consistent with numerous other studies across species and tissues, which increase our confidence in these results. However, the longitudinal stress-related associations in the present study are the first such results reported, and as such, replication studies are necessary prior to integration into established conceptualizations of the effects of adversity on NR3C1 methylation. Specifically, one potentially important unmeasured confounder in the current study is cell composition. Because our young children were not able to provide saliva via passive drool, we used a standardized method that involves using a sponge to collect pooled saliva beneath the tongue and at the intersection of the cheek and gum, so that in addition to various types of leukocytes, it is likely to yield some epithelial cells. Cell type within a tissue has been found to be the second biggest contributor to DNAm variation after tissue type (Farré et al. Reference Farré, Jones, Meaney, Emberly, Turecki and Kobor2015). However, despite this cell type issue, Smith et al. (Reference Smith, Kilaru, Klengel, Mercer, Bradley, Conneely and Binder2015) found DNAm from saliva may be more similar to that of brain tissues on average than DNAm from blood. Further, the best method to correct for interindividual variability in cell type is still highly debated (McGregor, Labbe, & Greenwood, Reference McGregor, Labbe and Greenwood2017), which is further complicated by the lack of well-characterized reference profiles for young child saliva. Regardless, future environmental epigenetics studies will benefit from controlling for this important type of variation (Jones et al., in press).
An additional limitation of the current study is the lack of accompanying gene expression data. Although there is prior evidence that NR3C1 promoter DNAm is associated with gene expression (Turecki & Meaney, Reference Turecki and Meaney2016), future longitudinal epigenetics research will benefit from assessing whether the DNAm changes observed alter gene expression (Jones et al., in press).
The current study also had several significant strengths that should be noted. First, the current study utilized a diverse sample of preschoolers exposed to a range of adversities including documented maltreatment while also controlling for potential population stratification in all primary analyses. Second, the current study is the first to explicitly model change in NR3C1 methylation signifying an important step in this line of research. Third, we report CpG site-specific effects given their utility in guiding basic molecular research on the biological and functional effects of DNAm differences observed (Daskalakis & Yehuda, Reference Daskalakis and Yehuda2014). Fourth, the inconsistent numbering of NR3C1 CpG sites have contributed to difficulty comparing results across multiple laboratories (Turecki & Meaney, Reference Turecki and Meaney2016), so we employed a universal CpG number system developed by Palma-Gudiel et al. (Reference Palma-Gudiel, Cordova-Palomera, Leza and Fananas2015)
In summary, epigenetic marks, such as DNAm, appear to functionally mold genetically guided developmental plasticity in response to early environmental experiences, thus providing a molecular basis for the enduring effects of early adversity exposures via biological embedding (Jones et al., in press). The current study is the first to examine maltreatment and other adversities inherent within environments of poverty as predictors of intraindividual change in NR3C1 methylation. Results from the current study highlight the complex nature of stress-related epigenetic variation during early development. Future studies that model early environmental adversity-induced DNAm changes into later developmental stages will continue to add to our understanding of the molecular mechanisms underlying the development of psychopathology.
