Harry Harlow's seminal experiments in rhesus macaques revealed that affective, soft touch is a motivating force within mammalian infants over and beyond the drive for food. Infant macaques were reared within cages containing one wire mesh surrogate mother offering milk and another surrogate mother covered with a soft terry cloth, not providing sustenance. To the surprise of the researchers, the infants merely fed from the wire mother, returning to the soft mother for the remainder of the day. In other experiments, infants used soft mothers as sources of security when facing novelty and exploring new surroundings. In the absence of soft mothers in these situations, the infants were frantic and distressed. Based on these experiments, Harlow concluded that contact comfort is overwhelmingly important for normative social development (Harlow, Reference Harlow1958).
In recent years, normative variation in tactile stimulation in rats during the first week of life was shown to have lasting consequences for behavioral and neuroendocrine responses to stress in offspring (Caldji et al., Reference Caldji, Tannenbaum, Sharma, Francis, Plotsky and Meaney1998; Liu et al., Reference Liu, Diorio, Tannenbaum, Caldji, Francis, Freedman and Meaney1997). Researchers leveraged differences in the frequency of licking and grooming behaviors (LG) of rat mothers to test whether differential exposure to tactile stimulation was responsible for the later manifestation of behavioral differences, or if the association was due to common genetic factors that distinguish high and low LG rats. By cross-fostering pups from high and low LG mothers, it became clear that maternal care in the first week of life had lasting consequences for stress reactivity phenotypes in pups regardless of genetic background (Liu, Diorio, Day, Francis, & Meaney, Reference Liu, Diorio, Day, Francis and Meaney2000; Meaney, Reference Meaney2010).
This work led to investigations of epigenetic mechanisms, modifications made to DNA and chromatin structure that provide an additional layer of information “upon” the genome. High LG levels in the first week of life increased expression of the glucocorticoid receptor (GR) in the hippocampus, associated with enhanced negative feedback regulation on corticotropin releasing hormone (CRH; a key regulator of the stress response) and behavioral responsiveness to stress (Burton et al., Reference Burton, Chatterjee, Chatterjee-Chakraborty, Lovic, Grella, Steiner and Fleming2007; Jutapakdeegul, Casalotti, Govitrapong, & Kotchabhakdi, Reference Jutapakdeegul, Casalotti, Govitrapong and Kotchabhakdi2003). An epigenetic mechanism was reported to account for this association, specifically, high LG levels were shown to have lower levels of methylation in a transcription factor, nerve growth factor 1 (NGF1), binding site within exon 17 of the GR gene, nuclear receptor subfamily 3, group C, member 1 (NR3C1; exon 1F in humans), whereas low LG levels related to higher methylation at this site (Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004).
DNA methylation (DNAm) describes the presence of the covalent chemical addition of a methyl to a cytosine base, most commonly adjacent to a guanine base (i.e., a CpG dinucleotide). DNAm of gene promoters generally blocks the binding of transcription elements that regulate gene activity, and thus typically higher DNAm at the promoter corresponds with silencing of gene activity, and lower DNAm with gene expression (Illingworth & Bird, Reference Illingworth and Bird2009; Saxonov, Berg, & Brutlag, Reference Saxonov, Berg and Brutlag2006). In the aforementioned rat studies, low LG and high DNAm of the promoter prevented the binding of the NGF1 transcription factor and expression of NR3C1 (Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004). In turn, lower expression of GR receptors corresponded with less negative feedback on the stress response, and increased stress reactivity in offspring.
Unlike DNA sequence, DNAm undergoes dramatic changes synchronized to developmental processes, guiding gene expression. From conception, DNAm plays a critical role in embryogenesis, with coordinated phases of methylation and demethylation guiding the fate of cells in the formation of the high-level organization of the body and brain (Wu & Zhang, Reference Wu and Zhang2014). In contrast to these persistent DNAm marks tightly linked to the expression profiles specific to cell lineage, DNAm sites across the genome are also malleable to environmental exposures (e.g., Gao, Jia, Zhang, Breitling, & Brenner, Reference Gao, Jia, Zhang, Breitling and Brenner2015), and exhibit predictable, progressive changes over the course of aging (Horvath et al., Reference Horvath, Zhang, Langfelder, Kahn, Boks, van Eijk and Ophoff2012). Thus, DNAm offers a molecular basis for the convergence of genetic and environmental information in the biological embedding of experiences over development (Boyce & Kobor, Reference Boyce and Kobor2015; Meaney, Reference Meaney2010), as shown in rats for the embedding of early tactile experiences, and other examples of early exposures (e.g., smoking and maternal mood in utero).
Given these intriguing attributes of DNAm in development and biological plasticity, numerous studies have tested whether early human postnatal exposures similarly relate to DNAm of NR3C1 in humans. Consistent with the LG rat model, a meta-analysis of existing literature concluded that adversities early in life relate to very small but consistent increased DNAm in NR3C1 exon 1F in humans, primarily in peripheral tissues (Turecki & Meaney, Reference Turecki and Meaney2016). For instance, children exposed to physical abuse demonstrated increased methylation of NR3C1 exon 1F in leukocytes (Romens, McDonald, Svaren, & Pollak, Reference Romens, McDonald, Svaren and Pollak2015). Although effects of early adversity in humans are in the same direction as the original finding in rats, this promoter region demonstrates very low DNAm levels, and proportionately very small effects relative to the originally reported effect of LG at this locus.
Other studies investigating DNAm signatures across the genome suggest that early institutional care (Esposito et al., Reference Esposito, Jones, Doom, MacIsaac, Gunnar and Kobor2016; Naumova et al., Reference Naumova, Lee, Koposov, Szyf, Dozier and Grigorenko2012) and child abuse (Suderman et al., Reference Suderman, Borghol, Pappas, Pinto Pereira, Pembrey, Hertzman and Szyf2014) correspond with DNAm at sites at later developmental time points, at least in peripheral tissues. However, the extent that epigenetic marks capture the specific experience of contact between caregiver and infant is unknown. In the following, we discuss the biological pathways of tactile contact in mammals, existing theory on how the experience of touch and broader care affects development and plasticity, and finally make a case for the investigation of the downstream DNAm correlates of neonatal experiences of tactile contact in humans.
Neurobiological Pathways of Touch
For mammals, tactile stimulation is critical to the formation of pair bonds, relying on the action of endogenous β-endorphin at mu-opiate receptors (MORs). β-endorphin is a neuropeptide produced in the brainstem, which serves as a neurotransmitter in neural systems regulating pain, emotional distress, reward, and physiology (Machin & Dunbar, Reference Machin and Dunbar2011) via action upon receptors distributed throughout subcortical and cortical brain areas. β-endorphin has a particularly high affinity for MORs (Bodnar, Reference Bodnar2014; Kieffer & Evans, Reference Kieffer and Evans2009), and action at these receptors has been shown to be the necessary mechanism that induces learning of maternal features and, consequently, bond formation in response to tactile stimulation in rats (Roth & Sullivan, Reference Roth and Sullivan2006) and sheep (Shayit, Nowak, Keller, & Weller, Reference Shayit, Nowak, Keller and Weller2003). Thus, touch elicits the release of endogenous β-endorphin, and binding at MORs is a biological requirement for attachment.
This action of endogenous β-endorphin has rewarding and analgesic effects that foster such a bond (Le Merrer, Becker, Befort, & Kieffer, Reference Le Merrer, Becker, Befort and Kieffer2009; Machin & Dunbar, Reference Machin and Dunbar2011). The rewarding and stress-reducing nature of β-endorphin action is apparent in animal models of separation distress: separation from attachment figures is aversive and induces emotional distress in unison with reduced MOR activation (Machin & Dunbar, Reference Machin and Dunbar2011), which is reversed by MOR agonists (Barr et al., Reference Barr, Schwandt, Lindell, Higley, Maestripieri, Goldman and Heilig2008; Machin & Dunbar, Reference Machin and Dunbar2011). The reward and physiological calm induced by the action of β-endorphin at MORs becomes associated with aspects of the caregiver, conditioning infants to prefer and enjoy contact, which is the basis for bond formation (Dunbar, Reference Dunbar2010; Machin & Dunbar, Reference Machin and Dunbar2011; Sullivan et al., Reference Sullivan, Taborsky-Barba, Mendoza, Itano, Leon, Cotman and Lott1991).
The psychological effects of soft touch and release of endogenous opioids has also been well demonstrated in humans. Soft touch induces rewarding and pleasurable affect as well as reduced physiological arousal and pain (Koepp et al., Reference Koepp, Hammers, Lawrence, Asselin, Grasby and Bench2009; McGlone & Spence, Reference McGlone and Spence2010) via activation of specialized C-tactile fibers in the hairy skin (Löken, Wessberg, Morrison, McGlone, & Olausson, Reference Löken, Wessberg, Morrison, McGlone and Olausson2009; Walker, Trotter, Swaney, Marshall, & McGlone, Reference Walker, Trotter, Swaney, Marshall and McGlone2017). Neuroimaging studies have demonstrated that individual differences in MOR receptor availability in the orbital frontal cortex and insula are related to pain tolerance (Mueller et al., Reference Mueller, Klega, Buchholz, Rolke, Magerl, Schirrmacher and Schreckenberger2010), and that inactivation of β-endorphin activity at MORs induces psychological states of sadness (Zubieta et al., Reference Zubieta, Smith, Bueller, Xu, Kilbourn, Jewett and Stohler2001) and reduced feelings of social connection (Inagaki, Ray, Irwin, Way, & Eisenberger, Reference Inagaki, Ray, Irwin, Way and Eisenberger2016). Based on the rewarding and stress-reducing actions of β-endorphin demonstrated in animals and humans, researchers have theorized that this system is critical to the formation as well as the maintenance of social bonds through the life course (Depue & Morrone-Strupinsky, Reference Depue and Morrone-Strupinsky2005; Dunbar, Reference Dunbar2010).
Oxytocin (OXT) is another neuropeptide critical to social bonding, primarily produced in the hypothalamus, and diffused broadly both centrally and peripherally to regulate emotional states, the hypothalamus–pituitary–adrenal axis, and the autonomic nervous system (Carter, Reference Carter2014; Stoop, Reference Stoop2014). Similar to β-endorphin, OXT is released in response to affiliative interactions (Crockford, Deschner, Ziegler, & Wittig, Reference Crockford, Deschner, Ziegler and Wittig2014). Rather than direct effects on reward, OXT has been shown to have stress-reducing effects, and approach-enhancing effects in particular contexts (Crockford et al., Reference Crockford, Deschner, Ziegler and Wittig2014; Taylor, Reference Taylor2006), via enhanced social cognition (Heinrichs & Domes, Reference Heinrichs and Domes2008). Collectively, this literature suggests that OXT and β-endorphins are critical to the formation of social bonds (Lim & Young, Reference Lim and Young2006), as together, these two peptides regulate reward learning and memory processes, and attenuate the stress response when released in response to bonding activities.
An intriguing possibility is that contact early in life can shape the development of these neurobiological pathways via experience-dependent neuroplasticity (Moore & Depue, Reference Moore and Depue2016). Early postnatal development is a period of heightened plasticity to outer experiences, encompassing a “critical period” or window in which experiences have accentuated and irreversible effects on developing neural circuitries (Fox, Levitt, & Nelson, Reference Fox, Levitt and Nelson2010). Critical periods of neural development may render heightened lability to the effects of the social environment, including the amount and quality of contact with caregivers. The experience of greater frequency of contact may enhance experiences of reward and comfort, strengthening circuitries and pathways establishing social reward and bonding. Moreover, the stress-reducing experience of contact may enhance regulatory capacity over stress, as caregivers become sources of support and comfort in the face of challenges. Thus, the safety, calm, and reward experienced due to high levels of care may perpetuate and engrain these experiences within the developing brain in a manner that facilitates subsequent social functioning and self-regulation. Next, we review the evidence supporting this possibility.
Evidence for Postnatal Programming
Experience-dependent plasticity in response to early care is evident in both animal models and human associations. In the LG rat model, the long-term consequences of high levels of tactile stimulation are far reaching. In addition to stress response pathways reviewed above, LG levels also appear to alter the expression of OXT receptors (Beery, McEwen, MacIsaac, Francis, & Kobor, Reference Beery, McEwen, MacIsaac, Francis and Kobor2016; Francis, Champagne, & Meaney, Reference Francis, Champagne and Meaney2000), brain-derived neurotrophic factor (BDNF; Macrì, Laviola, Leussis, & Andersen, Reference Macrì, Laviola, Leussis and Andersen2010), spatial learning (Liu et al., Reference Liu, Diorio, Day, Francis and Meaney2000), and the structure and function of dendrites in the neocortex and hippocampus (Pinkernelle, Abraham, Seidel, & Braun, Reference Pinkernelle, Abraham, Seidel and Braun2009; Smit-Rigter, Champagne, & van Hooft, Reference Smit-Rigter, Champagne and van Hooft2009; Takatsuru, Yoshitomo, Nemoto, Eto, & Nabekura, Reference Takatsuru, Yoshitomo, Nemoto, Eto and Nabekura2009). Similarly, a single exposure to OXT in early life has lasting effects on social behavior in voles (Carter, Boone, Pournajafi-Nazarloo, & Bales, Reference Carter, Boone, Pournajafi-Nazarloo and Bales2009), and variation in tactile stimulation relates to different cortical organization and connectivity in vole offspring (Seelke, Perkeybile, Grunewald, Bales, & Krubitzer, Reference Seelke, Perkeybile, Grunewald, Bales and Krubitzer2016).
In humans, evidence stemming from institutionally reared infants suggests that the absence of basic care is devastating to normative social and cognitive development. In institutional care, one caregiver is responsible for multiple infants, leading to the majority of a caregiver's time with one infant devoted to feeding and cleaning, leaving infants alone in their beds an average of 17.5 hr per day (van IJzendoorn et al., Reference van IJzendoorn, Palacios, Sonuga-Barke, Gunnar, Vorria, McCall and Juffer2011). Institutionally reared infants show delays in physical growth, atypical stress responses, and dysfunction in cognition and social behavior (van IJzendoorn et al., Reference van IJzendoorn, Palacios, Sonuga-Barke, Gunnar, Vorria, McCall and Juffer2011). Institutional rearing continues to link to physical and psychological outcomes across the life span (Humphreys et al., Reference Humphreys, Gleason, Drury, Miron, Nelson, Fox and Zeanah2015; Sigal, Perry, Rossignol, & Ouimet, Reference Sigal, Perry, Rossignol and Ouimet2003), and these effects appear to be reliant on critical periods of development. Institutionally reared children randomly assigned out of orphanages to foster homes prior to the age of 2 have substantial gains in cognitive and social outcomes relative to children who remained within orphanages until later ages (Zeanah, Gunnar, McCall, Kreppner, & Fox, Reference Zeanah, Gunnar, McCall, Kreppner and Fox2011).
Human studies on more normative variation in care typically rely on self-report and observational measures of caregiving sensitivity or quality. These measures can be very broad, typically including multiple components (e.g., sensitivity to distress and nondistress, positive regard/warmth, intrusiveness). In an effort to capture a continuous dimension of care in humans analogous to the LG studies done in rats, one study examined caregiver behavior during routine activities, including changing, dressing, and feeding at 9 months, because these activities are mutually challenging for both the infant and caregiver (Hane & Fox, Reference Hane and Fox2006; Hane & Philbrook, Reference Hane and Philbrook2012). Low quality of caregiving behavior related to concurrent infant biobehavioral profiles of fearfulness to novel cues, impaired attention in a social task, and greater stress reactivity, as well as subsequent social difficulties at 2–3 years of age (Hane, Henderson, Reeb-Sutherland, & Fox, Reference Hane, Henderson, Reeb-Sutherland and Fox2010). Moreover, these associations held after controlling for infant temperament, suggesting that caregiving quality may have a direct impact on the child stress response independent of temperamental factors.
Despite the high impact of the LG rat studies in child development research, limited human work has directly addressed the role of variation in the specific stimulus of tactile contact in shaping child biology. In one study, maternal depression related to decreased vagal withdrawal (corresponding with capacity for emotion regulation) and increased infant negative emotionality, but only in those infants who received low levels of maternal touch (Sharp et al., Reference Sharp, Pickles, Meaney, Marshall, Tibu and Hill2012). Similarly, postnatal soft touch was protective against the negative effects of prenatal mood disorder symptoms on child internalizing and externalizing (Pickles, Sharp, Hellier, & Hill, Reference Pickles, Sharp, Hellier and Hill2017; Sharp et al., Reference Sharp, Pickles, Meaney, Marshall, Tibu and Hill2012). Using the still face procedure (a stressful situation in which mothers cease to exhibit expressivity to their infants) and a vagal marker of experienced stress, another study indicated that infants who were not touched by mothers experienced greater stress than infants who received touch. Moreover, unlike the untouched infants, those receiving touch demonstrated full recovery to baseline stress levels after the procedure (Feldman, Singer, & Zagoory, Reference Feldman, Singer and Zagoory2010; Pratt, Singer, Kanat-Maymon, & Feldman, Reference Pratt, Singer, Kanat-Maymon and Feldman2015). Finally, frequency of maternal touch during a 10-min play session in 4- to 6-year-olds related to greater orientation to social cues (Reece, Ebstein, Cheng, Ng, & Schirmer, Reference Reece, Ebstein, Cheng, Ng and Schirmer2016). Thus, preliminary evidence in humans suggests that tactile contact could similarly have lasting effects on child biology and behavior.
The Current Study: DNAm as a Marker of the Biological Encoding of Tactile Contact
The present study assesses whether DNAm reflects the biological embedding of early postnatal contact in humans. DNAm can be altered by experience in a manner that is relevant to the expression of genes and, consequently, the emergence of biological phenotypes. DNAm can also serve as an indicator of epigenetic age given predictable changes in DNAm across the life course. Here, we explore these two functions of DNAm in light of possible links with early postnatal contact.
DNAm as a mechanism of plasticity
The brain develops in an activity-dependent manner, meaning that unfolding neural circuits rely on signals from the external world to inform what pathways to strengthen, and which connections to prune if not worthy of energy expenditure. These activity-dependent processes that guide brain development are in part orchestrated by dynamic alterations in DNAm (Baker-Andresen, Ratnu, & Bredy, Reference Baker-Andresen, Ratnu and Bredy2013). For instance, the regulation of neurogenesis involves transcription factor binding, depending on the presence and absence of DNAm marks (Wang, Tang, He, & Jin, Reference Wang, Tang, He and Jin2016). Moreover, neuronal activation in the mammalian brain corresponds with dynamic changes in DNAm within neuroregulatory and plasticity-related genes (Guo et al., Reference Guo, Ma, Mo, Ball, Jang, Bonaguidi and Song2011). Similarly, DNAm in the amygdala, hippocampus, and cortex is associated with memory formation and consolidation (Day & Sweatt, Reference Day and Sweatt2011; Yu, Baek, & Kaang, Reference Yu, Baek and Kaang2011). Thus, DNAm is involved in the activity-dependent machinery required for the formation and shaping of neural connections.
Given this function of DNAm, it is possible that variation in DNAm of particular genes related to the encoding of tactile stimulation may reflect differences in early life tactile experiences. Because the activity of β-endorphin at MORs and OXT in response to touch are critical to the formation of social bonds, it is possible that DNAm at genes encoding MOR, μ-opioid receptor M1 (OPRM1), and the OXT receptor (OXTR) is relevant to early contact experiences. Likewise, a third gene, BDNF, is critically involved in neuronal synaptogenesis, demonstrating drastic increases during the early postnatal period (Karpova, Reference Karpova2014). Furthermore, dynamic changes in DNAm at BDNF was previously linked to fear conditioning (Lubin, Roth, & Sweatt, Reference Lubin, Roth and Sweatt2008) and extinction (Baker-Andresen, Flavell, Li, & Bredy, Reference Baker-Andresen, Flavell, Li and Bredy2013). DNAm at BDNF could similarly reflect associative conditioning to tactile contact that is critical to bond formation in early life.
In combination with experiences that may affect DNAm at particular genes, genetic variation is also a key influence of DNAm (Henikoff & Greally, Reference Henikoff and Greally2016). Genetic variants (e.g., single nucleotide polymorphisms [SNPs]) affect nearby DNAm at CpG sites, referred to as methylation quantitative trait loci (mQTLs; Shoemaker, Deng, Wang, & Zhang, Reference Shoemaker, Deng, Wang and Zhang2010). In an assessment of the relative contribution of SNPs versus prenatal exposures to DNAm at sites across the genome, it was found that the combined, interactive effects of genetic factors and prenatal exposures were the most relevant to DNAm patterns (Teh et al., Reference Teh, Pan, Chen, Ong, Dogra, Wong and Holbrook2014). Moreover, the val/met SNP in the human BDNF gene (rs6265; Val66Met) appears to have relevance for the relationship between early exposures and DNAm. Specifically, the strength of the association between antenatal maternal anxiety and DNAm patterns relevant to amygdalar and hippocampal volumes was found to depend on this BDNF variant (Chen et al., Reference Chen, Pan, Tuan, Teh, MacIsaac, Mah and Holbrook2015). Although not explored in relation to DNAm, genetic polymorphisms of OPRM1 and OXTR are similarly reported to relate to differences in responsiveness to the social environment (Gong et al., Reference Gong, Fan, Liu, Yang, Zhang and Zhou2017; Troisi et al., Reference Troisi, Frazzetto, Carola, Di Lorenzo, Coviello, D'Amato and Gross2011, Reference Troisi, Frazzetto, Carola, Di Lorenzo, Coviello, Siracusano and Gross2012; Walum et al., Reference Walum, Lichtenstein, Neiderhiser, Reiss, Ganiban, Spotts and Westberg2012). Thus, in the present study we assess the dual role of common genetic variants of OPRM1, OXTR, and BDNF alongside postnatal contact in molding DNAm within these genes.
DNAm and the epigenetic clock
Across the life span, environmental and stochastic changes accumulate in DNAm throughout the genome and across tissues, and are believed to contribute to increasing diversity over time (Teschendorff, Reference Teschendorff2013). Well-powered studies assessing the DNAm landscape across the genome have identified normative patterns of cell changes associated with age (Teschendorff, West, & Beck, Reference Teschendorff, West and Beck2013). Based on the robustness and consistency of age-related DNAm marks, epigenetic “clocks” have been developed, which encompass sets of specific CpGs that are highly accurate predictors of chronological age (Hannum et al., Reference Hannum, Guinney, Zhao, Zhang, Hughes, Sadda and Zhang2013; Horvath et al., Reference Horvath, Zhang, Langfelder, Kahn, Boks, van Eijk and Ophoff2012; Jones, Goodman, & Kobor, Reference Jones, Goodman and Kobor2015). Epigenetic clocks serve as intriguing, emerging biomarkers of health and disease, with the deviation between epigenetic age and chronological age (i.e., age acceleration) associating with morbidity, mortality, and experiences of psychological stress in adult and aging cohorts (Chen et al., Reference Chen, Marioni, Colicino, Peters, Ward-Caviness, Tsai and Horvath2016; Marioni et al., Reference Marioni, Shah, McRae, Chen, Colicino, Harris and Deary2015; Zannas, Reference Zannas2016).
Because epigenetic clocks were calibrated to primarily adult and aging samples, the interpretation of age acceleration is less straightforward in pediatric samples, where proportionally smaller deviations might be developmentally relevant. In a recent longitudinal study, DNAm age was estimated at birth, age 7, and age 17 and correlated with developmental markers, with conflicting results (Simpkin et al., Reference Simpkin, Howe, Tilling, Gaunt, Lyttleton, McArdle and Relton2017). Specifically, age acceleration positively associated with changes in weight, body mass index, and Tanner stage, but negatively associated with changes in height and fat mass. Given the error of existing clocks (e.g., Horvath clock error is 3.6 years), and that differences in weeks are significant to developmental outcomes for gestational age, specialized clocks that estimate gestational age in cord blood were developed (Bohlin et al., Reference Bohlin, Håberg, Magnus, Reese, Gjessing, Magnus and Nystad2016; Knight et al., Reference Knight, Craig, Theda, Bækvad-Hansen, Bybjerg-Grauholm, Hansen and Smith2016). Gestational age was subsequently found to associate with both maternal risk factors and pregnancy and birth outcomes (Girchenko et al., Reference Girchenko, Lahti, Czamara, Knight, Jones, Suarez and Räikkönen2017). The directions of the associations were mixed: prenatal risk factors predicted both gestational age acceleration (e.g., maternal age at delivery and maternal preeclampsia) and deceleration (e.g., gestational diabetes). Infant lower birth size was also related to age acceleration. The authors concluded that deviations in epigenetic age, in either direction, may mark exposure to risk during pregnancy.
In summary, the biological relevance and implications of deviations in epigenetic age in pediatric samples remains understudied. We assess the relationship between early postnatal care and this DNAm gauge of epigenetic age deviation from chronological age in the present study.
Rationale for current study
Like rats, macaques, and other mammalian animal models, humans rely on tactile contact with caregivers for healthy social and cognitive development. The current study tests whether normative variation in contact in early postnatal life in humans associates with a DNAm biological signature 4–5 years later, which has yet to be addressed given the focus in the human literature on early adversity. We consider multiple levels of child biology that may be relevant to DNAm, including genotype and infant distress levels.
Applying a daily diary strategy, caregivers reported in real time on the frequency and duration of contact with infants, as well as infant distress behaviors, 24-hr per day during the fifth week of life, followed by buccal swab collection of buccal epithelial cells (BECs) when children were 4–5 years of age for DNAm quantification. Because infant crying behavior demonstrates an “n” shaped curve in the first 3 months of life, peaking between 1 and 2 months and demonstrating substantial interindividual variability (Barr, Reference Barr2010; Barr & Gunnar, Reference Barr, Gunnar, Barr, Hopkins and Green2000), infant distress was controlled for and investigated as a moderator of the effects of maternal contact.
Although this investigation is exploratory in nature, we specifically target the methylation of four candidate genes relevant to the neurobiological stress and reward pathways of touch and to neuroplasticity: NR3C1, OPRM1, OXTR, and BDNF. For NR3C1, we specifically tested CpG sites previously shown to associate with LG in rats and early adversity in humans. Given the relevance of genetic variability to DNAm patterns, and the prominence of gene–environment interaction in explaining DNAm patterns (Teh et al., Reference Teh, Pan, Chen, Ong, Dogra, Wong and Holbrook2014), we assessed common SNP variants within OPRM1, OXTR, and BDNF to identify mQTLs and potential interactions with contact and child distress. We also apply a whole epigenome approach to DNAm beyond these four genes to identify novel differentially methylated regions (DMRs) between low and high contact. Finally, we assess the impact of contact on epigenetic age deviation.
Method
Participants
We recruited 1,279 mother–infant dyads from maternity wards in the greater Vancouver area in British Columbia, Canada, through the Parents Helping Infants study (PHI; Barr et al., Reference Barr, Barr, Fujiwara, Conway, Catherine and Brant2009). To qualify, infants had to be born at greater than 34 weeks gestation, not have serious medical conditions, and have parents who could speak and read English. All PHI participants who completed at least 3 complete days of the Baby's Day Diary were included (n = 1,055 dyads) and analyzed for physical contact time. Following the criteria described below, 155 dyads were identified as high contact and 152 dyads as low contact, and were followed up with a written letter to invite participation in a follow-up study involving the collection of biological samples 4–5 years following participation in the PHI study. A phone call followed the initial invitation package 2 weeks later to determine participant interest, answer questions, review eligibility criteria, and schedule appointments to obtain informed consent and sample collection. Children were eligible if they were not currently suffering from any major health problems. BEC and saliva samples were obtained from a total of 94 children either at home or at the BC Women & Children's Health Centre. A demographic questionnaire was administered to both parents at the same time.
The final sample included 39 females and 55 males, 31% identified as North or South Asian, 60% as White, 2% as Black, and 7% as other. The average age of children at follow up was 4.52 years (SD = 0.52).
Measures
Baby's Day Diary
For 4 consecutive days when infants were 5 weeks old, caregivers completed the Baby's Day Diary, an e-diary that is a widely used, reliable, and valid tool for recordings of infant behaviors and interactions in the home (Barr, Kramer, Boisjoly, McVey-White, & Pless, Reference Barr, Kramer, Boisjoly, McVey-White and Pless1988; St. James-Roberts & Plewis, Reference St. James-Roberts and Plewis1996). The ease of use of the Baby's Day Diary is ideal for obtaining accurate gauges of the frequency and duration of behaviors of interest (Lam et al., Reference Lam, Barr, Catherine, Tsui, Hahnhaussen, Pauwels and Brant2010). Although the diaries were primarily completed by mothers, any caregiver who was caring for the baby for this period of time recorded care behaviors. Continuously over 24 hr for 4 days, caregivers recorded six infant behaviors (awake and content, fussing, crying, unsoothable crying, feeding, and sleeping), and any caregiving carrying or holding that involved body contact.
Contact
A continuous measure of physical contact time between infants and caregivers was calculated by adding all caregiver diary-reported time spent holding or carrying infants in hours, and then taking the averaging over the number of days in which complete diary entries were obtained from the caregiver. The average contact time per day across the sample (n = 1,055 dyads) was 9 hr, 7 min (SD = 2 hr, 52 min), and contact ranged broadly from 3 to 23 hr per day. The distribution was unimodal and normal with a slight positive skew (see Figure 1), compatible with that obtained for maternal LG in rodent studies. The approach applied to the LG rat model was used to define high and low contact groups (Champagne, Francis, Mar, & Meaney, Reference Champagne, Francis, Mar and Meaney2003). Specifically, infants were classified as receiving high or low care if their continuous score was one standard deviation above or below, respectively, the mean for daily physical contact time. These groups differed by 6 hr of mean daily physical contact time. In the final recruited sample, high (n = 55) and low (n = 39) groups differed by 8.93 hr.
Figure 1. (Color online) Distribution of tactile contact time in minutes per day. Average contact time in minutes per day for 1,055 mother–infant dyads who completed at least 3 complete days of the Baby's Day Diary. Of these dyads, 155 were identified as high contact and 152 as low contact, represented by the dark gray bars at the tails of the distribution.
Infant distress
In the Baby's Day Diary, caregivers recorded the frequency and duration of infant behaviors including fussing, crying, and unsoothable crying. To account for this important individual difference variable, a principal component analysis (PCA) was performed using the prcomp package in R to condense these six correlated variables to a set of independent composite variables that account for maximal variability in the data. All six variable distributions were normalized and fully standardized prior to the PCA calculation as nonnormal distributions and unscaled data can bias PCA results. The first PC represented crying and unsoothable crying duration and frequency, accounting for 45% of the variance, and the second PC demonstrated strong loadings for fussing duration and frequency, accounting for 25% of the variance. As a comprehensive measure of infant distress, the first PC was used in subsequent analyses to account for maximal variability without running multiple tests.
Genome-wide DNAm assay
BEC collection and DNA extraction
BEC samples were collected using Isohelix Buccal Swabs (Cell Projects, Harrietsham, Kent, England) and stabilized per manufacturer's instructions for storage at room temperature. Genomic DNA was isolated using Isohelix Buccal DNA Isolation Kits (Cell Projects, Harrietsham, Kent, England) and purified and concentrated using DNA Clean & Concentrator (Zymo Research, Irvine, CA).
Illumina infinium methylation assay
A sample of 1000 ng of genomic DNA was used for bisulfite conversion using the EZ DNA Methylation Kit (Zymo Research). Bisulfite deaminates unmethylated cytosine to uracil, while methylated cytosine is protected from deamination, thus converting epigenetic information to sequence-based information. Quantitative DNAm measurements of purified bisulfite converted DNA were performed with the Infinium HumanMethylation450 BeadChip Assay per manufacturer's instructions (Illumina, San Diego, CA). The 450k platform quantifies the DNAm status at over 450,000 CpG sites, covering 99% of available reference sequence genes, representing approximately 1%–2% of total CpG sites in the genome (Bibikova et al., Reference Bibikova, Barnes, Tsan, Ho, Klotzle, Le and Shen2011). A sample of 160 ng of bisulfite converted DNA was whole-genome amplified and fragmented by an enzymatic process and hybridized to BeadChip arrays. Illumina GenomeStudio software (San Diego, CA, ISA) was used to calculate β values by dividing the methylation probe signal intensity by the sum of methylated and unmethylated probe signal intensities. The β values therefore range from 0 (site completely unmethylated across cells) to 1 (fully methylated). Technical replicates across runs achieved r > .99, indicating high reproducibility. The β values were extracted and imported into R statistical software for the remainder of processing and analysis. For all statistical analyses, logit-transformed β values, called M values, are used as these are less heteroscedastic, although β values (which are biologically meaningful) are used for visualization purposes (Du et al., Reference Du, Zhang, Huang, Jafari, Kibbe, Hou and Lin2010).
Preprocessing of DNAm data
A subset of probes were initially removed, included 65 SNP probes, probes located on the X or Y sex chromosomes, and probes with missing values or detection ps > .01 in five or more samples, polymorphic probes, and those shown to cross hybridize (Price et al., Reference Price, Cotton, Lam, Farré, Emberly, Brown and Kobor2013). After filtering, 441,576 probes remained to be carried forward to normalization using the beta mixture quantile method (Teschendorff, Marabita, et al., Reference Teschendorff, Marabita, Lechner, Bartlett, Tegner, Gomez-Cabrero and Beck2013). Because the Illumina platform includes two probe types (Type I and II) with differing distributions, this normalization step is necessary for the two types to have similar, comparable distributions of β values. Next, known technical variation having to do with the specific chip and position on the chip was accounted for using “ComBat” from the R “sva” package (Johnson, Li, & Rabinovic, Reference Johnson, Li and Rabinovic2007). Finally, computational deconvolution was used to extract approximate proportions of cell types in each sample (method derived from Smith et al., Reference Smith, Kilaru, Klengel, Mercer, Bradley, Conneely and Binder2015). It was found that in addition to BECs, samples also contained minority portions of CD34+ leukocyte markers. DNAm at each CpG was regressed onto predicted proportion of BECs for each individual, and the residuals were used to account for this biological variation. To confirm that these technical and biological confounds were removed from the β value matrix, at each step a PCA of DNAm values was performed and correlated with confounds and variables of interest to ensure that variability in DNAm was no longer strongly attributable to batch effects and cell type (Figure 2). Technical batch effects were no longer significantly associated with any of the top DNAm PCs; however, there were still residual effects of CD34+ monocytes. Because these residual effects are independent of contact, we carried forward with analyses. Sex, ethnicity, and age were highly correlated with DNAm variability (Figure 2), and were accounted for in all statistical models.
Figure 2. (Color online) Associations between top DNAm principal components and study variables before and after preprocessing. To ensure preprocessing steps were adequately accounting for unwanted variation in DNAm data, we observed associations among the principal components with technical and biological variables and variables of interest in the study. (a) Prior to correction via preprocessing, there are clear effects of technical (chip and chip position) and biological variables (sex, ethnicity, age, and BEC cell type proportions) on DNAm. (b) After preprocessing, technical batch effects and cell type proportions have been corrected. The remaining variables are included as covariates in statistical models. Chip is the specific array the samples were run on, and chip position is the specific position samples were placed on the chip.
NR3C1 DNAm and candidate genotype pyrosequencing
As we were interested in replicating the original LG findings, we investigated DNAm at NR3C1 exon 1F in more depth using pyrosequencing because the 450k array only represents three CpGs in this region. DNAm at the remaining genes was explored at sites quantified by the 450k array, with the option of pyrosequencing any regions showing significant associations with contact in greater depth. To enhance our ability to detect a contact DNAm signal at OPRM1, OXTR, and BDNF, we additionally incorporated genotype, which may moderate the effects of contact.
Human NR3C1 promoter pyrosequencing primers were designed by EpigenDx Inc. (Catalog No. ADS749), and were used to quantify seven sites within exon 1F, two of which overlap with sites quantified on the 450k array.
PyroMark Assay Design 2.0 (Qiagen Inc.) software was used to design the pyrosequencing assay covering the target candidate gene regions. HotstarTaq DNA polymerase kit (Qiagen Inc.) was used to amplify the target regions using the biotinylated primer set with the following PCR conditions: 15 min at 95 °C, 45 cycles of 95 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s, and a 5-min 72 °C extension step. Streptavidin-coated sepharose beads were bound to the biotinylated-strand of the PCR product and then washed and denatured to yield single-stranded DNA. Sequencing primers were introduced to allow for pyrosequencing (Pyromark™ Q96 MD pyrosequencer, Qiagen Inc.). All primer sequences used to genotype OPRM1, OXTR, and BDNF, and for NR3C1 promoter DNAm are available from the corresponding author upon request.
OPRM1
The A118G SNP (rs1799971) is located in exon 1 of OPRM1 and causes the normal amino acid at residue 40, asparagine to be replaced by aspartic acid. In vitro studies on the effect of A118G or G allele on the affinity of the MOR receptor have been mixed, but it is generally shown to be a gain-of-function allelic variant (Barr et al., Reference Barr, Schwandt, Lindell, Higley, Maestripieri, Goldman and Heilig2008; Bond et al., Reference Bond, LaForge, Tian, Melia, Zhang, Borg and Yu1998; Deb, Chakraborty, Gangopadhyay, Choudhury, & Das, Reference Deb, Chakraborty, Gangopadhyay, Choudhury and Das2010). The G allele is associated with enhanced social attachment in animal models (Barr et al., Reference Barr, Schwandt, Lindell, Higley, Maestripieri, Goldman and Heilig2008; Curley, Reference Curley2011; Higham et al., Reference Higham, Barr, Hoffman, Mandalaywala, Parker and Maestripieri2011; Moles, Kieffer, & D'Amato, Reference Moles, Kieffer and D'Amato2004) and possibly enhanced reactivity to the social environment in humans (Bertoletti, Zanoni, Giorda, & Battaglia, Reference Bertoletti, Zanoni, Giorda and Battaglia2012; Troisi et al., Reference Troisi, Frazzetto, Carola, Di Lorenzo, Coviello, D'Amato and Gross2011; Way, Taylor, & Eisenberger, Reference Way, Taylor and Eisenberger2009). Thus, the G allele may correspond with greater biological responses to social experience, and thus a stronger relationship between postnatal contact and OPRM1 DNAm.
OXTR
Two SNPs within OXTR, rs7632287 (intergenic A-to-G substitution) and rs53576 (intronic A-to-G substitution), have been linked in well-powered meta-analyses to social behavior (Bakermans-Kranenburg & van IJzendoorn, Reference Bakermans-Kranenburg and van IJzendoorn2014; Gong et al., Reference Gong, Fan, Liu, Yang, Zhang and Zhou2017; LoParo & Waldman, Reference LoParo and Waldman2015), although the functional effects of these variants are not known and may be in linkage disequilibrium with functional SNPs. Given that the G allele of rs7632287 and the G allele of rs53576 are associated with higher social responsiveness (Gong et al., Reference Gong, Fan, Liu, Yang, Zhang and Zhou2017; Walum et al., Reference Walum, Lichtenstein, Neiderhiser, Reiss, Ganiban, Spotts and Westberg2012), we expected these alleles to correspond with greater effects of postnatal contact on OXTR DNAm.
BDNF
BDNF protein may enhance plasticity in the brain by supporting neuron growth, differentiation, and survival, especially in hippocampal regions (Huang & Reichardt, Reference Huang and Reichardt2001). The val/met SNP in the human BDNF gene (rs6265; Val66Met) involves a coding change from valine to methionine and has been linked to plasticity processes, with the met allele corresponding with reduced activity-dependent secretion of BDNF in hippocampal neurons (Egan, Weinberger, & Lu, Reference Egan, Weinberger and Lu2003). The general pattern reported in the literature is that the met allele corresponds with heightened effects of environmental exposures on DNAm and on neural phenotypes (Chen et al., Reference Chen, Pan, Tuan, Teh, MacIsaac, Mah and Holbrook2015; Gunnar et al., Reference Gunnar, Wenner, Thomas, Glatt, McKenna and Clark2012; Hayden et al., Reference Hayden, Klein, Dougherty, Olino, Dyson, Durbin and Singh2010; Juhász, Földi, & Penke, Reference Juhász, Földi and Penke2011; Mata, Thompson, & Gotlib, Reference Mata, Thompson and Gotlib2010; Willoughby, Mills-Koonce, Propper, & Waschbusch, Reference Willoughby, Mills-Koonce, Propper and Waschbusch2013). Thus, this polymorphism may moderate any potential associations between postnatal contact and DNAm at BDNF.
Analytic strategy
To identify if DNAm relates to early postnatal care, we applied three different analyses.
Candidate analysis
The first analysis targeted CpGs from candidate genes of interest, which were separated out from the remaining sites on the 450k array. For these candidate genes, we were interested in replicating a former association of contact with NR3C1, and in testing associations between contact and CpGs on genes related to the encoding of tactile experiences (OPRM1 and OXTR) and neuroplasticity (BDNF). For the latter three genes, we incorporated common genetic variants as potential moderators of contact and child distress.
We selected CpGs from the 450k array identified as located proximal to or within candidate genes of interest according to the Illumina annotation (“UCSC_RefGene_name”), which ranged between 17 and 90 CpGs (Table 1). This isolation of CpGs from candidate regions limits the number of statistical tests, enhancing power for these target sites. To further reduce the number of tests, we filtered down candidate CpGs to only sites that were variable among individuals in the sample. Invariable probes were identified and removed based on interquartile range (IQR). Specifically, the IQR was calculated between the 90th and 10th sample percentiles at each CpG, to account for outliers (Lemire et al., Reference Lemire, Zaidi, Ban, Ge, Aïssi, Germain and Hudson2015). Probes below an IQR threshold of 5% were removed, meaning that the range of DNAm at that CpG was below 5% across samples, leaving 7–19 CpGs per candidate gene for a total of 47 sites across the four genes.
Table 1. Number of probes and CpG islands before and after variability filtering
It should be noted that many invariable sites are located within CpG islands that tend to correspond with gene promoters. Gene promoters are commonly unmethylated and thus more likely to be invariable. By limiting to variable probes, a lower proportion of sites tested are identified as located within CpG islands (Table 1). Methylation status of each of these CpGs serves as the dependent variable in the final models. For all candidate sites (including those investigated below as mQTLs), we ran both a main effects model of contact and an interaction model allowing the impact of contact to depend on child distress:
Main effects model:
$$\eqalign{ {\rm DNAm} =\;& {\rm \beta} _{\rm 0} + {\rm \beta} _{\rm 1} \times {\rm Contact} + {\rm \beta} _{\rm 2} \times {\rm Distress} + {\rm \beta} _{\rm 3} \cr & \times {\rm Sex} + {\rm \beta} _{\rm 4} \times {\rm Ethnicity} + {\rm \beta} _{\rm 5} \times {\rm Age}{\rm.}} $$
Interaction model:
$$\eqalign{ {\rm DNAm} =\;& {\rm \beta} _{\rm 0} + {\rm \beta} _{\rm 1} \times {\rm Contact} + {\rm \beta} _{\rm 2} \times {\rm Distress} + {\rm \beta} _{\rm 3} \cr &\! \times {\rm Sex} + {\rm \beta} _{\rm 4} \times {\rm Ethnicity} + {\rm \beta} _{\rm 5} \times {\rm Age} + {\rm \beta} _{\rm 6} \cr & \! \times {\rm Contact} \times {\rm Distress} + {\rm \beta} _{\rm 7} \times {\rm Contact} \times {\rm Sex \,+}\, {\rm \beta} _{\rm 8} \cr & \! \times {\rm Contact} \times {\rm Ethnicity} + {\rm \beta} _{\rm 9} \times {\rm Contact} \times {\rm Age} \cr & \! + {\rm \beta} _{{\rm 10}} \times {\rm Distress} \times {\rm Sex} + {\rm \beta} _{{\rm 11}} \times {\rm Distress} \cr & \!\times {\rm Ethnicity} + {\rm \beta} _{{\rm 12}} \times {\rm Distress} \times {\rm Age}{\rm.}} $$
These same models were performed on each of the sites that were pyrosequenced for NR3C1 separately. Note, interaction terms with covariates are included in all interaction models for proper adjustment (Hull, Tedlie, & Lehn, Reference Hull, Tedlie and Lehn1992; Yzerbyt, Muller, & Judd, Reference Yzerbyt, Muller and Judd2004).
To identify mQTLs, methylation at each probe was regressed onto genotype for each candidate gene, and those with significant associations following multiple test correction (Benjamini–Hochberg method) were considered mQTLs. For identified mQTL sites, the following two statistical models were applied for each site, to account for possible interactions between genotype and both contact and child distress:
$$\eqalign{ {\rm DNAm} =\;& {\rm \beta} _{\rm 0} + {\rm \beta} _{\rm 1} \times {\rm Contact} + {\rm \beta} _{\rm 2} \times {\rm Genotype} + {\rm \beta} _{\rm 3} \cr & \times {\rm Sex} + {\rm \beta} _{\rm 4} \times {\rm Ethnicity} + {\rm \beta} _{\rm 5} \times {\rm Age} + {\rm \beta} _{\rm 6} \cr &\times {\rm Contact} \times {\rm Genotype} + {\rm \beta} _{\rm 7} \times {\rm Contact} \times {\rm Sex} \cr & + {\rm \beta} _{\rm 8} \times {\rm Contact} \times {\rm Ethnicity} + {\rm \beta} _{\rm 9} \times {\rm Contact} \cr & \times {\rm Age} + {\rm \beta} _{{\rm 10}} \times {\rm Genotype} \times {\rm Sex} + {\rm \beta} _{{\rm 11}} \cr & \times {\rm Genotype} \times {\rm Ethnicity} + {\rm \beta} _{{\rm 12}} \times {\rm Genotype \times Age}{\rm.} \cr {\rm DNAm} =\;& {\rm \beta} _{\rm 0} + {\rm \beta} _{\rm 1} \times {\rm Distress} + {\rm \beta} _{\rm 2} \times {\rm Genotype} + {\rm \beta} _{\rm 3} \cr & \times {\rm Sex} + {\rm \beta} _{\rm 4} \times {\rm Ethnicity} + {\rm \beta} _{\rm 5} \times {\rm Age} + {\rm \beta} _{\rm 6} \cr & \times {\rm Distress} \times {\rm Genotype} + {\rm \beta} _{\rm 7} \times {\rm Distress} \times {\rm Sex} \cr & + {\rm \beta} _{\rm 8} \times {\rm Distress} \times {\rm Ethnicity} + {\rm \beta} _{\rm 9} \times {\rm Distress} \cr & \times {\rm Age} + {\rm \beta} _{{\rm 10}} \times {\rm Genotype} \times {\rm Sex} + {\rm \beta} _{{\rm 11}} \cr & \times {\rm Genotype} \times {\rm Ethnicity} + {\rm \beta} _{{\rm 12}} \times {\rm Genotype} \times {\rm Age}{\rm.} } $$
We did not include three-way interactions (Genotype × Contact × Distress) given the limited sample size, and correspondingly sparse data within the groups created from such higher order interactions.
Epigenome-wide association analysis
The second analysis was an epigenome-wide association analysis, in which we tested whether probes across the genome related to high versus low contact (although the 450k array only quantifies ~1.6% of CpGs in the genome). This analysis assessed whether sites across the genome associated with contact. Given the limited sample size and large number of tests for this exploratory analysis, we focused on simple, main effects of contact, and applied a fairly relaxed multiple test correction threshold of 0.1.
After removing candidate sites, 441,398 CpGs remained for epigenome-wide analysis. We used a package in R, “DMRcate” (Peters et al., Reference Peters, Buckley, Statham, Pidsley, Samaras, Van Lord and Molloy2015), to identify DMRs between low and high contact groups. We only used the main effects model, given the large number of tests and the loose threshold applied for filtering significant regions. This computational tool uses a Gaussian kernel smoothing of DNAm across the genome, analyzing highly correlated CpG sites together as regions rather than individually, boosting statistical power.
For each of the above analyses, multiple test correction was performed using the Benjamini–Hochberg method (Hochberg & Benjamini, Reference Hochberg and Benjamini1990), using a threshold of .05 for candidate analyses, and .1 for the remaining, whole epigenome analysis performed in DMRcate. In all results, we report adjusted p values.
Epigenetic age deviation
The third analysis assessed the relationship between contact and epigenetic age deviation, the difference between epigenetic age, calculated using the Horvath epigenetic clock (Horvath, Reference Horvath2013), and chronological age. To assess whether postnatal contact related to epigenetic age deviation, we applied the following main effects and interaction models. In these models, age was not included as a covariate as it is accounted for in the calculation of epigenetic age deviation (which is the difference between biological and chronological age). These models do account for years of education as a proxy for socioeconomic status (Bradley & Corwyn, Reference Bradley and Corwyn2002) given prior associations between socioeconomic status and epigenetic age (Zannas, Reference Zannas2016):
Main effects model:
$$\eqalign{ \hbox{DNAm age deviation} =\; &{\rm \beta} _{\rm 0} + {\rm \beta} _{\rm 1} \times {\rm Contact} + {\rm \beta} _{\rm 2} \cr & \!\! \times {\rm Distress} + {\rm \beta} _{\rm 3} \times {\rm Sex} + {\rm \beta} _{\rm 4} \cr & \!\!\times {\rm Ethnicity} + {\rm \beta} _{\rm 5} \times {\rm Education}{\rm.}\,\,\,\,\,\,\,\,\,\,} $$
Interaction model:
$$\eqalign{ \hbox{DNAm age deviation} =\;& {\rm \beta} _{\rm 0} + {\rm \beta} _{\rm 1} \times {\rm Contact} + {\rm \beta} _{\rm 2} \times {\rm Distress} \cr & \!+ {\rm \beta} _{\rm 3} \times {\rm Sex} + {\rm \beta} _{\rm 4} \times {\rm Ethnicity} + {\rm \beta} _{\rm 5} \cr & \! \times {\rm Education} + {\rm \beta} _{\rm 6} \times {\rm Contact} \cr & \!\times {\rm Distress} + {\rm \beta} _{\rm 7} \times {\rm Contact} \times {\rm Sex} \cr & \!+ {\rm \beta} _{\rm 8} \times {\rm Contact} \times {\rm Ethnicity} + {\rm \beta} _{\rm 9} \cr & \!\times {\rm Contact} \times {\rm Education} + {\rm \beta} _{{\rm 10}} \cr & \!\times {\rm Distress} \times {\rm Sex} + {\rm \beta} _{{\rm 11}} \times {\rm Distress} \cr & \!\times {\rm Ethnicity} + {\rm \beta} _{{\rm 12}} \times {\rm Distress} \cr & \!\times {\rm Education}{\rm.}} $$
Results
Descriptive statistics and correlations among variables are displayed in Table 2 and Table 3. Genotype allele frequencies varied across different ethnic groups, consistent with previous literature (Bond et al., Reference Bond, LaForge, Tian, Melia, Zhang, Borg and Yu1998; Montag et al., Reference Montag, Brockmann, Lehmann, Müller, Rujescu and Gallinat2012; Pivac et al., Reference Pivac, Kim, Nedić, Joo, Kozarić-Kovacić, Hong and Muck-Seler2009). To further explore the association between contact and infant distress, we plotted this relationship (Figure 3). The moderate correlation can be partly explained by infants at the extreme levels of distress, with those infants showing the highest levels of distress falling into the high contact group, and those with the lowest levels in the low contact group.
Figure 3. Positive correlation between infant distress and continuous contact scores for final sample.
Table 2. Descriptive statistics and correlations among variables
Note: Bold indicates significance at p < .01. For all variables, n = 94, except for years of education (n = 93).
Table 3. Descriptive statistics for sex, ethnicity, and genotype groups
Candidate gene analyses
The goals of the following analyses were to replicate the original finding in the rat LG model of an association between postnatal contact and DNAm of the promoter region of NR3C1 in BECs, and to explore three candidate genes involved in neurobiological processes relevant to the biological embedding of postnatal contact.
NR3C1 pyrosequencing and array
To quantify DNAm at the exon 1F site, and this promoter region in general, seven CpG loci were pyrosequenced. The position of these sites and relation to those sites quantified by the 450k array is displayed in Figure 4, which shows that DNAm at the promoter sites is for the most part close to zero (in all cases, less than 5% methylated). Correspondence between the three sites that were quantified by both pyrosequencing and the array can be found in online-only supplementary Figure S.1.
Figure 4. (Color online) Pyrosequencing results for DNAm of NR3C1 for high and low contact groups. Sites within the red box are the analogous CpGs reported to demonstrate substantial differences between high and low licking and grooming behaviors in the original rat study of NR3C1 exon 17.
These overlapping sites were essentially invariable in our sample, partly due to the extremely low DNAm levels, and did not pass the variability filter we applied for narrowing down candidate CpG sites. Because we were interested in replicating the prior findings for this specific region of the promoter, we kept these CpGs for comparison to pyrosequencing, and for testing associations with contact. As shown in Figure 4, there were no significant differences between high and low contact at the promoter sites. Similarly, in the sites passing our variability filter quantified by the array, there were no differences (Figure 5a). In the interaction model, which allows the effect of contact to vary by infant distress, there were also no significant differences between high and low contact for any NR3C1 sites on the array or that were pyrosequenced (Figure 5b). There was adequate correspondence between values estimated through pyrosequencing and the 450k array, which is shown in online-only supplementary Figure S.1.
Figure 5. The 450k array results for DNAm of NR3C1, OPRM1, OXTR, and BDNF including estimates and confidence intervals for main effects and interaction models.
OPRM1, OXTR, and BDNF
We also tested associations between contact and DNAm at three candidate genes related to the neurobiological encoding and plasticity processes involved in early life contact. The sites quantified on the array within these gene candidates demonstrated comparable levels of variability to NR3C1, showing standard deviations from 2% to 7% methylation within variable probes (relative to 2% to 12% in variable NR3C1 sites). We did not find any significant associations between contact and candidate site DNAm in main effect models, testing main associations between contact and DNAm at candidate sites (Figure 5a). We also did not find support for significant interactions between contact and infant distress in terms of DNAm outcomes in these candidate regions (Figure 5b).
mQTLs analysis
Given the combined action of genetic factors and environment in shaping DNAm, we also used common genetic variants to identify mQTLs local to the candidate genes of interest, and tested whether associations of genetic variants with DNAm depended on contact or infant distress levels. Of the 11 variable DNAm sites of BDNF, the A allele significantly associated with increased methylation at two CpG sites (cg07238832 and cg10635145). One other mQTL was identified within OXTR, such that the intergenic A allele (rs7632287) related to a small decrease in DNAm at cg17036624. The results for statistical models testing whether genetic effects on DNAm are dependent on contact or on infant distress can be found in online-only supplementary Table S.1. There were no significant interactions between genotype and contact or infant distress.
DMR analysis
In this epigenome-wide association analysis, we tested whether patterns of DNAm across stretches of the genome related to high versus low contact. Using a false discovery rate threshold of 0.1, and filtered down results to only those DMRs containing at least two CpGs with at least a 5% differences in DNAm (in terms of β values) between high and low contact groups, DMRcate yielded a total of 44 CpGs within 5 DMRs, ranging from 2 to 22 CpGs (Table 4). Further details for each individual CpG can by found in the online-only supplementary Table S.2. DNAm within these identified regions is plotted for both high and low contact groups in Figure 6.
Figure 6. (Color online) Differentially methylated regions identified in epigenome-wide association analysis between high and low contact groups. TSS, transcription start site.
Table 4. Summary of differentially methylated regions
Epigenetic age deviation
This third analysis assessed the relationship between contact and epigenetic age deviation, accounting for differences in infant distress levels. The main effects model, testing for a main association between contact and epigenetic age deviation accounting for covariates, yielded no significant association. In contrast, the interaction model, allowing the relationship between contact and epigenetic age deviation to depend on infant distress levels, supported a significant interaction between contact group and infant distress (Figure 7; results in online-only supplementary Table S.3). Specifically, for the low contact group, there was a significant, negative association between distress and epigenetic age deviation, whereas for the high contact group, infant distress demonstrated a moderate and nonsignificant association with epigenetic age deviation. Higher infant distress in the context of low care thus related to lower epigenetic age relative to chronological age. Similarly, low infant distress in the context of low care related to older epigenetic age relative to chronological age. This result suggests that epigenetic age deviation is relevant to normative exposure to contact in early life, at least with consideration to infant distress.
Figure 7. (Color online) Relationship between infant distress and epigenetic age deviation for high and low contact groups.
Discussion
Infants thrive on close contact with caregivers in early postnatal life, when they are completely dependent on others for care and emotion regulation. Although the amount of contact during early critical periods has been the focus of intensive research using animal models, very little is known about the biological consequences of normative variations in early contact in humans. The current study applied a unique approach by directly addressing the role of frequency and duration of physical contact between caregiver and infant to quantify natural variation in human contact in the fifth week of postnatal life. We paid consideration to multiple levels of child factors that may be relevant to the association between contact and DNAm, including child genotype and distress levels. Although we did not find any differences between high and low contact groups in DNAm at candidate gene sites, including NR3C1, we did identify DMRs distinguishing between high and low contact. Moreover, we found that epigenetic age deviation was predicted by contact. In particular, high infant distress corresponded with reduced epigenetic age relative to chronological age, but only in conditions of low contact. Whether or not functional, these results suggest that experiences of contact from the earliest moments of postnatal life have lasting DNAm signatures.
Using a naturalistic daily diary approach, we assessed normative variations in infant contact, similar to the original LG studies in rats. In the LG studies, observations were made of maternal behavior every 4 min for 100-min observation periods over the first 10 days postpartum (Liu et al., Reference Liu, Diorio, Day, Francis and Meaney2000). The greatest differences in LG behavior was observed at Days 3–4, in which the high and low LG groups differed by 10% of observations. In contrast, for average daily contact over 4 days, our human high and low contact groups, differed by 8.93 hr, a difference of 37% in terms of daily contact. These groups were thus substantially different in postnatal contact during this early critical period.
Candidate genes
We first explored DNAm within four candidate genes involved in the neurobiological encoding of tactile contact. Given the link between DNAm and neuroplasticity, we reasoned that tactile contact levels may “get under the skin” by adjusting DNAm within these biological pathways that encode the effects of soft touch. DNAm at the first gene, NR3C1, and particularly at the exon 1F promoter, has been reported to relate to high versus low LG in the rat model, as well as early adversity in human studies in peripheral tissues (Turecki & Meaney, Reference Turecki and Meaney2016). With the present design, tailored to capture the human equivalent of LG in rats, we failed to replicate this previously reported association at exon 1F in humans. Moreover, we did not find DNAm at any sites located at NR3C1 to relate to contact.
It is noteworthy that to our knowledge it remains unclear whether the originally reported association between LG in rats and DNAm of exon 17 is robust, with one study reporting null results in a direct replication experiment (Pan, Fleming, Lawson, Jenkins, & McGowan, Reference Pan, Fleming, Lawson, Jenkins and McGowan2014). Although many human studies have targeted the analogous site in humans in relation to early adversity (Turecki & Meaney, Reference Turecki and Meaney2016), the effect sizes in these studies are very small, with the effects reported in most studies not exceeding differences in β values of 5% (Palma-Gudiel, Córdova-Palomera, Leza, & Fañanás, Reference Palma-Gudiel, Córdova-Palomera, Leza and Fañanás2015). An effect of at least 5% is a common, yet somewhat arbitrary, threshold for biological significance (although see Breton et al., Reference Breton, Marsit, Faustman, Nadeau, Goodrich, Dolinoy and Murphy2017, for discussion on interpreting small magnitude effects in DNAm studies).
Small effects at promoters, and in DNAm studies in general (Breton et al., Reference Breton, Marsit, Faustman, Nadeau, Goodrich, Dolinoy and Murphy2017), are related to the biological properties of DNAm. A single CpG site is present on each chromosome pair, or two times per cell; however, DNAm is quantified as an average across thousands or millions of cells in a biological sample, leading to a 0%–100% range of DNAm β values. Sites that are primarily methylated, or unmethylated, thus exhibit very high and very low β values, respectively. Promoter sites tend to be unmethylated, as seen in our data with the majority of our biological samples showing less than 5% DNAm at the NR3C1 promoter. The β values at the low and high extremes pose a particular problem for detecting effects, as the variability is highly compressed (Du et al., Reference Du, Zhang, Huang, Jafari, Kibbe, Hou and Lin2010). The three CpGs within the NR3C1 promoter that were quantified on the array did not surpass our variability filter. Moreover, all three CpGs were identified as invariable sites to be filtered for data reduction based on all publicly available 450k array data (Edgar, Jones, Robinson, & Kobor, Reference Edgar, Jones, Robinson and Kobor2017). Biologically significant effects of at least 5%, that are more likely to surpass the technical error of DNAm quantification are nearly impossible to observe at these invariable sites, which is why their removal might be preferable. In summary, despite the numerous studies reporting increased DNAm at NR3C1 exon 1F, it is worth considering whether these differences are meaningful, especially in comparison to the original, unreplicated finding reporting striking differences of 90% in rats (Weaver et al., Reference Weaver, Cervoni, Champagne, D'Alessio, Sharma, Seckl and Meaney2004).
Two candidate genes, OPRM1 and OXTR, were selected based on prior literature implicating their contribution to both reward and stress-attenuating neurobiological pathways that foster bonds through tactile exchanges with caregivers. We did not find any associations between DNAm at sites across these two genes and early postnatal contact. We also did not find any interactions between contact or infant distress and a SNP at an mQTL site identified in OXTR. At least in the current sample, and in BECs, it appears that DNAm of OPRM1 and OXTR is independent of variation in contact or infant distress, and that OXTR DNAm may be minimally driven by a cis genetic variant.
Prior work in OXTR DNAm in rats reported that low levels of LG related to lower DNAm of OXTR in blood but not in brain tissue (Beery et al., Reference Beery, McEwen, MacIsaac, Francis and Kobor2016). Given the suggestion that BECs are preferable surrogate tissues over blood for capturing DNAm variability relevant to the brain (Lowe et al., Reference Lowe, Gemma, Beyan, Hawa, Bazeos, Leslie and Ramagopalan2013; Smith et al., Reference Smith, Kilaru, Klengel, Mercer, Bradley, Conneely and Binder2015), our null result for OXTR may be consistent with this prior finding. Although there are studies suggesting that substance use corresponds with hypermethylation of OPRM1 in peripheral tissue (Cecil, Walton, & Viding, Reference Cecil, Walton and Viding2015), DNAm of OPRM1 in relation to postnatal contact has yet to be explored. The present study is the first to report that a social exposure, low versus high contact, does not relate to OPRM1 DNAm in peripheral BECs.
Animal work has indicated that DNAm of BDNF dynamically changes in key brain areas in response to learning (Lubin et al., Reference Lubin, Roth and Sweatt2008), and that DNAm at BDNF in the hippocampus and blood in rat models, and in cord blood in humans, relates to in utero exposure to the toxin bisphenol A (Kundakovic et al., Reference Kundakovic, Gudsnuk, Herbstman, Tang, Perera and Champagne2015). We did not find evidence in the current study for relations between early postnatal contact in relation to BDNF DNAm. It is possible that, if DNAm at BDNF in BECs is dynamic, no signatures of early postnatal contact were found because modifications at the level of DNAm were temporary. It is also possible that postnatal contact is just not a sufficiently potent exposure to have a lasting signature on BDNF DNAm. In humans, connections between prenatal exposures and DNAm quantified across the genome with the 450k array was found to depend on the BDNF val/met genotype (Chen et al., Reference Chen, Pan, Tuan, Teh, MacIsaac, Mah and Holbrook2015). Given our small sample, we only examined whether the association between maternal contact and DNAm at two mQTLs depended on the val/met genotype, rather than DNAm across the genome. We did not find any interactions between BDNF genotype and contact or infant distress. Overall, we did not find evidence that early postnatal contact relates to DNAm at this gene, regardless of BDNF genotype.
In summary, our analysis of differential DNAm in four candidate genes identified on the basis of their functional roles in encoding the reward and stress-reducing pathways of tactile contact produced null results. Our failure to find any significant associations between contact and DNAm at these candidate genes highlights concerns raised regarding the broader candidate gene literature, namely, whether scientists are capable of predicting relevant genetic targets given our rudimentary understanding of the boundless complexity of the biological processes at play in human development (Border & Keller, Reference Border and Keller2017; Duncan, Pollastri, & Smoller, Reference Duncan, Pollastri and Smoller2014). Although it is our opinion that informed candidate approaches are not entirely obsolete (Moore, Reference Moore2017), our null results, and the problem of invariability of DNAm at gene promoters, advises a broadening of scope in the epigenetics literature in psychology, which is heavily weighted toward popular candidate gene promoters.
Differentially methylated genomic regions
In contrast to our candidate gene approach, our analysis of DNAm assayed across the epigenome identified five regions in which differential methylation varied by high or low maternal contact during infants' fifth week of life. These results were yielded by an unbiased approach identifying region-specific clusters of sites that differ between high and low contact groups. The first DMR was within the lactase dehydrogenase A like 6A gene (LDHAL), encoding a protein involved in metabolic pathways. High contact was associated with higher DNAm. The second and third DMRs are intergenic (located in stretches of non–protein-coding DNA between genes), with the second showing high contact to be related to reduced DNAm at two sites, and the third showing high contact to be associated with greater DNAm consistently across 22 CpGs. The fourth DMR falls within the gene body of major histocompatibility complex, class II, DR beta 5 (HLA-DRB5), and demonstrates a flip between high and low contact groups, such that high contact is associated with greater DNAm at the first two sites, and then reduced DNAm at the last three sites. The HLA-DRB5 molecule plays a major role in the immunologic mediation of T cell responses. Finally, the fifth DMR is located just after the transcription start site of the protein-coding gene zinc finger AN1-type containing 2A (ZFAND2A), showing a substantial climb from low to high methylation, and with high contact associated with increased DNAm across the region. Although preliminary, these results suggest that normative variation in maternal contact in early postnatal life has an epigenetic signature in childhood.
Beyond the boost in power, our regional approach also had the advantage of yielding consistent differences between groups across a number of CpGs, filtered by a 5% β value threshold (for at least two sites). These consistent differences across stretches of DNA bolster confidence that our results are robust and possibly functional. However, it is important to caution that inferring functional effects of DNAm differences is not straightforward, given that these signatures were identified in BECs and that the relationship between DNAm and gene expression is complex. DNAm is linked to the regulation of gene expression, but may be either the cause or the consequence of expression levels, depending on the direction of the DNAm change and the location relative to the gene (Gutierrez-Arcelus et al., Reference Gutierrez-Arcelus, Lappalainen, Montgomery, Buil, Ongen, Yurovsky and Dermitzakis2013; Jones, Reference Jones2012; Wagner et al., Reference Wagner, Busche, Ge, Kwan, Pastinen and Blanchette2014). For instance, DNAm at gene promoters generally represses expression, but relationships between DNAm and expression at sites within gene bodies or between genes (i.e., intergenic), where we located our DMRs, are not well defined. With these caveats in mind, we speculate that these differential regions may reflect a component of the plasticity mechanisms involved in the biological embedding of early postnatal contact via cellular reprogramming.
Epigenetic age deviation
In our final analysis in this study, we assessed whether early life contact was associated with epigenetic age deviation, which is the difference between an epigenetic index of age and chronological age. We report no main association between maternal contact and epigenetic age deviation. However, the relationship between contact and epigenetic age deviation was contingent upon infant distress levels, such that for the high contact group, infant distress demonstrated a moderate, nonsignificant association with epigenetic age deviation, and for the low contact group, there was a significant negative association. Higher infant distress in the context of low care thus corresponded with lower epigenetic age relative to chronological age.
The combination of high distress in infants and low contact is not likely to be positive for child development, as infants during this early period of helplessness are completely reliant on caregivers. Acknowledging that interpretation of these results is preliminary given the current nascent literature applying the epigenetic clock to pediatric samples, we hypothesize that in this younger age range, age acceleration may not have the same implications as acceleration in adult samples, which relates to morbidity, mortality, and stress exposure (Chen et al., Reference Chen, Marioni, Colicino, Peters, Ward-Caviness, Tsai and Horvath2016; Marioni et al., Reference Marioni, Shah, McRae, Chen, Colicino, Harris and Deary2015; Zannas, Reference Zannas2016). Rather, in pediatric samples, age acceleration may be an indicator of progressive development through milestones, and deceleration an indicator of delays in maturation due to the costs of heightened stress. In other words, it is very tempting to speculate that the very same cellular processes that reflect morbidity and mortality at later development stages may also be involved in the acquisition of function and thriving in the young organism.
In terms of infant distress levels, this interaction is reminiscent of the notion of sensitivity to the environment, and specifically patterns in which infants higher in temperamental negative emotionality respond more strongly to social environments in terms of biological and behavioral outcomes (Bush & Boyce, Reference Bush and Boyce2016; Ellis, Boyce, Belsky, Bakermans-Kranenburg, & van IJzendoorn, Reference Ellis, Boyce, Belsky, Bakermans-Kranenburg and van IJzendoorn2011). However, it is important to note that infant distress measured in the current study is not likely to correspond with temperament, as differences in crying investigated prospectively appear to be transiently related to stress reactivity, and not predictive of more stable variation in emotions and regulatory capacity that composes temperamental measures (Barr & Gunnar, Reference Barr, Gunnar, Barr, Hopkins and Green2000; St. James-Roberts & Plewis, Reference St. James-Roberts and Plewis1996; White, Gunnar, Larson, Donzella, & Barr, Reference White, Gunnar, Larson, Donzella and Barr2000). Rather than infant distress marking sensitivity to high and low care, it is possible that those infants demonstrating the lowest levels of distress (which were all in the low contact group) are thriving with minimal care, whereas those infants high in distress are facing greater biological consequences in the face of low contact levels. If further investigations indicate the functional relevance of epigenetic age deviation as a link between early exposures and child social and psychopathology outcomes (Miller, Chen, & Cole, Reference Miller, Chen and Cole2009), then very early intervention aimed to increase the amount of contact between infant and caregiver may be critical for those infants with the highest levels of distress.
Implications of findings for psychopathology
Taken together, our results linking postnatal contact to DMRs across the genome and to epigenetic age deviation for higher distress infants suggests that the biological encoding of contact occurs early in life. This link between levels of contact and epigenetic signatures may implicate the specific stimulus of soft tactile contact in shaping development through plasticity mechanisms. Although further work is necessary to establish the specific brain and behavioral impacts of early postnatal contact, we speculate that the experience of touch might have particular relevance for broader social behavior and dysfunction. The ability to establish social connections relies upon an individual's history of social bonds, and particularly those with caregivers initiated by tactile contact. Deficits in forming close and secure social bonds is a common feature of psychopathology, and may be a precursor to problem behavior. For instance, social competence at age 4 predicts subsequent internalizing and externalizing symptoms (Bornstein, Hahn, & Haynes, Reference Bornstein, Hahn and Haynes2010) and the quality of mother–infant bonds at 2 weeks predicts child behavior problems at age 5 (Fuchs, Möhler, Reck, Resch, & Kaess, Reference Fuchs, Möhler, Reck, Resch and Kaess2016). Likewise, the establishment of close relationships due to enriched early experiences of contact may enhance resilience, as social support has been shown to buffer the effects of stress (Hostinar, Sullivan, & Gunnar, Reference Hostinar, Sullivan and Gunnar2014).
Influences of maternal contact
In humans there are a number of broader reasons that might drive differences in the amount of tactile contact infants receive from caregivers. First, our results suggest that evocative processes are at play, and specifically that infant distress elicits greater maternal contact. Second, external factors, such as the presence of other young children, or socioeconomic burden, can affect the amount of time caregivers can devote to their infants. At least in our sample, we did not find evidence that differences in years of education, a proxy for socioeconomic status, had a DNAm signal (Figure 2). We also did not find evidence that differences in education could explain the effect of infant distress on epigenetic age in the low contact group. Third, another possible driver of contact differences that we were unable to explore in the present investigation is caregiver emotional distress. It has been documented that mothers with depression touch and respond to their infants less (Field, Reference Field1984, Reference Field2010). It is worth exploring further if touch is one of the mechanistic drivers of associations between postnatal caregiver depression and distress and child DNAm (Cicchetti, Hetzel, Rogosch, Handley, & Toth, Reference Cicchetti, Hetzel, Rogosch, Handley and Toth2016), and social and regulatory outcomes (Choe, Olson, & Sameroff, Reference Choe, Olson and Sameroff2013; Pemberton et al., Reference Pemberton, Neiderhiser, Leve, Natsuaki, Shaw, Reiss and Ge2010; Székely et al., Reference Székely, Lucassen, Tiemeier, Bakermans-Kranenburg, van IJzendoorn, Kok and Herba2014). Fourth, more normative individual differences in attachment style and social traits may affect contact. Contact may be more frequent for caregivers who are, for instance, securely attached or highly agreeable, and thus find intimate exchanges more rewarding (Depue & Morrone-Strupinsky, Reference Depue and Morrone-Strupinsky2005).
Limitations
There are a number of limitations to the current study that deserve adequate consideration. First, our sample size of 94 is small; however, sampling from the extreme levels of the contact distribution enhanced statistical power, and we further attempted to maintain power by limiting the number of tests in our DNAm analyses wherever possible. It is still important for these results to be replicated. Second, this exploratory study did not include child behavioral or biological outcomes at follow-up beyond DNAm, which makes the interpretation of the association of epigenetic age deviation with infant distress in low contact environments difficult. Third, we quantified DNAm in BECs. Although saliva and BEC samples are proposed to be more relevant than blood cells for DNAm patterns in the brain (Lowe et al., Reference Lowe, Gemma, Beyan, Hawa, Bazeos, Leslie and Ramagopalan2013; Smith et al., Reference Smith, Kilaru, Klengel, Mercer, Bradley, Conneely and Binder2015), it is still not possible for us to make inferences about the identified DMRs in neural tissue.
Conclusion
Early postnatal tactile stimulation has shown robust associations with biological and behavioral outcomes in animal models, but the extent that the equivalent experience shapes human development is rarely addressed. Although we did not find evidence for expected associations between contact and DNAm at candidate genes related to the neurobiological encoding of tactile contact, we did find evidence overall that variation in tactile contact during early postnatal life associated with DNAm signatures in children. We report a particularly intriguing association between infant distress levels and epigenetic age deviation for those infants who experienced below average levels of contact, potentially indicating that epigenetic age deceleration, rather than acceleration, marks development risk in pediatric samples. Although these findings must be followed up with further longitudinal study of child neurobehavioral outcomes, we found initial support for the lasting biological embedding of postnatal contact at the level of DNAm. These results highlight the biological relevance of the experience of close and comforting contact critical to the formation of social bonds.
Supplementary Material
To view the supplementary material for this article, please visit https://doi.org/10.1017/S0954579417001213.
