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

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.

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.

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.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.