Associations of peripheral blood DNA methylation and estimated monocyte proportion differences during infancy with toddler attachment style

ABSTRACT

Attachment is a motivational system promoting felt security to a caregiver resulting in a persistent internal working model of interpersonal behavior. Attachment styles are developed in early social environments and predict future health and development outcomes with potential biological signatures, such as epigenetic modifications like DNA methylation (DNAm). Thus, we hypothesized infant DNAm would associate with toddler attachment styles. An epigenome-wide association study (EWAS) of blood DNAm from 3-month-old infants was regressed onto children’s attachment style from the Strange Situation Procedure at 22-months at multiple DNAm Cytosine-phosphate-Guanine (CpG) sites. The 26 identified CpGs associated with proinflammatory immune phenotypes and cognitive development. In post-hoc analyses, only maternal cognitive-growth fostering, encouraging intellectual exploration, contributed. For disorganized children, DNAm-derived cell-type proportions estimated higher monocytes –cells in immune responses hypothesized to increase with early adversity. Collectively, these findings suggested the potential biological embedding of both adverse and advantageous social environments as early as 3-months-old.

KEYWORDS:

• Attachment styles
• DNA methylation
• monocytes
• early social environment
• APrON Study

1. Introduction

1.1 Attachment development and early life environments

Early life relationships are integral to who we will later become, both in our behaviors and in our biology. This early life impact is perhaps best exemplified by attachment theory (Bowlby, 1979). Attachment is a learned behavioral pattern that children display with their parental caregivers, usually mothers, in which children variously can (or cannot) use their caregiver as a secure base and safe haven from stress. It represents a brain-based biological evolutionary system promoting infant survival (Ainsworth et al., 1978; Chisholm, 1996) through strategies that maintain proximity of infants to caregivers (Ainsworth et al., 1978; Shah et al., 2010). Four attachment styles have been identified: secure (Type B), insecure avoidant (Type A), insecure ambivalent/resistant (Type C), and disorganized (Type D) (Ainsworth et al., 1978; Main & Solomon, 1986). The infant is expected to be securely attached (Type B) and able to use their caregiver as a secure base for exploration and safe haven from stress when the caregiver is sensitively responsive (Ainsworth et al., 1978).

In contrast, insecure avoidant attachment (Type A) is observed in infants whose parents are withdrawn or under-stimulating; in this context, infants are observed to inhibit their displays of emotion, comfort, and proximity seeking to parents. Insecure ambivalent/resistant attachment (Type C), often referred to as anxious attachment in adulthood, is observed in infants whose parents are unpredictably, but intermittently responsive when comfort and safety is needed; in this context, infants are observed to exaggerate their emotional displays (e.g. anger, fear, and desire for comfort) in their relationships with parents (Sroufe, 2005; Sroufe et al., 1999). Disorganized attachment (Type D) is observed in infants whose parents provide fearful or fear-provoking stimuli or who present with unresolved trauma (Main & Hesse, 1990); in this context, infants are observed to present contradictory, misdirected, or freezing and/or fear behaviors in interactions with parents. These behaviors are thought to represent an inability to organize around their caregiver’s distressing behavior and occur more often in stressful and abusive family settings (Cyr et al., 2010; Madigan et al., 2007; Van Ijzendoorn et al., 1999).

Attachment development occurs within the first years of life, as the child gains both biological and cognitive maturity. By approximately 8–12 months children reach the “clear-cut attachment” phase, which is preceded by the “attachment-in-the-making” phase from approximately 2–8 months (Bowlby, 1969, 1982). It is during this phase that infants learn from their environment and influential caregiver interactions. Though attachment is still developing during this time, infants’ attention and behaviors show differences indicative of future attachment patterns (Conradt & Ablow, 2010; Sherman et al., 2015). Once an attachment style has developed, this internal working model of relationships is usually maintained, though flexible to some modification, into and throughout adulthood (Carlson et al., 2004; Chopik et al., 2019; Fraley, 2019; Groh et al., 2014; Hazan & Shaver, 1994; Khan et al., 2020; Sroufe, 2005; Waters & Waters, 2006). Many aspects of attachment are also echoed in our understanding of epigenetics, specifically DNA methylation (DNAm) (Aristizabal et al., 2019; Merrill et al., 2019), a molecular marker that may be an ideal avenue to examine physiological differences during attachment formation in very early life.

1.2 Evidence for the link between DNA methylation and attachment style

DNAm is the covalent addition of a methyl group to DNA Cytosine-phosphate-Guanine (CpG) dinucleotides. It is relatively stable, mitotically heritable, potentially associates with gene transcription, is variable during the first 5 years of life, and does not change the DNA sequence (Aristizabal et al., 2019; Jones et al., 2018; Nagy & Turecki, 2012; Pérez et al., 2019). DNAm allows researchers to study how early life environments can “get under the skin” (Boyce & Kobor, 2015) to potentially result in dynamic variation in physiology, behavior, and health outcomes – ideal for studying as a potential biomarker of early attachment experiences. Understandably then, there have been a plethora of studies examining maternal sensitivity links with offspring DNAm in both model organisms and humans (Dall’ Aglio et al., 2020; Mcgowan et al., 2011; Moore et al., 2017; Provenzi et al., 2020), as well as several studies examining the association between DNAm and attachment, specifically (Darling Rasmussen & Storebø, 2021; Robakis et al., 2020).

While the majority of these maternal sensitivity and attachment studies found associations using candidate gene approaches, especially the oxytocin (OXTR) and glucocorticoid (NR3C1) receptors (Darling Rasmussen & Storebø, 2021; Provenzi et al., 2020), one epigenome-wide study found infant attachment style may associate with up to 12% of DNAm variation in their sample (Garg et al., 2018). Additionally, DNAm is vital to cell type identity, and can, thus, be used to estimate cell type proportions (Houseman et al., 2012). Cell type proportion differences have been associated with differences in early life environment, especially exposure to social stressors (Esposito et al., 2016; Lam et al., 2012; Miller et al., 2011). Current evidence suggests a replicable and robust association of maternal sensitivity, and potentially attachment style, with DNAm, especially at genes in developmental pathways and as an intermediary between early adversity and DNAm (Darling Rasmussen & Storebø, 2021; Provenzi et al., 2020). Thus, while the presence of an association of caregiving to children’s DNAm during early life is consistent across the literature, to date, most of the field has focused on candidate gene investigations of psychopathology as well as the experiences of adversity and extreme caregiving environments in relation to characterizing parent-child relationships (Darling Rasmussen & Storebø, 2021; Robakis et al., 2020). However, there is currently a gap in our understanding of biological correlates in human attachment development, a period when both attachment and DNAm could be hypothetically reactive to the parental care environment, within the context of typical maternal caregiving.

1.3 Current study

