Introduction
Parental sensitivity, i.e. the responsiveness to children's signals and communications, is an important predictor of developmental outcomes across the behavioral, emotional, and cognitive domains (Kok et al., Reference Kok, Linting, Bakermans-Kranenburg, van IJzendoorn, Jaddoe, Hofman and Tiemeier2013; Madigan et al., Reference Madigan, Prime, Graham, Rodrigues, Anderson, Khoury and Jenkins2019; Thomas, Letourneau, Campbell, Tomfohr-Madsen, & Giesbrecht, Reference Thomas, Letourneau, Campbell, Tomfohr-Madsen and Giesbrecht2017). Low sensitivity of primary caregivers – typically mothers – has been associated with a host of negative outcomes, including higher risk for child psychopathology (Haltigan, Roisman, & Fraley, Reference Haltigan, Roisman and Fraley2013; Kimbrel, Nelson-Gray, & Mitchell, Reference Kimbrel, Nelson-Gray and Mitchell2007), externalizing and internalizing problems (Kok et al., Reference Kok, Linting, Bakermans-Kranenburg, van IJzendoorn, Jaddoe, Hofman and Tiemeier2013; Rijlaarsdam et al., Reference Rijlaarsdam, Stevens, Jansen, Ringoot, Jaddoe, Hofman and Tiemeier2014), and lower cognitive abilities (Bernier, Carlson, Deschênes, & Matte-Gagné, Reference Bernier, Carlson, Deschênes and Matte-Gagné2012). This influence can be long-lasting, as shown by prospective human studies (Raby, Roisman, Fraley, & Simpson, Reference Raby, Roisman, Fraley and Simpson2015; Stams, Juffer, & van IJzendoorn, Reference Stams, Juffer and van IJzendoorn2002) and experimental work in animals (Meaney, Reference Meaney2001). Yet, the molecular mechanisms underlying the enduring effects of maternal care on neurodevelopmental and behavioral outcomes in humans remain unclear.
Previous studies have provided initial support for DNA methylation (DNAm) – an epigenetic modification regulating gene expression – as a mechanism of interest for these processes (Mulder, Rijlaarsdam, & Van IJzendoorn, Reference Mulder, Rijlaarsdam, Van IJzendoorn, K. and R.2017; Szyf, Reference Szyf2013; Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004). DNAm involves the addition of a methyl group to DNA base pairs, primarily to the 5-carbon of cytosine nucleotides, resulting in 5-methylcytosine. DNAm is sensitive to both environmental and genetic influences (Ladd-Acosta & Fallin, Reference Ladd-Acosta and Fallin2016; Smith et al., Reference Smith, Kilaru, Kocak, Almli, Mercer, Ressler and Conneely2014; Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004), with the latter being evidenced by changes in the methylome due to DNA variation, named methylation quantitative trait loci (mQTLs) (Gaunt et al., Reference Gaunt, Shihab, Hemani, Min, Woodward, Lyttleton and Relton2016). Further, DNAm plays an essential role in healthy development and functioning by modulating the programming of wider biological systems (e.g. neural and immune functioning) and by coordinating key cellular processes (e.g. tissue differentiation) (Carey, Reference Carey2012). DNAm might thus represent a mechanism by which genetic and environmental factors, including the early caregiving environment, jointly and/or independently predict developmental outcomes (Ladd-Acosta & Fallin, Reference Ladd-Acosta and Fallin2016).
Most evidence of maternal care effects on DNAm comes from animal models. In a seminal study by Weaver et al. (Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004), low levels of maternal care in the first week of life – operationalized as the frequency of licking/grooming and arched-back nursing behaviors – altered hippocampal DNAm patterns in offspring at the glucocorticoid receptor (gr, also known as nr3c1), a key regulator of stress response (Geer et al., Reference Geer, Marchler-Bauer, Geer, Han, He, He and Bryant2010). Importantly, these epigenetic changes were long-lasting, but could be reversed via cross-fostering or chemical interventions, leading to a normalization of physiological and behavioral responses to stress (Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004, Reference Weaver, Champagne, Brown, Dymov, Sharma, Meaney and Szyf2005). These findings generated widespread interest, as they indicated (i) a causal role of maternal care on offspring's epigenetic dysregulation and downstream phenotypes, independent of genetic liability, and (ii) the possibility of influencing developmental trajectories through environmental interventions, mediated by DNAm. Since this initial work, other studies have replicated the effects of maternal care on gr methylation in rodents (Turecki & Meaney, Reference Turecki and Meaney2016) and extended findings to demonstrate DNAm changes in other tissues and genes (Beery, McEwen, MacIsaac, Francis, & Kobor, Reference Beery, McEwen, MacIsaac, Francis and Kobor2016; Blaze et al., Reference Blaze, Asok, Borrelli, Tulbert, Bollinger, Ronca and Roth2017; Doherty, Forster, & Roth, Reference Doherty, Forster and Roth2016) [e.g. brain-derived neurotrophic factor (bdnf) and oxytocin receptor (oxtr)] as well as in other species such as rhesus macaques (Provençal et al., Reference Provençal, Suderman, Guillemin, Massart, Ruggiero, Wang and Szyf2012).
Although rodents and primates widely differ from humans in their caregiving, a number of similarities in maternal–infant relationships have been observed across mammalian species (Feldman, Reference Feldman2016; Knop, Joëls, & van der Veen, Reference Knop, Joëls and van der Veen2017). Parallels at the sensory, hormonal, behavioral, and brain circuit levels have been noted (Feldman, Reference Feldman2016; Glynn & Baram, Reference Glynn and Baram2019; Knop et al., Reference Knop, Joëls and van der Veen2017), including the touch-based behavior characterizing rodents, primates, and humans in the early caregiving and the involvement of the limbic network in maternal–infant relationships (Feldman, Reference Feldman2016). Guided by the animal literature, a growing number of studies have sought to determine the extent to which different forms of caregiving and adversities affect DNAm in humans.
Human studies have focused on different forms of adversities (Daskalakis & Yehuda, Reference Daskalakis and Yehuda2014) including poly-victimization (Marzi et al., Reference Marzi, Sugden, Arseneault, Belsky, Burrage, Corcoran and Caspi2018), and on extreme deviations in the early caregiving environment, such as maltreatment (Cecil et al., Reference Cecil, Smith, Walton, Mill, McCrory and Viding2016; Gouin et al., Reference Gouin, Zhou, Booij, Boivin, Côté, Hébert and Vitaro2017; Mehta et al., Reference Mehta, Klengel, Conneely, Smith, Altmann, Pace and Binder2013; Stenz et al., Reference Stenz, Prados, Courtet, Prada, Nicastro, Adouan and Perroud2016; Weder et al., Reference Weder, Zhang, Jensen, Yang, Simen, Jackowski and Kaufman2014), institutionalization (Naumova et al., Reference Naumova, Lee, Koposov, Szyf, Dozier and Grigorenko2012), and maternal psychopathology (Oberlander et al., Reference Oberlander, Weinberg, Papsdorf, Grunau, Misri and Devlin2008). Generally, literature focusing on the caregiving environment has provided preliminary support in line with animal findings, identifying, for example, similar increases in GR methylation in both postmortem hippocampal tissue and peripheral tissues from individuals exposed to childhood maltreatment or early-life stress (Turecki & Meaney, Reference Turecki and Meaney2016). Studies also indicate that epigenetic patterns associated with the caregiving environment extend beyond GR, implicating other genes related to, among other processes, neurodevelopment and stress, such as OXTR and BDNF. Moreover, by leveraging epigenome-wide DNAm, novel genes were identified (e.g. KCNQ2, miR124-3) in relation to maltreatment and child abuse in individuals with post-traumatic stress disorder (Mehta et al., Reference Mehta, Klengel, Conneely, Smith, Altmann, Pace and Binder2013), borderline personality disorder (Stenz et al., Reference Stenz, Prados, Courtet, Prada, Nicastro, Adouan and Perroud2016), and depression (Weder et al., Reference Weder, Zhang, Jensen, Yang, Simen, Jackowski and Kaufman2014).
While these results are promising and suggest a role of the caregiving environment in the human methylome, the current evidence in humans is limited in a number of key ways. First, since research has mostly focused on extreme deviations in the caregiving environment in selected samples, little is known about how typical variation in maternal sensitivity associates with offspring DNAm in the general population. Second, while studies on extreme deviations in maternal care have leveraged epigenome-wide approaches, the literature on normative variation in maternal care has solely focused on candidate genes. This has impeded the identification of novel DNAm loci associated with maternal sensitivity, which might instead be detected with a hypothesis-free approach by performing an epigenome-wide association study (EWAS). Third, studies have typically relied on cross-sectional designs, in which the early caregiving environment is measured retrospectively via the use of questionnaires, raising doubts about the directionality of observed associations and about the validity of measurements, which may be prone to recall bias (Baldwin, Reuben, Newbury, & Danese, Reference Baldwin, Reuben, Newbury and Danese2019; Reuben et al., Reference Reuben, Moffitt, Caspi, Belsky, Harrington, Schroeder and Danese2016). Moreover, previous studies rarely investigated whether the identified associations may be confounded by genetic background shared between parents and offspring. The examination of the relationship between maternal care and DNAm might indeed capture intergenerational genetic transmission. Lastly, the influence on offspring DNAm of factors preceding postnatal maternal care, including the prenatal environment, remains unexplored.
To address these gaps, we firstly examined how typical variation in observed maternal sensitivity prospectively associates with epigenome-wide DNAm patterns in a general population of children. Secondly, with a series of follow-up analyses, we explored the extent to which associations reflected genetic influences as well as confounding by ‘baseline’ DNAm levels at birth, which precede exposure to postnatal maternal care and might constitute a biological indicator of the prenatal environment as well as of genetic effects on the methylome.
Materials and methods
Participants
The present research was embedded in the Generation R Study, a prospective population-based cohort study from fetal life onwards in Rotterdam, The Netherlands (Kooijman et al., Reference Kooijman, Kruithof, van Duijn, Duijts, Franco, van IJzendoorn and Jaddoe2016) (online Supplementary Information 1). Ethical approval was obtained from the Medical Ethics Committee of Erasmus MC, University Medical Center Rotterdam. For the purposes of this study, children within the Generation R Study with data on maternal sensitivity (at 3 and/or 4 years) and DNAm (at 6 years) were selected (N = 235). Since 5 sibling-pairs were present, we later excluded one sibling per pair (N = 230) to ensure genetic relatedness did not impact results.
Participant characteristics are shown in online Supplementary Table S1. Participants with data on both maternal sensitivity and DNAm (age 6) differed from participants invited to the age 6 assessment in gestational age at birth [M subsample = 40.3 weeks (s.d. = 1.4), M fullsample = 39.8 (s.d. = 1.9), t = 5.6, p = 6.50 × 10−08], but not other covariates.
Measures
Maternal sensitivity
Maternal sensitivity was assessed at ages 3 and 4 years through observations of mother–child interactions during teaching tasks too complex for the age of the child. These involved (i) building a tower and (ii) completing an etch-a-sketch drawing. Mother–child interactions were recorded and subsequently coded, according to the revised Erickson seven-point rating scales (Egeland, Erickson, Clemenhagen-Moon, Hiester, & Korfmacher, Reference Egeland, Erickson, Clemenhagen-Moon, Hiester and Korfmacher1990), based on two interdependent subscales: intrusiveness (IN) and supportive presence (SP), which together form the maternal sensitivity construct. Inter-coder reliability amounted to 0.81 at age 3 and 0.84 at age 4 (Kok et al., Reference Kok, Thijssen, Bakermans-Kranenburg, Jaddoe, Verhulst, White and Tiemeier2015).
Eight measures of maternal sensitivity (i.e. IN and SP scales × two tasks × two time-points) were available. IN scores were reversed, and both IN and SP scores were standardized. An overall maternal sensitivity score was calculated, for participants with data at age 3 and/or 4, by averaging such standardized measures (Cents et al., Reference Cents, Kok, Tiemeier, Lucassen, Székely, Bakermans-Kranenburg and Berg2014). This was done in line with previous literature (Kok et al., Reference Kok, Thijssen, Bakermans-Kranenburg, Jaddoe, Verhulst, White and Tiemeier2015), due to the stability of the maternal sensitivity scores between age 3 and 4 years (Kok et al., Reference Kok, Linting, Bakermans-Kranenburg, van IJzendoorn, Jaddoe, Hofman and Tiemeier2013), the temporality of these assessments, which both precede DNAm at age 6, and to maximize our sample size. Cronbach's α reliability of the obtained measure was acceptable (Cronbach's α = 0.70) (Cortina, Reference Cortina1993).
DNA methylation
DNAm in whole blood at age 6 was used for our epigenome-wide analyses. This was selected due to it being the closest DNAm assessment after maternal sensitivity observations (age 3 and 4 years), and to test the prospective association of maternal sensitivity with DNAm. Based on previous studies in animals, which found maternal care to have long-lasting influences on the methylome (Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004), we expected for maternal care effects to endure in early childhood.
To obtain DNAm data, DNA extraction and bisulfite conversion via the EZ-96 DNA Methylation kit (Shallow) (Zymo Research Corporation, Irvine, California, USA) were performed, and samples were processed with the Illumina Infinium HumanMethylation450 BeadChip (Infinium 450 K), which measures 485 577 CpGs. The incorporating control probe adjustment and reduction of global correlation pipeline (Lehne et al., Reference Lehne, Drong, Loh, Zhang, Scott, Tan and Chambers2015) was employed for the preparation and normalization of the data using R. Firstly, the minfi package (Aryee et al., Reference Aryee, Jaffe, Corrada-Bravo, Ladd-Acosta, Feinberg, Hansen and Irizarry2014) in R was used to read the idat files. Probes that had a detection p value above background (based on the sum of methylated and unmethylated intensity values) ≥1 × 10−16 were set to missing per array. Next, the intensity values were stratified by autosomal and non-autosomal probes and quantile normalized for each of the six probe-type categories separately: type II red/green, type I methylated red/green, and type I unmethylated red/green. For each probe, DNAm levels were indexed by β values (i.e. the ratio of methylated signal divided by the sum of the methylated and unmethylated signal [M/(M + U + 100)]). Quality control procedures were additionally performed (e.g. check for sex mismatch). Only arrays with a call rate above 95% per sample were considered for additional processing. DNAm data were winsorized (>3 s.d.) to reduce the influence of potential outliers. In total, we obtained information on 457 872 autosomal sites in 493 6-year-olds.
We additionally used DNAm data collected at birth in cord blood for a follow-up analysis. This was subject to the same pipeline as the DNAm data at age 6 and was also measured based on the Infinium 450 K BeadChip. Only CpGs identified as significant or within DNAm significant regions were selected for these analyses.
Covariates
All analyses were adjusted for a key set of covariates guided by previous literature (Birney, Smith, & Greally, Reference Birney, Smith and Greally2016; Breton et al., Reference Breton, Byun, Wenten, Pan, Yang and Gilliland2009; Liang & Cookson, Reference Liang and Cookson2014; Rakyan, Down, Balding, & Beck, Reference Rakyan, Down, Balding and Beck2011), including batch effects (plate number), sex, gestational age at birth, maternal smoking during pregnancy (never smoked, smoked until pregnancy known, continued during pregnancy), and estimated cell-type proportions (Houseman et al., Reference Houseman, Accomando, Koestler, Christensen, Marsit, Nelson and Kelsey2012) (online Supplementary Information 1). We additionally adjusted for two sets of covariates: (i) maternal education (highest level completed) as a proxy for socioeconomic status, and postnatal maternal psychopathology (Brief Symptom Inventory), and (iii) DNAm levels at birth (cord blood tissue), together with respective cell-type and batch effect adjustments (online Supplementary Information 1).
Statistical analyses
Analyses were performed in R (version 4.0.0) and are described in-depth in online Supplementary Information 1. A probe-level EWAS (multiple linear regression models) was run with the CpGassoc R package (Barfield, Kilaru, Smith, & Conneely, Reference Barfield, Kilaru, Smith and Conneely2012), to test for associations of maternal sensitivity with each DNAm site individually (Bonferroni epigenome-wide significance threshold: p < 1.09 × 10−07). To account for potential bias and inflation, the BACON R package (van Iterson, van Zwet, & Heijmans, Reference van Iterson, van Zwet and Heijmans2017) was used.
Moreover, to capture correlations across CpGs, reduce data dimensionality, and attenuate the multiple testing burden, a regional-level EWAS was performed by using the R package DMRff (Suderman et al., Reference Suderman, Staley, French, Arathimos, Simpkin and Tilling2018). This estimates correlations across nominally-significant probes within a 500 bp window (default setting) and combines the EWAS summary statistics of such neighboring CpGs to identify differentially methylated regions while accounting for multiple testing with a Bonferroni procedure in both gene regions and sub-regions (Suderman et al., Reference Suderman, Staley, French, Arathimos, Simpkin and Tilling2018).
A candidate gene look-up was also performed to maximize comparability with previously reported DNAm–maternal care associations. To date, DNAm levels of four genes have been associated with maternal care in humans (Bosmans, Young, & Hankin, Reference Bosmans, Young and Hankin2018; Conradt et al., Reference Conradt, Hawes, Guerin, Armstrong, Marsit, Tronick and Lester2016; Provenzi et al., Reference Provenzi, Fumagalli, Giorda, Morandi, Sirgiovanni, Pozzoli and Montirosso2017; Unternaehrer et al., Reference Unternaehrer, Meyer, Burkhardt, Dempster, Staehli, Theill and Meinlschmidt2015), by at least one study: GR, BDNF, the serotonin receptor (SLC6A4), and OXTR. We looked-up the EWAS results for probes located within these genes, as annotated in the HumanMethylation450 v1.2 Manifest File. Following previous studies (Cecil et al., Reference Cecil, Walton, Jaffee, O'Connor, Maughan, Relton and Barker2017; Marzi et al., Reference Marzi, Sugden, Arseneault, Belsky, Burrage, Corcoran and Caspi2018), gene-level Bonferroni correction was used as a significance threshold (i.e. p < 0.05/number of annotated probes).
To identify enriched biological pathways, we performed an in-house gene ontology (GO) analysis (Cecil et al., Reference Cecil, Walton, Jaffee, O'Connor, Maughan, Relton and Barker2017, Reference Cecil, Walton, Pingault, Provençal, Pappa, Vitaro and McCrory2018; Hannon et al., Reference Hannon, Dempster, Viana, Burrage, Smith, Macdonald and Mill2016) on sites with p < 0.001 in the probe-level EWAS, in line with previous literature (Cecil et al., Reference Cecil, Walton, Jaffee, O'Connor, Maughan, Relton and Barker2017, Reference Cecil, Walton, Pingault, Provençal, Pappa, Vitaro and McCrory2018; Mulder et al., Reference Mulder, Walton, Neumann, Houtepen, Felix, Bakermans-Kranenburg and Cecil2020; Roberts et al., Reference Roberts, Suderman, Zammit, Watkins, Hannon, Mill and Fisher2019). We performed p value adjustments based on default procedures (Hannon et al., Reference Hannon, Dempster, Viana, Burrage, Smith, Macdonald and Mill2016). Enriched pathways were confirmed by an independent GO approach from the missMethyl R package (Phipson, Maksimovic, & Oshlack, Reference Phipson, Maksimovic and Oshlack2016) (p < 0.05).
Finally, a series of follow-up analyses were run. Firstly, the influence of genetic factors on our top hits (i.e. Bonferroni-significant sites or sites within Bonferroni-significant DNAm regions) was assessed by drawing on an mQTL database (Gaunt et al., Reference Gaunt, Shihab, Hemani, Min, Woodward, Lyttleton and Relton2016) (www.mqtldb.org). We examined whether hits were associated with known mQTLs during childhood, based on the results from a genome-wide complex trait conditional analysis. Secondly, we explored the robustness of top hits to additional adjustments for (i) postnatal maternal education and maternal psychopathology (N = 223) and (ii) pre-exposure DNAm (N = 226). The latter was done to account for the effect of DNAm at birth on DNAm at age 6 and to capture potential pre-existing influences (e.g. intrauterine exposures) on DNAm in childhood. Spearman correlations between DNAm at birth and age 6 were also calculated, per CpG. Thirdly, based on a list of our CpG hits, the in-house GO analysis and missMethyl validation were run, with the same procedures as the main GO analysis specified above. Finally, to understand the relevance of our findings to the brain, which is linked to the caregiving environment (Kok et al., Reference Kok, Thijssen, Bakermans-Kranenburg, Jaddoe, Verhulst, White and Tiemeier2015; Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004), we looked-up brain–blood concordance values for our top hits using the BECon online tool (https://redgar598.shinyapps.io/BECon/) (Edgar, Jones, Meaney, Turecki, & Kobor, Reference Edgar, Jones, Meaney, Turecki and Kobor2017).
Results
Probe-level EWAS
Maternal sensitivity was not associated with any single CpGs at age 6, after genome-wide correction (p < 1.09 × 10−07) (Fig. 1, online Supplementary Table S2). BACON analysis revealed a normal λ (λ = 1.00), minimal bias (Bayesian estimate of bias = −0.002), and deflation in the test results – indicative of low power (Bayesian inflation factor = 0.925) (online Supplementary Fig. S1). Following BACON correction for deflation, one intergenic CpG reached genome-wide significance: cg25628898 (estimate = −0.008; s.e. = 0.002; p = 1.03 × 10−07) (online Supplementary Table S2).
Fig. 1. Manhattan plot of CpG sites associated with maternal sensitivity. Note. The Manhattan plot displays the log p values for each site tested in association with maternal sensitivity in the EWAS, across autosomal chromosomes. No genome-wide significant association was observed (p < 1.06 × 10−07).
Regional-level EWAS
With a regional-level EWAS, we identified 13 DNAm regions associated with maternal sensitivity (p < 1.09 × 10−07; α = 0.05) (Table 1, Fig. 2, online Supplementary Table S3), spanning 143 CpGs. The top three DNAm regions coincided with the ANKMY1, RNF39, and ZBTB22 and TAPBP genes (Table 1). The largest estimates were shown at regions encompassing COLEC11 and DOCK4. None of the CpGs within our significant regions was related to prenatal maternal smoking, based on previous research in neonates and children (Joubert et al., Reference Joubert, Felix, Yousefi, Bakulski, Just, Breton and London2016; Rzehak et al., Reference Rzehak, Saffery, Reischl, Covic, Wahl, Grote and Koletzko2016), suggesting adjustments in the EWAS accounted for its confounding role. When siblings (N = 230) were excluded, all but one region (annotated to RNF5P1, RNF5, AGPAT1) remained significantly associated with maternal sensitivity.
Fig. 2. Miami plot of DNA methylation regions associated with maternal sensitivity. Note. The Miami plot displays the log p values and estimates direction for each DNA methylation region tested in association with maternal sensitivity, across autosomal chromosomes. Thirteen regions were Bonferroni significant, three of which showed a positive relation with maternal sensitivity and 10 a negative one.
Table 1. DNA methylation regions significantly associated with maternal sensitivity from the regional-level EWAS
DNAm region location: genomic location of the DNA methylation region (chromosome, start position, and end position); Annotated gene(s): gene(s) annotated to the CpGs within the DNA methylation region; N CpGs included: number of CpGs included in the DNA methylation region; Estimate: estimate for the association of maternal sensitivity with DNA methylation at a region; Standard error: standard error for the association of maternal sensitivity with DNAm at a region; Raw p value: unadjusted p value for the association of maternal sensitivity with DNA methylation at a region; Bonferroni adj. p value: p value adjusted for multiple testing with Bonferroni correction.
*This region was not genome-wide significant when siblings were excluded from the sample.
Candidate gene look-up
The candidate gene look-up showed that, of the four selected genes (NR3C1, BDNF, SLC6A4, OXTR), which included 14–74 sites, no CpG met Bonferroni-adjusted gene-wide significance in association with maternal sensitivity (Table 2, online Supplementary Table S4). Only three sites reached nominal significance (p < 0.05).
Table 2. The association between maternal sensitivity and DNA methylation: candidate gene look-up
Gene: candidate gene; Chr: chromosome; N probes: number of probes annotated to the gene (based on the Infinium 450 K); Gene-level sign.: gene-level Bonferroni significance in any of the probes annotated to the candidate gene (p < 0.05/number of annotated probes); Nominal sign.: nominal significance in any of the probes annotated to the candidate gene (p < 0.05); Estimate range: range of estimates for the probes annotated to the candidate genes; % Positive associations: percentage of probes with a positive association with maternal sensitivity; % Negative associations: percentage of probes with a negative association with maternal sensitivity.
Gene ontology
The in-house GO analysis, based on sites with p < 0.001 in the probe-level EWAS, revealed enrichment for 148 pathways. Yet, this threshold might have been overinclusive. Thirty-nine of the 148 pathways were confirmed by the missMethyl GO method (p < 0.05) (online Supplementary Table S5). Both methods indicated enrichment for, among others, calcium ion channels functioning, phosphorylation, and tissue and cell polarity.
Follow-up analyses
Firstly, an mQTL search revealed that five of the 13 significant DNAm regions contained at least one CpG associated with one or more known SNPs (Table 3, online Supplementary Table S6). Eight regions, including ZBTB22/TAPBP (one of our top regions), did not present any mQTLs. Of the 143 sites within the 13 significant regions, 22% (n = 31) associated with one or more known SNPs. All associations were in cis.
Table 3. mQTLs within the statistically significant DNA methylation regions
DNAm region location: genomic location of the DNA methylation region (chromosome, start position, and end position); Annotated gene(s): gene(s) annotated to the DNA methylation region; N CpGs included: number of CpGs included in the DNA methylation region; N mQTL associations: number of SNPs–DNA methylation associations at a region; N CpGs with mQTLs: number of CpGs presenting one or more mQTL(s) at a region; % CpGs with mQTLs: percentage of CpGs presenting one or more mQTL(s) at a region.
Secondly, after additional adjustments for socioeconomic status and maternal psychopathology, associations attenuated at seven regions (median: −1%, range: −44% to 13%). Regions which did not decrease in effect were TAPBP, RNF39, two non-annotated regions, ANKMY1, and ALOX12P2 (online Supplementary Table S7). When adjusting for pre-exposure DNAm levels (online Supplementary Table S8), associations attenuated at 10 regions (median: −45%, range: −97% to 17%), with RNF39 being the most affected. Regions whose estimates did not decrease were ZBTB12, FBXO44/FBXO2, and a non-annotated region (chromosome 7). The median correlation between each CpG DNAm level at birth and age 6 was of ρ = 0.43 (range: 0.11–0.86) (online Supplementary Table S9).
Thirdly, in a follow-up GO analysis, based on the sites within the significant DNAm regions (n = 143), enrichment was found at 63 pathways (in-house method). Of these, 33 were validated by missMethyl (p < 0.05). Both methods indicated enrichment for, among others, several lipoprotein processes (e.g. particle remodeling), and peptide binding (online Supplementary Table S10).
Lastly, of the 13 significant DNAm regions, six contained half or more sites with greater than average blood–brain tissue concordance (Edgar et al., Reference Edgar, Jones, Meaney, Turecki and Kobor2017) in at least one brain tissue (for BA7 r > |0.36|, for BA10 r > |0.40|, for BA20 r > |0.33|), for a total of 67 sites (online Supplementary Table S11) (not empirically tested).
Discussion
This is the first epigenome-wide study investigating the prospective association between typical variation in maternal sensitivity (observed) and offspring DNAm, in a general population of children. Genome-wide significant associations were observed at 13 DNAm regions, four of which did not contain mQTLs and were minimally affected by adjustments for postnatal confounders and by pre-exposure DNAm levels, thus showing robustness in associations.
Summary of key findings
Our first aim was to examine the prospective relationship between maternal sensitivity and child DNAm using complementary approaches. Firstly, no individual CpG was identified in the probe-level EWAS after genome-wide correction. This might indicate that associations at site-level are subtle and challenging to identify, especially considering this study assessed typical variation in maternal care as opposed to extreme deviations (e.g. abuse). The high multiple testing correction burden that probe-level EWASs entail may also impede the detection of single sites of small effect, which could be uncovered with larger samples. For instance, with our sample (N = 235) and model (multiple linear regression, 10 predictors), 80% power, and a genome-wide threshold, only moderate estimates (as small as 0.27) could be detected.
When employing a regional approach, which can detect weaker but more widespread signals by accounting for correlations across CpGs, 13 DNAm regions were significantly associated with maternal sensitivity (p < 1.06 × 10−07, α = 0.05). These findings support the presence of offspring methylomic signatures of maternal care, which may be best uncovered through hypothesis-free approaches with methods capturing the correlational patterns of DNAm. Yet, replication of these findings is needed, and the possibility of false-positive findings should not be excluded. Notably, when considering a more stringent significance threshold (p < 2.18 × 10−09; α = 0.001), as suggested to reduce false-positive rates (Colquhoun, Reference Colquhoun2014), most of the regions (77%, N = 10) remained significantly related to maternal sensitivity.
Further, we failed to detect an association between maternal sensitivity and DNAm variation at candidate genes previously identified by studies of maternal care in humans (Bosmans et al., Reference Bosmans, Young and Hankin2018; Conradt et al., Reference Conradt, Hawes, Guerin, Armstrong, Marsit, Tronick and Lester2016; Provenzi et al., Reference Provenzi, Fumagalli, Giorda, Morandi, Sirgiovanni, Pozzoli and Montirosso2017; Unternaehrer et al., Reference Unternaehrer, Meyer, Burkhardt, Dempster, Staehli, Theill and Meinlschmidt2015). Inconsistencies may reflect several factors, including differences in sample characteristics (e.g. psychiatric v. population-based samples), maternal care assessments (retrospective v. prospective reports), and analysis (e.g. gene regions covered by pyrosequencing v. Infinium 450 K). Lastly, candidate gene studies may be particularly vulnerable to false positives, as shown in the genetic field (Sullivan, Reference Sullivan2007).
As a second aim, we explored whether identified maternal sensitivity–DNAm associations may be influenced by genetic factors, based on mQTL mapping. Twenty-two percent of the sites in our significant regions were linked to known SNPs. This suggests that associations for those sites may be in part confounded by genetic factors and corroborates previous research highlighting DNAm responsiveness to both external exposures and genetic variation (Ladd-Acosta & Fallin, Reference Ladd-Acosta and Fallin2016). However, the presence of mQTLs alone does not preclude environmental effects. Indeed, recent studies have found that interindividual variability in DNAm is primarily explained by gene–environment combinations (additive and interactive effects) (Czamara et al., Reference Czamara, Eraslan, Page, Lahti, Lahti-Pulkkinen, Hämäläinen and Binder2019; Teh et al., Reference Teh, Pan, Chen, Ong, Dogra, Wong and Holbrook2014). Moreover, mQTLs were identified based on a publicly available database, as our sample was underpowered to directly test for genetic confounding. Future studies employing genetically-sensitive designs could more precisely quantify the effect of maternal sensitivity on DNAm by directly modeling genetic influences.
When exploring the robustness of findings to additional adjustments, we observed attenuations at half of the regions, after controlling for socioeconomic status and maternal psychopathology. When considering pre-exposure DNAm levels, estimates attenuated at most regions. Although neonatal methylomic patterns were measured in cord blood at birth and not in peripheral blood (used at age 6), which may lead to additional differences, these findings indicate that associations partly reflected pre-existing DNAm levels. This was clearly exemplified by RNF39, a region strongly associated with sensitivity, robust to postnatal confounders, and genetic influences. After adjustments, its estimate reduced by 97%, showing that associations did not result from postnatal caregiving, as they were already present at baseline (birth). These findings cast doubts on previous studies of caregiving which did not consider pre-exposure DNAm levels, and raise questions on the directionality of associations between maternal care and DNAm, as well as on the potential role of other factors affecting child DNAm (at birth and in childhood) and maternal sensitivity (e.g. shared genetics, maternal distress).
Here, we highlight four ‘high-confidence’ associations with maternal caregiving, which were not linked to any mQTLs, and were most robust to adjustments for confounders and pre-exposure DNAm levels. These spanned (i) ZBTB22/TAPBP, (ii) ZBTB12, (iii) DOCK4, and (iv) a non-annotated region in chromosome 4. All four genes are protein-coding (Geer et al., Reference Geer, Marchler-Bauer, Geer, Han, He, He and Bryant2010). DOCK4 is implicated in neuronal processes, such as neuronal migration, and dendritic arborization (Shi, Reference Shi2013) and its DNAm region presented higher than average blood–BA10 concordance in this study. ZBTB22 and ZBTB12 are involved in transcriptional regulation and nuclear chromatin localization (Agapite et al., Reference Agapite, Albou, Aleksander, Argasinska, Arnaboldi, Attrill and Yook2020). These two genes, together with TAPBP, are within the major histocompatibility complex (MHC). While these associations should be carefully interpreted as the MHC is characterized by extensive linkage disequilibrium (Price et al., Reference Price, Weale, Patterson, Myers, Need, Shianna and Reich2008), this genomic region plays an important role in immune functioning and has been implicated in neuronal plasticity (Shatz, Reference Shatz2009; Sobue et al., Reference Sobue, Ito, Nagai, Shan, Hada, Nakajima and Yamada2018). TAPBP specifically is involved in MHC class I protein complex assembly, gene expression regulation, and immunodeficiency (Agapite et al., Reference Agapite, Albou, Aleksander, Argasinska, Arnaboldi, Attrill and Yook2020). In this study, enrichment for MHC class I protein assembly and peptide binding was found for maternal sensitivity, potentially suggesting that such exposure might enact on TAPBP-related functions via DNAm.
Generally, our high-confidence genes have been previously associated with psychological and developmental problems, inflammation, and stress responses. Molecular changes were shown at TAPBP for major depressive disorder and suicide (Murphy et al., Reference Murphy, Crawford, Dempster, Hannon, Burrage, Turecki and Mill2017), TAPBP and DOCK4 for schizophrenia (Alkelai et al., Reference Alkelai, Lupoli, Greenbaum, Kohn, Kanyas-Sarner, Ben-Asher and Lerer2012; Lee, Kim, & Song, Reference Lee, Kim and & Song2013; Zhang et al., Reference Zhang, You, Li, Long, Zhu, Teng and Zeng2020), ZBTB22 for intellectual disability (Agapite et al., Reference Agapite, Albou, Aleksander, Argasinska, Arnaboldi, Attrill and Yook2020) and psychopathologies following hypercortisolism (Glad et al., Reference Glad, Andersson-Assarsson, Berglund, Bergthorsdottir, Ragnarsson and Johannsson2017), and DOCK4 for autism and dyslexia (Liang et al., Reference Liang, Wang, Zou, Wang, Zhou, Sun and Tomoda2014; Maestrini et al., Reference Maestrini, Pagnamenta, Lamb, Bacchelli, Sykes, Sousa and Monaco2010). Enrichment for pathways including Dock4 has been repeatedly associated with stress-related responses in mice (Lee et al., Reference Lee, Chang, Yeom, Kim, Choi, Shim and Hahm2005; Lisowski et al., Reference Lisowski, Juszczak, Goscik, Wieczorek, Zwierzchowski and Swiergiel2011; Papale, Madrid, Li, & Alisch, Reference Papale, Madrid, Li and Alisch2017), while ZBTB12 DNAm is related to markers of inflammation (e.g. white blood cell counts) (Noro et al., Reference Noro, Gianfagna, Gialluisi, De Curtis, Di Castelnuovo, Napoleone and Izzi2019).
Limitations and suggestions for future research
Our findings should be interpreted in light of several limitations. Firstly, identified associations may have been influenced by additional parental factors that we could not control for in the present study, either because this information was not available (e.g. parental temperament, parental genotype) or due to the low number of cases (e.g. maternal medication and substance use in pregnancy). Nevertheless, we did control for the most important maternal confounders (smoking during pregnancy, socioeconomic status, psychopathology). Secondly, if unmeasured changes in maternal sensitivity and covariates occurred during the 2–3-year time-lag between our exposure and outcome, noise would be introduced in the identified associations. A prospective design, as opposed to a cross-sectional one, remains however preferable due to the possibility to better understand the directionality of associations. Nonetheless, repeated postnatal measurements of both DNAm and maternal sensitivity would be ideal to longitudinally examine how associations change over time and disentangle directionality. Thirdly, we did not have information on whether the mothers included in this study were primary or secondary caregivers (at 4 years only). Yet, within Generation R, most mothers are primary caregivers (White et al., Reference White, Muetzel, El Marroun, Blanken, Jansen, Bolhuis and Tiemeier2018). Additionally, while the use of the Infinium 450 K provided novel insights into the genes affected by maternal sensitivity, future research should employ, when possible, the EPIC 850 K array due to its wider and more diverse genomic coverage (Illumina, 2020). Lastly, our investigation solely focused on the association of maternal sensitivity with the child methylome. Related molecular signatures, such as transcription changes and epigenetic clocks, could be examined in future research to better understand the biological consequences of maternal care.
In conclusion, this exploratory population-based study suggests a prospective association of typical variation in maternal sensitivity with epigenome-wide DNAm in children. We highlight four DNAm regions that showed the strongest associations with maternal sensitivity as well as limited evidence of genetic and pre-exposure influences, and which should thus be prioritized in future confirmatory research. These results permit further delineation of the relationship between DNAm and maternal care in humans and warrant corroboration by future research with large, longitudinal, and genetically-sensitive studies.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291720004353
Acknowledgements
We thank children and parents, midwives, general practitioners, hospitals, and pharmacies in Rotterdam for their contribution to Generation R, the Genetic Laboratory of the Department of Internal Medicine at Erasmus MC for the generation and management of the Infinium 450 K, Verbiest, Higgins, Jhamai, Dr Stolk, and Verkerk for their aid in the EWAS database creation, Dr Teumer for the contribution to quality control and normalization, and Dr Suderman for answering all our queries on the dmrff method.
Financial support
The Generation R study is supported by Erasmus MC, Erasmus University Rotterdam, the Rotterdam Homecare Foundation, the Municipal Health Service Rotterdam area, the Stichting Trombosedienst & Artsenlaboratorium Rijnmond, the Netherlands Organization for Health Research and Development (ZonMw), and the Ministry of Health, Welfare and Sport. DNAm data were funded by the Netherlands Genomics Initiative/Netherlands Organization for Scientific Research (NWO), Netherlands Consortium for Healthy Aging (project 050-060-810), the National Institute of Child and Human Development (R01HD068437), and by the Department of Internal Medicine (Genetic Laboratory) at Erasmus MC. This project received funding from the European Union's Horizon 2020 program (project 733206). The authors are supported by the Dutch Ministry of Education, Culture, and Science and the NWO (project 024.001.003 for AN, HT, MJB-K, and MHvI), the Canadian Institutes of Health Research (AN), an NWO-VICI grant (NWO-ZonMW: 016.VICI.170.200 for HT, LDA), an NWO VENI grant (project 91618147 for JR), the European Joint Programming Initiative ‘A Healthy Diet for a Healthy Life’ (project 529051022; JFF), the European Research Council grant (ERC AdG 669249, MJB-K), the Netherlands Organization for Scientific Research Spinoza Prize (MHvI), and the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant (project 707404 and grant agreement 848158 EarlyCause Project for CAMC).
Conflict of interest
None.
Ethical standards
The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.