We examined epigenome-wide DNAm associations with attachment style within a low-risk community sample of mothers and infants. Specifically, we assessed the connection between blood DNAm in infant offspring during the sensitive “attachment-in-the-making” phase of attachment development at 3-months with attachment style at 22-months. While the research strongly supports a link between early life context, especially maternal care, and attachment development, these data presented a unique opportunity to examine potentially lasting biological associations during this development period that had heretofore remained a “black box in attachment theory” (Bosmans et al., 2020).

First, we hypothesized that DNAm levels at CpGs in 3-months-olds would associate with toddler attachment style. Second, we hypothesized that birth characteristics and parent-child interaction quality may be effect modifiers that contribute to DNAm associations at 3-months-old with toddler attachment style.

2. Materials and methods

2.1. Participants

2.1.1. Alberta pregnancy outcomes and nutrition (APrON) cohort

The APrON cohort is a Canadian longitudinal cohort of 2139 families with children that enrolled women between 2009 and 2013 during their first or second trimester of pregnancy. Mothers voluntarily consented both during enrollment and at the time of collection to provide DNA samples from their children at 3 months of age. A representative subsample of children from mothers who reported no smoking, alcohol consumption, or steroid use during pregnancy donated blood samples (Table 1). Sample demographics were obtained through mothers’ survey self-reports or through hospital birth records. Infants were mostly white and represented equal numbers of males and females (Table 1). In total, data and samples from 93 dyads were employed in this analysis. Approximately half of the dyads were classified as securely attached (Type B) and the observed percentages aligned with the literature on attachment style ratios in normative, community populations (Mesman et al., 2016) (Table 1). There was a proportionally equal distribution of sexes, ethnicities, and ages across the attachment styles.

Table 1. Participant descriptive statistics.

	# Individuals Available Data	n	%	Mean	sd	Range
Secure Attachment (B)	93	47	50%	-	-	-
Insecure-Avoidant Attachment (A)	93	21	23%	-	-	-
Insecure-Ambivalent Attachment (C)	93	14	15%	-	-	-
Insecure-Disorganized Attachment (D)	93	10	11%	-	-	-
Female Infants	93	50	54%	-	-	-
Caucasian Infants	92	81	88%	-	-	-
Infant Age at Sample Collection (Weeks)	93	-	-	12.70	1.05	5.14
Infant Age at SSP Measurements (Months)	93	-	-	21.27	3.51	17
Infant Gestational Age at Birth (Months)	93	-	-	39.37	1.44	8.14
Infant Birth Weight (kg)	93	-	-	3425.41	520.60	2644
Maternal Age at Birth (Years)	93	-	-	31.08	4.02	22
Maternal Postpartum Depression Score at 3 Months	90	-	-	5.18	4.18	19
Mother Married at 3 Months	93	86	92%	-	-	-
Combined Household Income of $70 K or more	92	56	61%	-	-	-
Maternal Caregiver Total PCITS Score at 6 Months	86	-	-	37.07	4.66	23
Infant Child Total PCITS Score at 6 Months	86	-	-	16.21	3.16	14

2.2 Measures

2.2.1 Strange situation procedure (SSP)

At 22 months of age, maternal-infant pairs completed the SSP, the gold standard assessment for attachment styles in infancy (Ainsworth et al., 1978). This procedure contains eight brief episodes designed to activate infant’s attachment behaviors by evoking mild levels of stress in children through a series of mother–child separations and reunions, and interactions between the infant and a “stranger” (a research confederate), where the child’s behaviors were observed through a two-way mirror and video-recorded for coding by trained and reliable coders. All SSP videotapes were coded for styles of attachment using standard categories of secure (Type B) and insecure with subtypes avoidant (Type A) and ambivalent/resistant (Type C), and disorganized (Type D) (Ainsworth et al., 1978; Main & Solomon, 1986). A random 15% of the recordings were also re-coded by an independent expert coder. Cohen’s kappa for interrater reliability was .73. Due to the relatively small group sample sizes of insecure category subtypes the groups were dichotomized into infants classified as secure (B) versus not secure (A, C, and D) and disorganized (D) versus not disorganized (A, B, C). Disorganized (D) behavior was coded using the Scale for Disordered/Disoriented Infant Behavior (1 = low, 5 = moderate, 9 = extreme; with 5 = threshold) (Mary Main & Solomon, 1990; Mary Main & Stadtman, 1981), which assesses the degree of disorganized behavior regarding conflict (e.g. aggressive outbursts) and/or disruptive behaviors (e.g. immobilized, disoriented, sudden disordered activities, uninterpretable noises or movements) during the SSP.

2.2.2 Parent–child interaction teaching scale (PCITS)

The PCITS is a well-validated observational tool for assessing parent-child interaction quality before 36 months of age (Oxford, M. & Finlay, D, 2013) and conducted with APrON mother-child dyads when the child was 6 months old (Letourneau et al., 2018). The 76-item (scored 0 = no, 1 = yes) tool assesses interaction quality in four maternal subscales: sensitivity to cues, response to child’s distress, social-emotional growth fostering, and cognitive growth fostering, as well as two child subscales: clarity of cues and responsiveness to caregiver. Subscale scores are summed into caregiver, child, and overall totals with higher totals indicating more positive parent–child interaction quality.

2.2.3 Maternal edinburgh postnatal depression scale (EPDS)

Maternal postnatal depressive symptoms were assessed when the infants were approximately 3-months-old using the EPDS. This self-report instrument consists of 10-items rated on a 4-point scale that has been well-validated (Affonso et al., 2000).

2.3 Biological sample and data processing

Overall, 93 samples and 96,339 variable probes were included in the analyses (Table 2). Probe variability was determined using an inter-quantile range filter subset to probes where the DNAm β value varies by at least 5% across samples in the 5th and 95th percentile (Edgar, Jones, Robinson, et al., 2017). This reduced the testing space by removing probes that did not vary across individuals in our sample to better meet the Benjamini-Hochberg False Discovery Rate (FDR) control method assumptions of equal significance likelihood (Benjamini & Hochberg, 1995; Korthauer et al., 2019). For additional biological sample handling, data processing, cell type proportion estimation, and estimated cell type principal component calculation, please see Appendix.

Table 2. Cohort blood samples and DNA methylation probes were normalized and processed for quality control. A table of participant samples and Cytosine-phosphate-Guanosine (CpG) probes from the Illumina 450 K Beachip array in each step of preprocessing. Upon importing raw IDAT files, control probes were clustered to determine if any sample had abnormal values (Fortin et al., 2017). One sample had control probes that failed and was removed. Samples would also have been removed as outliers if failing additional outlier methods: detectOutlier, pcout, locFDR (Du et al., 2018; Efron, 2007; Efron et al., 2001; Hannum et al., 2013; Pidsley et al., 2013). Using these criteria, no additional samples were removed as outliers. Bad probes, defined as having a low detection p-value (≤1 × 10^–16), fewer than 3 reads, or an NA in more than 1% of samples, were removed. Probes measuring SNPs or X and Y chromosomes (Price et al., 2013) were removed. Poorly designed probes, defined as polymorphic probes measuring a CpG that contains a difference in genetic sequence (Fortin et al., 2017) or multiple binding probes predicted to non-specifically bind (Pidsley et al., 2016; Price et al., 2013) were removed. Sex mismatches were discovered by employing K-means clustering with probes associated with sex and comparing this to the reported sex. One sample was removed as a sex mismatch, and one sample of the technical replicate pair was removed. Then, probes which were invariable in the cohort, a probe with <5% methylation β-value range between the 5th and 95th percentile, were removed. Finally, only samples with both DNAm at 3-months-old and Strange Situation Procedure (SSP) attachment style data at 22-months-old were included in the analysis for a total of 93 samples and 96,339 variable probes included in the analysis.

450 K Blood Samples QC Step	Samples Removed	Probes Removed	Samples Remaining	Probes Remaining
Sample Outlier	1	0	116	485,577
Bad Probes	0	8,323	116	477,189
SNP Probes	0	65	116	485,512
XY Probes	0	11,051	116	466,138
Polymorphic CpGs	0	18,697	116	447,441
Multiple Binding Probes	0	37,499	116	409,942
Sex Mismatch	1	0	115	409,942
Technical Replicates Removed	1	0	114	409,942
Variability Filter	0	313,603	114	96,339
Available SSP Data	21	0	93	96,339

2.4 Cell type proportion analysis

2.4.1 Estimated cell type proportion association

Cell type differences are a major driver of inter-individual DNAm variability (Jones et al., 2018). After estimating cell type proportions, we tested for potential DNAm predicted cell type differences associated with attachment style controlling for sex, age, and ethnicity in an ANCOVA model. As they represent the two most extreme phenotypes, both secure (Type B) versus insecure (Non-Type B) attachment style and disorganized (Type D) versus organized (Non-Type D) attachment style comparisons were tested. Because the constituency of the attachment styles is unequal, for significant associations 10,000 randomized permutations were run with attachment styles classification shuffled and assigned randomly (ter Braak, 1992). The resulting value is the permutation p-value that represents the likelihood that comparisons were observed due to chance. This comparison was also replicated with cell proportions estimated from multiple reference panels (Figure 1).

Figure 1. Higher predicted monocyte proportion in infant blood of toddlers with a disorganized attachment pattern than toddlers with a secure attachment pattern across all references used. Boxplots of toddler attachment pattern organization (Type D vs Non-Type D) and estimated monocyte proportion in infant buffy coat blood samples. Each plot shows estimated monocyte proportions calculated with a different cell type reference: (A) Identifying Optimal Libraries (IDOL) adult reference (Koestler et al., 2016), (B) umbilical cord blood reference (Gervin et al., 2019), (C) mixed reference combining both IDOL and cord blood, and (D) mixed reference combining both IDOL and cord blood, but also constraining predicted nucleated red blood cell proportion, a cell type only found in umbilical cord and newborn blood, to zero (Hermansen, 2001). The bar and asterisk indicate a significant difference (*p < .05 and **p < .01) between Type D and Non-Type D attachment patterns. Color refers to each of the four measured attachment styles from the Strange Situation Procedure (Ainsworth et al., 1978; Mary Main & Solomon, 1990).

2.5 DNA methylation association analysis

2.5.1 Biological DNA methylation effect size (Δβ) calculation

Change in DNAm value (or β value, indicating the proportion of methylated DNA molecules from zero to one) is the DNAm biological effect size, referred to as Δβ (Jones et al., 2018). Δβs were calculated by subtracting the average β of Type A,C,D attachment cases from the Type B group on a per-CpG basis. The technical Δβ threshold of ǀ≥0.03ǀ was determined to likely be greater than technical noise using the root mean squared error of the technical replicates (RMSE = 0.015) and the maximum standard error across all 96,339 variable CpGs for all individuals (SE = 0.027). To obtain a clearer understanding of biological differences in CpGs associated with attachment styles, significant CpG’s Δβs were also calculated for each style’s comparison by subtracting the mean β value of one type from the other.

2.5.2 Epigenome-wide association study (EWAS) analysis

For the a priori candidate analysis, biologically relevant candidates were chosen from recent publications reporting DNAm associations on attachment (Darling Rasmussen & Storebø, 2021) and early maternal caregiving (Provenzi et al., 2020). Of the potential candidates from these findings, 132 variable CpGs (13 genes) were tested a priori (Supplemental Table 1).

The exploratory EWAS was conducted by running an ANCOVA model on β values of all variable probes (n = 96,339) (Supplemental Table 2). Both candidate and exploratory approaches employed the same statistical model. The variable of interest was secure (Type B) versus insecure (Non-Type B) attachment styles. Model construction was tested for multicollinearity (Supplemental Figure 2) and relations were assumed to be linear and normally distributed. Variables (sex, age, ethnicity, and three PCs cell type proportions) that contribute to DNAm differences (Jones et al., 2018) were included as covariates in the model:

$CpG\mathrm{\setminus }AttachmentPattern+Age+Sex+Ethnicity+PredictedCellProportionPCs+\epsilon$

Any cell types found to associate with attachment style differences (monocytes) were excluded from these cell type proportion PCs.

All regressions were multiple-test corrected using the Benjamini-Hochberg (BH) FDR control method on the full testing space (Benjamini & Hochberg, 1995). Medium-confidence statistical threshold of FDR<.10 and a technical effect size threshold of Δβ≥|0.03| were used to identify significant associations with attachment style. All reported associations were from homoscedastic models as determined by Bartlett’s test. To examine if any detected CpGs were sensitive to the exclusion of cell types (monocytes) from the PCs in the first model, a second EWAS model including the PCs of all cell type proportions was run and results compared (Supplemental Table 3).

Table 3. Twenty-six CpGs were identified to associate with a secure attachment style. Represented are the details of each CpG discovered to be significantly associated with secure attachment style in infant buffy coat blood, including chromosome, map position, ∆β (Type B – Non-Type B). The FDR and p-value of each site of the primary model (not adjusted for estimated monocyte proportion) are listed, while CpGs that were not significant in the secondary model adjusting for all estimated cell type proportions (including monocytes) were marked with an asterisk. This is followed by descriptive information on each CpG: UCSC genome browser gene association, if one is present; genetic region; if there is a highly associated SNP, or methylation quantitative trait loci (mQTL), within 50kb of the CpG; and, if the CpG is part of a co-methylated region (CMR) (Gatev et al., 2020) with the number of CpGs in that CMR. CMR lengths marked with a dagger indicate there is another medium-confidence CpG also within that same CMR. Finally, the correlation and variability of each CpGs DNAm β value between blood and averaged over three areas of brain tissue (Broadman areas 10, 20, and 7) as provided by the Blood-Brain Epigenetic Concordance (BECon) tool (Edgar, Jones, Meaney, et al., 2017).

CpG	Chr	Position	∆β	p-value	FDR	Gene	Gene Feature	CpG Island	mQTL SNP	CMR Length	Average Brain-Blood Corr	Average Brain Var
cg12691488	1	243,053,673	−0.0361	1.33E-07	0.0027	-	intergenic	chr1:243,053,597–243,054,072	-	-	0.3873	0.1832
cg27049094	2	75,067,716	−0.0314	1.24E-05	0.0138	HK2	intragenic	-	-	-	0.3863	0.0349
cg07586008	19	49,636,535	0.0475	2.59E-05	0.0183	PPFIA3	intragenic	chr19:49,636,224–49,636,694	-	3†	0.1578	0.0779
cg18032014	7	95,402,431	0.0414	0.0001	0.0308	DYNC1I1	intragenic	chr7:95,401,691–95,402,432	-	4†	−0.1314	0.1817
cg18984103	14	104,163,025	0.0646	0.0001	0.0353	KLC1	intragenic	chr14:104,165,051–104,165,278	-	-	0.1157	0.0858
cg20709110	22	19,946,873	0.0371	0.0001	0.0353	COMT	intragenic	-	-	-	0.0931	0.1083
cg14640466	8	135,056,370	−0.0305	0.0001	0.0358	-	intergenic	-	-	-	−0.0059	0.1374
cg24913825	2	198,134,173	−0.0313	0.0002	0.0405	-	intergenic	-	-	-	−0.1196	0.0430
cg04742282	6	28,911,873	0.0301	0.0002	0.0405	LINC01556	intragenic	-	-	13†	0.1255	0.0433
cg17889682	7	95,402,733	0.0651	0.0002	0.0409	DYNC1I1	intragenic	chr7:95,401,691–95,402,432	-	4†	−0.3412	0.1450
cg06421590	4	3,043,752	−0.0359	0.0002	0.0416	-	intergenic	chr4:3,043,018–3,043,547	-	4†	0.4510	0.1920
cg01807241	4	3,043,527	−0.0368	0.0002	0.0444	-	intergenic	chr4:3,043,018–3,043,547	-	4†	−0.0471	0.0153
cg07273376	11	23,919,652	−0.0301	0.0003	0.0473	-	intergenic	-	-	-	0.0627	0.1790
cg10281478	7	95,402,528	0.0314	0.0003	0.0490	DYNC1I1	intragenic	chr7:95,401,691–95,402,432	-	4†	0.1559	0.0320
cg25005112	6	28,911,938	0.0318	0.0004	0.0591	LINC01556	intragenic	-	-	13†	0.1755	0.1026
cg22744079	4	155,410,860	0.0360	0.0005	0.0635	DCHS2	intragenic	chr4:155,410,623–155,413,393	-	3	0.3961	0.1404
cg09812376	17	74,270,190	−0.0875	0.0005	0.0675	QRICH2	three plus	-	-	-	0.3147	0.0878
cg18443437	5	3,326,731	−0.0607	0.0005	0.0686	-	intergenic	chr5:3,324,644–3,324,973	-	-	0.4137	0.1647
cg01265860	21	36,256,316	0.0303	0.0006	0.0720	RUNX1	intragenic	chr21:36,258,952–36,259,472	-	-	−0.0696	0.0798
cg17292758	19	49,636,594	0.0354	0.0006	0.0731	PPFIA3	intragenic	chr19:49,636,224–49,636,694	-	3†	0.3176	0.0917
cg11070193	1	26,527,870	0.0678	0.0007	0.0784	CATSPER4	intragenic	-	-	3	−0.2147	0.1068
cg04487205	17	36,997,563	0.0344	0.0008	0.0841	C17orf98	promoter	chr17:36,997,448–36,997,661	-	7	0.3627	0.1002
cg13293729*	6	28,911,902	0.0363	0.0009	0.0900	LINC01556	intragenic	-	-	13†	0.0147	0.0263
cg12115302*	12	30,323,676	0.0337	0.0010	0.0908	-	intergenic	chr12:30,322,763–30,323,555	-	-	−0.3578	0.0473
cg14616234*	9	84,305,683	−0.0314	0.0010	0.0917	-	intergenic	chr9:84,302,139–84,305,001	-	-	0.1186	0.1463
cg03592051*	16	27,899,678	−0.0349	0.0010	0.0988	GSG1L	promoter	-	-	-	0.4490	0.1410

Identified CpGs were annotated to genes, CpG islands, genomic regions, and map locations with the UCSC genome browser (Price et al., 2013). To connect the CpGs to potential biological functions, the GeneCards (www.genecards.org) database and both PheWAS (https://atlas.ctglab.nl/PheWAS) and EWAS Atlases (http://bigd.big.ac.cn/ewas) were used. Chromatin state was gaged via ChromeHMM (Ernst & Kellis, 2017) based on histone modifications. Additionally, the Blood-Brain Epigenetic Concordance (BECon) tool (https://redgar598.shinyapps.io/BECon/) was used to investigate the concordance and variability of CpGs in blood and brain regions (Broadman areas 7, 10, 20) (Edgar, Jones, Meaney, et al., 2017).

2.5.3 Methylation quantitative trait loci (mQTL) analysis

Single nucleotide polymorphism (SNP) allelic variants can influence DNAm patterns of nearby CpGs, referred to as mQTLs (Bell et al., 2011; Cheung et al., 2017; Gertz et al., 2011). Therefore, significant CpGs were also analyzed against SNPs to examine potential contributions of the genomic context. To increase mQTL discovery likelihood, only variable SNPs within a 50kb window from the CpGs were included in the analysis using a CpG ~ SNP ANOVA model (Supplemental Table 4) (Islam et al., 2019; Smith et al., 2014).

Table 4. Immune system pathways were the top identified in a KEGG pathway enrichment analysis for secure attachment style exploratory EWAS results. Rank ordered values from the pathway terms that passed a medium-confidence FDR<.10 from the methylglm enrichment analysis, accounting for the number of CpGs in each term (Ren & Kuan, 2019). Provided p-values were from the secondary EWAS model accounting for all estimated cell type proportions (including monocytes). This analysis used the KEGG (Kyoto Encyclopedia of Genes and Genomes) database resource to discover pathway associations (Ganter et al., 2005). The number of genes in each term, p-value, and FDR value were included.

KEGG ID	Description	Gene Set Size	p-value	FDR
4660	T cell receptor signaling pathway	108	6.45E-11	2.84E-09
4062	Chemokine signaling pathway	189	4.08E-05	0.0007
4630	JAK-STAT signaling pathway	155	4.84E-05	0.0007
4380	Osteoclast differentiation	128	0.0001	0.0012
4060	Cytokine-cytokine receptor interaction	265	0.0003	0.0023
4650	Natural killer cell mediated cytotoxicity	136	0.0006	0.0044
4670	Leukocyte transendothelial migration	116	0.0019	0.0120
5142	Chagas disease	104	0.0025	0.0136
4110	Cell cycle	124	0.0031	0.0151
5160	Hepatitis C	134	0.0055	0.0242
4810	Regulation of actin cytoskeleton	213	0.0166	0.0613
4360	Axon guidance	129	0.0167	0.0613
4514	Cell adhesion molecules	133	0.0262	0.0888
4141	Protein processing in endoplasmic reticulum	165	0.0293	0.0920

2.5.4 Co-Methylated region (CMR) analysis

Because DNAm patterns at CpGs located close to one another are highly correlated, and may act together, we examined if associated CpGs were within CMRs as determined using the Co-Methylation with genomic CpG Background (CoMeBack) method (Gatev et al., 2020). 21,803 CMRs were constructed from this cohort and interrogated here to define potential regions of co-methylation for CpGs significantly associated with future secure attachment development.

2.5.5 Gene ontology (GO) analysis

A GO analysis was performed on EWAS results using methylGSA (Ren & Kuan, 2019) accounting for CpG probe number length bias (Mi et al., 2012). The secondary EWAS model accounting for all cell types was chosen in order to not overrepresent the possible effect of differing monocyte proportions in potential pathway associations. Ranked CpGs with a p < .0001 were included in analyses of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, GO, and Reactome databases, as was implemented in recent EWAS literature (Croft et al., 2011; Ganter et al., 2005; Mooney et al., 2020; Ren & Kuan, 2019).

2.5.6 Contribution analysis

We also interrogated the potential influence of maternal-infant relationships and birth characteristics known to associate with DNAm on significant CpGs (Mayne, 2012). To investigate these contributions, a method used previously in a study examining DNAm was implemented (Roberts et al., 2018). Percentage contribution for each variable of interest was calculated as:

CpG ~ β_base Attachment Style + Age + Sex + Ethnicity + Predicted Cell Proportion PCs

CpG ~ β_adj Attachment Style + Contributor + Age + Sex + Ethnicity + Predicted Cell Proportion PCs

% Contribution = 100 $\left[\frac{{\beta }_{base}-{\beta }_{adj}}{{\beta }_{base}}\right]$

3. Results

3.1 Higher monocyte proportions at 3-months were associated with toddler type D classification

While there was no difference in infant estimated monocyte proportions between future secure (Type B) and future non-secure attachment styles (Type A,C,D) (1,87) F = 1.685 p = 0.2010, 90% CI [0.00, 0.32], infants who would later be categorized as disorganized (Type D) had a significantly higher predicted proportion of monocytes at 3-months-old than those with later organized attachment styles (A,B,C) (Figure 1) (1,87) F = 7.838 p = 0.0063, 90% CI [0.12, 0.48] Cohen’s f = 0.30. A randomized permutation test (p = 0.011) indicates this finding was not due to chance or imbalance in the groups. Additionally, this significant infant cell type difference was retained between future Type D and future organized attachment styles when using or including an umbilical cord blood reference for cell type estimation ((1,87) F = 6.365 p = 0.0135, Cohen’s f = 0.27) (Figure 1). Future Type D infants were consistently estimated to have a higher estimated infant monocyte proportion than future organized attachment infants. Due to the unique importance of cell type on DNAm studies, this estimated cell type proportion difference was accounted for in the future secure (Type B) versus non-secure (Type A,C,D) attachment style EWAS.

3.2 Twenty-six medium-confidence CpGs were associated with toddler type B classification

With roughly equal groups of infants who would later be categorized as secure (Type B) and non-secure (Type A,C,D), we examined a model to explore potential DNAm precursors to future secure attachment style development. To operationally define medium-confidence CpGs, we adopted a statistical threshold of FDR<0.10 and a technical effect size threshold of |Δβ≥.03| reminiscent of thresholds used in EWAS literature (Bush et al., 2018; Essex et al., 2013), where none of the CpGs included in the a priori candidate tests passed; but, two CpGs in TNF passed the FDR threshold alone (Supplemental Table 1). Using 96,339 variable CpGs in the exploratory EWAS, we ran an ANCOVA model of attachment style on DNAm (Figure 2; Supplemental Table 2) where 26 CpGs were significantly associated with future secure attachment style (Type B vs Non-Type B) (Figures 3–5; Table 3). Four CpGs met our medium-confidence criteria only in the primarily model, which did not account for monocyte cell type proportion, indicating their association is sensitive to estimated monocyte proportion (Table 3; Supplemental Figure 3; Supplemental Table 3).

Figure 2. Exploratory EWAS model results indicate a DNAm signal of secure attachment style. A p-value histogram and volcano plot of the primary EWAS model results. This model interrogates differences between Type B and Non-Type B attachment styles adjusting for sex, age, ethnicity, and three estimated cell type proportion principal components calculated excluding predicted monocyte proportions. This is because predicted estimated monocyte proportion correlated with attachment style. (A) A right-skewed histogram of the p-value distribution from each CpG result in the exploratory EWAS. (B) A volcano plot depicting medium-confidence CpGs identified as being significantly associated with toddler secure attachment style at an FDR<0.10 and ǀΔβ≥0.03ǀ, which are represented by the solid and dashed lines, respectively.

Figure 3. Sixteen medium-confidence CpGs’ DNAm differ with secure versus insecure attachment styles. Boxplots of secure versus insecure attachment style (Type B vs Not Type B) DNAm β values show differences for medium-confidence CpGs. β values range from zero to one and are on a variable y-axis scale of 0.5. Color refers to attachment style. Eight CpGs exhibited decreased β values and eight exhibited increased β values with secure attachment.

Figure 4. Seven medium-confidence CpGs associated with secure attachment style are within co-methylated regions (CMRs) with multiple significant CpGs. Three plots depicting seven medium-confidence CpGs associated with secure attachment within a CMR. Co-methylated regions were grouped if CpG β values were correlated by 40% and within 1kb of the next probe accounting for the interspersed density of background CpG sites within 400bp (Gatev et al., 2020). At the top of each plot is the chromosome and genomic position of each CMR, under which is the relative position of each CpG. CpGs marked with an asterisk were identified as medium-confidence CpGs associated with secure attachment style. Beneath those are gene associations based on the UCSC genome track (Kent et al., 2002) and the NCBI genome track (Agarwala et al., 2018), followed by the chromatin state of the region as determined by ChromeHMM using profiled histone marks (Ernst & Kellis, 2017). Finally, a line plot of DNAm β value, on a scale from zero to one, at each CpG in the CMR. Lines indicate average values, the surrounding opaque bars display confidence intervals, and colors indicate secure (Type B) and insecure (Non-Type B) attachment style. (A) CMR within LINC01556 and a tRNA cluster containing three medium confidence CpGs and 13 CpGs in total. (B) Intergenic CMR containing two medium-confidence CpGs and four CpGs in total. (C) CMR within PPFIA3 containing two medium-confidence CpGs and three CpGs in total.

Figure 5. Three medium-confidence CpGs associated with secure attachment style are within the four CpG DYNC1I1 co-methylated region (CMR). A plot which at the top is presented a map of chromosome seven and the CMR’s genomic position, under which is the relative position of each CpG. CpGs marked with an asterisk were identified as medium-confidence CpGs associated with secure attachment style. Co-methylated regions were grouped if CpG β values were correlated by 40% and within 1kb of the next probe accounting for the interspersed density of background CpG sites within 400bp (Gatev et al., 2020). Beneath those is the DYNC1I1 gene association based on the UCSC genome track (Kent et al., 2002), which contains the entire CMR, followed by the repressed chromatin state of the region determined by ChromeHMM histone marks (Ernst & Kellis, 2017). Lastly displayed are boxplots of DNAm β value, on a scale from zero to one, at each CpG in the CMR. Lines indicate average values and colors refer to attachment style (Type A, B, C, or D). There were differences greater than ǀ∆β ≥ 0.03ǀ across attachment style at the three medium-confidence CpGs.

Referencing BECon, 85% of hits (22/26) were variable in both blood and brain by |≥0.03| (Table 3; Supplemental Figure 4) (Edgar, Jones, Robinson, et al., 2017). Additionally, using the genotyping data, no mQTLs within 50kb were found to associate with any medium-confidence CpGs. A sensitivity analysis examined the ∆βs across each attachment type comparison for the detected medium-confidence CpGs (Supplemental Table 4). The largest number of comparisons within the 26 identified medium-confidence CpGs that passed the technical effect size threshold of ǀ∆β ≥ 0.03ǀ were between Type B and Type C, while the least were between Type A and Type C.

3.3 Thirteen of the 26 medium-confidence CpGs were within CMRs

Several of the identified CpGs were proximally located to CpGs with correlated DNAm patterns. Four CMRs contained multiple hits within the medium-confidence CpGs, containing 38% of the identified CpGs (Table 3). The thirteen CpG LINC01556 CMR was within an active promoter region, the three CpG PPFIA3 CMR was within a heterochromatin region, the four CpG intergenic regulatory CMR was within a strong enhancer (Figure 4), while the four CpG DYNC1I1 CMR was within a repressed chromatin state region (Figure 5). DCHS2, C17orf98, and CATSPER4 CMRs each contained one medium-confidence CpG, respectively (Table 3).

3.4 Immune system action was enriched in the infant DNAm associations with toddler type B classification

To better understand the biological pathways represented by DNAm associations with Type B, a GO analysis via logistic regression model (Mi et al., 2012) was conducted on the results of the secondary EWAS model that accounted for all cell types (Supplemental Table 3). Using the KEGG pathway reference, the top seven pathways represented immune system or cell activity and signaling (Table 4). While not immediately apparent in its relevancy, the eighth pathway, chagas disease, is characterized by increased monocyte proportions (Gómez-Olarte et al., 2019). Similarly, using additional references (Reactome and GO), the top pathways were again immune-related with B cell activation, interleukin signaling, and T-cell and leukocyte activity (Supplemental Tables 5 and 6).

3.5 Caregiver score for parent–child interaction quality contributed to the association between infant medium-confidence CpGs and future type B development

None of the birth characteristic factors contributed to the discovered CpG associations with Type B, nor did maternal postnatal depression or the PCITS child scales for parent–child interaction quality (Supplemental Table 7). However, the PCITS caregiver cognitive growth fostering subscale contributed to approximately 25% of the effect at 2 CpGs: cg17889682 within the DYNC1I1 CMR (F(7, 79) = 2.727, FDR = 0.0560, R² = 0.19, 21.05% contribution) and cg14640466 (F(7, 79) = 28.11, FDR = 0.0392, R² = 0.71, 23.99% contribution). Both CpGs are variable in brain and are concordant .20 and .52, respectively, with DNAm in the prefrontal cortex (Supplemental Figure 4). No other subscale contributed to the medium-confidence CpGs (Supplemental Table 7).

4. Discussion

These findings demonstrated the association between increased estimated infant monocyte proportion and DNAm at 3-months with Type D attachment style categorization at 22-months, though there was no difference in estimated infant monocyte proportions between future secure and non-secure categorization. Additionally, the majority of secure attachment (Type B) development associated CpGs were variable in, and correlated with, brain DNAm (Edgar, Jones, Meaney, et al., 2017) and within CMRs (Gatev et al., 2020). Though birth characteristics and maternal postnatal depression did not contribute to these associations, maternal cognitive growth fostering contributed a quarter of the effect to two CpGs. While no mQTLs were found to associate with the medium-confidence CpGs, this does not discount previous findings of a high genetic contribution to overall DNAm variation explained by attachment (Garg et al., 2018). The posthoc analysis of ∆βs across specific attachment type comparisons revealed, understandably, that the reported CpGs are a signal of future secure attachment style and mostly not informative of biological differences between insecure styles. Together, these results suggest the presence of a precursor biomarker of later attachment style development among individuals as young as 3-months-old.

4.1 Many medium-confidence DNAm associations with toddler secure attachment classification were within genes related to cognitive function and neurobiology

Positive early social environments, like those associated with a secure attachment style, can be a crucial resource in cognitive development and buffering later cognitive decline (Brummelte, 2016; Meins, 1997; Walsh et al., 2019). This may be related to both the vital importance of infants’ social environment and the biological plasticity of the early developmental period (Jones et al., 2018; Kolb et al., 2012). Evidence shows securely attached children develop improved cognitive skills in comparison to insecure children (Bergman et al., 2010; Edelstein, 2006; Meins, 1997; Moss & St-Laurent, 2001; De Ruiter & van Ijzendoorn, 1993; Sroufe, 2005; West et al., 2013). For example, secure attachment style at 15-months-old associated with better performance on both teacher ratings and all tasks of executive function at 5-years-old (Bernier et al., 2015).

Perhaps, these benefits of secure attachment to cognitive function are reflected in the contribution of caregiver cognitive growth fostering – the only detected contributor to medium-confidence secure attachment style CpG effects – as maternal cognitive growth fostering previously associated with higher cognitive function scores in 9-months-olds (Fuller et al., 2010). Creating an early environment where safe exploration and critical thinking are encouraged by caregivers may result in detectable biological differences associated with cognition. This contribution may suggest the parental initiative towards nurturing cognitive growth is important to the process of future secure attachment development and its associated biomarkers, which should be further investigated in prospective work. Both CpGs contributed to by cognitive growth fostering were variable in brain and correlated with Broadman area 10 (BA10) DNAm. One of these CpGs was within DYCN1I1, which, when expressed in brain, regulates mRNA localization along axons and synapses (Gershoni-Emek et al., 2016). Similarly, axon guidance was significant with secure attachment style in the KEGG pathway, replicating previous findings of DNAm and attachment (Garg et al., 2018) and echoing potential differences in cognition and plasticity.

Several medium-confidence CpGs have been associated with cognitive outcomes or reside within genes that do. DNAm in HK2, a gene highly expressed by neurons (Lundgaard et al., 2015), has associated with cognitive processing (Marioni et al., 2018). PPFIA3, a gene important for neurotransmitter release across the synapse, and four intergenic medium-confidence CpGs associated with DNAm differences in Trisomy-21, potentially related to its neurodiverse cognitive development phenotype (Haertle et al., 2019; Henneman et al., 2018). GSG1L expression negatively regulates the strength of AMPA receptor-mediated synaptic transmission, a critical mechanism for synaptic plasticity, and plays a pivotal role in information processing (Aichurok Kamalova et al., 2020; Gu et al., 2016; Jacobi & Engelhardt, 2021; Mao et al., 2017). Finally, COMT is well-characterized for its role in catecholamine neurotransmitter release and regulation. Both SNP variation and DNAm within COMT have been associated with cognitive outcomes, especially in association with early life experience (Bruder et al., 2005; Malhotra et al., 2002; Ursini et al., 2011; Vidrascu et al., 2019; Walton et al., 2014). Though this study only examined blood DNAm, the high concordance and variability in brain at the majority of CpGs, especially within the prefrontal cortex (BA10), (Table 3; Supplemental Figure 4) and established relation of attachment and cognition, lend credence to these CpGs as potential precursor biomarkers of cognitive and synaptic advantages associated with future secure attachment development.

4.2 Infant proinflammatory immune biomarkers differed with toddler secure and disorganized attachment styles

Both the increased estimated monocyte cell type proportion in the blood of disorganized attachment classification (Type D) infants and the medium-confidence CpG associations with secure attachment classification (Type B), indicate the possibility of proinflammatory immune differences during early attachment development. Monocytes are white blood cells that are important in both the adaptive and innate immune systems and normally increase in response to an immune challenge (Prame Kumar et al., 2018). Estimated higher basal proportions of monocytes in blood, like Type D infants in this study, may suggest the presence of chronic inflammation in response (Jong et al., 2017). Thus, it is reasonable to suppose future Type D infants may exhibit an increased inflammatory profile at 3-months-old reflected in identified DNAm differences.

These results echo several existing attachment findings. There is substantial evidence for dysregulated hypothalamic-pituitary-adrenal (HPA) axis and immune responses with attachment insecurity (Pietromonaco & Powers, 2015). Children categorized as Type D at 17-months-old had higher inflammation than organized attachment styles (Measelle et al., 2017). Children who experienced insecure attachment at 1-year-old had increased risk for inflammation-based illness 31 years later (Puig et al., 2013). Moreover, children with Type D attachments during infancy had higher levels of C-reactive protein (CRP), a marker of inflammation induced by white blood cell cytokines (Samson et al., 2019), than Type B children at 5-years-old (Bernard et al., 2019).

While not directly tested in this study, the literature implies a potential mutual origin of association for attachment insecurity and early life adversity through chronic stress. For example, higher estimated monocyte proportion in infancy for Type D infants, who are in the most stressful social environment of the four attachment types (Cyr et al., 2010), aligns with the prediction made by the biological embedding model of early adversity (Miller et al., 2011). This model posits exposure to early life stress may relate to aggressive monocyte cells nonresponsive to inhibition, thus leading to a proinflammatory phenotype and increased health risks (Ehrlich et al., 2016). One study directly testing this theory in adolescent girls extracted and stimulated monocytes, measuring IL-6 production. They found early life adversity related to increased IL-6 and a proinflammatory phenotype (Ehrlich et al., 2016). In our study, DNAm at CpGs within QRICH2 have also been previously associated with adverse childhood experiences in young adult blood (Artigas et al., 2020) and may be associated with immune cell changes (Maltby et al., 2018).

Moreover, two candidate CpGs originally identified via an association with parenting and inflammation within TNF, a gene that encodes for cytokine tumor necrosis factor alpha (TNF-α) (Louis et al., 1998), statistically associated with secure attachment. There were also three medium-confidence CpGs within DYNC1I1, a gene in the immune super-pathway which encodes a protein that upregulates IL-6 through nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) in vitro (Gong et al., 2019). The DNAm differences discovered in DYNC1I1 CpGs was the exact opposite of monocyte cell type proportions, with Type D DNAm the lowest and Type B DNAm the highest. Because the chromatin state of this region is repressed, it is reasonable to speculate these CpGs are normatively more methylated and DNAm could inversely relate to expression, though we did not measure gene expression in this study (Mendenhall & Bernstein, 2008). Additionally, the GO analysis resulted in a profusion of immune-related pathways, with the biological action of leukocytes (e.g. monocytes) as prevalent across all references, again granting more credence to proinflammatory biological differences during secure attachment development.

Stress experienced early in life may affect adult immune health through this biological embedding (Boyce & Kobor, 2015; Hertzman, 1999, 2012). Early-life adversity, including poverty, psychosocial stress, parental psychopathology, and negative caregiver experiences predict inflammation markers in late childhood (Miller & Chen, 2010; Slopen et al., 2013) and adulthood (Austin et al., 2018; Luecken & Lemery, 2004; Repetti et al., 2002). For example, one study demonstrated that a composite measure of proinflammatory cytokines (e.g. IL-6, TNF-α) and CRP were significantly associated with different forms of early adversity (Hoye et al., 2020).

However, there is substantial evidence that secure attachment style provides, and is a core component of, resilience to adversity (Darling Rasmussen et al., 2019). For instance, early maternal warmth can buffer the effect of early adversity proinflammatory signaling in adulthood (Chen et al., 2011). Also, increased proinflammatory cytokines and CRP associated with maternal depression and were moderated by Type B in 17-month-old infants (Measelle & Ablow, 2018). These results suggest the possibility that biomarkers of proinflammatory immune differences in response to both positive, resulting in secure attachment, and adverse, resulting in insecure or disorganized attachment, early life social environments can be observed as early as 3-months-old and may have implications for future attachment and adult health outcomes (Fraley, 2019). These data support an association between toddler secure attachment development and infant DNAm that suggest differences in early life proinflammatory immune action and cognitive development. Our findings illustrate both the potential consequences of early adversity and resilience of positive parenting through biological associations within the first 3-months of life.

4.3 Limitations

We acknowledge these analyses lack statistical power due to a small sample, increasing the potential for false negatives. This may have contributed to the lack of associations found among candidate CpGs and potential mQTLs. As a trade-off, we applied a moderate threshold when multiple-test correcting throughout all analyses (FDR<0.10) along with the more stringent technical effect size threshold of |Δβ≥0.03|. Thus, these associations need to be replicated and inferences made cautiously, especially as no gene expression data was reported. Similarly, venous blood is comprised of multiple cell types (Jones et al., 2018), including subpopulations not measured in references, and the associations of DNAm within specific cell types remain unknown. Our sample is homogeneous, primarily Caucasian, and higher socioeconomic status, reducing the ability to extrapolate these findings to high-risk or unrepresented groups. We also note there is a temporal disjoint between the DNAm measurement (3-months) and attachment style measurement (22-months), though this chronological difference allowed us to hypothesize about potential biological associations with the process of future attachment style development. Finally, though these findings provide support for the association of DNAm at 3-months-old with attachment style development, these data do not support causal mechanisms; however, we hope this report does provide the foundation for future studies into the biomarkers and processes of secure attachment formation.

Disclosure of potential conflicts of interest

The authors of this study have no affiliations with nor involvement in any organization or entity with any financial interests or non-financial interests in the subject matter or materials discussed in this manuscript.

9.1.1 Genotyping Microarray and Processing

A subset of the infants with DNAm profiling (n = 106) also had single nucleotide polymorphism (SNP) genotyping performed in addition. Extracted DNA was run on the Illumina Global Screening Array (GSA) measuring 654,027 SNPs. Quality control of SNPs was performed in plink2 (Chang et al., 2015) using genome build GRCh37/hg19. SNPs were filtered by a minor allele frequency of at least 1% and for missing genotype rates greater than 10%. Next, SNPs with <97% allelic invariability were removed, leaving a total of 106 individuals and 302,132 SNPs.

9.1.2 DNA Methylation Microarray and Processing

Genomic DNA was bisulfite converted using EZ-96 DNA Methylation kits (Zymo Research, Irvine, CA) and run on Infinium HumanMethylation450 (450K) BeadChips (Illumina) producing 485,577 data points for 482,421 Cytosine-phosphate-Guanosine (CpG) sites. Raw IDAT files were imported to RStudio where beta (β) values (signifying the amount of DNA molecules methylated at a site on a range from 0-1) were produced from intensity measures (Du et al., 2018).

Functional normalization was used for background subtraction, color correction and probe type normalization (Aryee et al., 2014; Fortin et al., 2017). Samples were removed as outliers if detected in any of three methods. The detectOutlier function (Du et al., 2018) checks if the sample mean is too far from the cohort mean. The outlyx function (Pidsley et al., 2013) checks if a sample does not cluster with the cohort using principal component analysis (PCA). The locfdr package (Efron, 2007; Efron et al., 2001; Hannum et al., 2013) checks if the z-score PCA of a sample is different from the cohort at a threshold FDR ≤ 0.2. Additionally, any sex mismatches were identified by comparing the reported sex with DNAm predicted sex (Table 1). We then removed probes on the X and Y chromosomes, with fewer than three beads contributing to signal, with NAs in more than 2% of samples, with poor detection p values (<1×10⁻¹⁶), or previously shown to have unreliable probe design (Price et al., 2013) (Table 1). Variation associated with batch was assessed using PCA and then removed using surrogate variable analysis (ComBat) (Johnson et al., 2007).

An inter-quantile range filter subset to probes where the DNAm β value varies by at least 5% across samples in the 5th and 95th percentile (Edgar, Jones, Robinson, et al., 2017). This reduced the testing space by removing probes that did not vary across individuals in our sample to better meet the Benjamini-Hochberg False Discovery Rate (FDR) control method assumptions of equal significance likelihood (Benjamini & Hochberg, 1995; Korthauer et al., 2019). In total, 114 samples and 96,339 probes were deemed variable and of good quality. Of the 114 quality blood samples, 93 individuals also had SSP data at 22-months and were included in the analyses (Table 1), resulting in an analysis with approximately 40% power as determined by the pwrEWAS package (Graw et al., 2019).

9.1.3 Cell Type Proportion Estimation

A major function of DNAm is in cell type differentiation (Aristizabal et al., 2019). Since buffy coat blood samples represent a mix of cell types we computationally derived six estimated cell type proportions using the Identifying Optimal Libraries (IDOL) adult reference set with the Houseman method (Houseman et al., 2012; Koestler et al., 2016): neutrophils, natural killer cells, monocytes, CD4 T cells, CD8 T cells, and B cells (Supplemental Figure 1). Because these samples are from infants, we also estimated cell types using umbilical cord blood and mixed IDOL adult and cord blood references (Gervin et al., 2019) (Supplemental Figure 1), though the adult reference is more appropriate due to the lack of nucleated red blood cells (nRBCs) in 3-month infant venous blood (Hermansen, 2001).

9.1.4 Estimated Cell Type Proportion Principal Component Analysis (PCA)

We used principal components to account for estimated cell type variability. For the primary model excluding any cell types with significant attachment style associations (monocytes), a classic PCA was performed using the prcomp function in the stats package on the remaining five cell types (Holland, 2019). Three principal components (PCs), representing 93% of that cell type variability were included as covariates. For the secondary model including all cell types, these data were compositional and add to a constant: one. Therefore, we used robust isometric logratio (ilr) PCA, which is appropriate for composite data (Filzmoser et al., 2009). Four PCs, representing 99% of the estimated cell type proportion variability, were included in the second model.

Acknowledgments

We are extremely grateful to all the families who took part in this study and the whole APrON team (http://www.apronstudy.ca/), investigators, research assistants, graduate and undergraduate students, volunteers, clerical staff and mangers. This cohort was established by an interdisciplinary team grant from Alberta Innovates Health Solutions (formally the Alberta Heritage Foundation for Medical Research) and additional funding from the Alberta Children’s Hospital Foundation and Kids Brain Health Network (formerly National Centre of Excellence NeuroDevNet) assisted with the collection and analysis of data presented in this manuscript.

Data Availability Statement

Raw data were generated at the University of British Columbia. Derived data and analysis scripts supporting the findings of this study are available from the corresponding author, SMM, on request.

9. Appendix

9.1 Biological Sample and Data Preprocessing Details

Venipuncture blood draws were performed by a certified pediatric phlebotomist at the Behavioural Research Unit in Alberta Children’s Hospital as dictated in the University of Calgary approved REB. Whole blood samples were collected and briefly stored at −80°C before being centrifuged at 3000 RCF for 15 minutes to separate and pipette. Samples were then shipped on dry ice to Dr. Kobor’s Lab at the University of British Columbia for DNA extraction and measurement.