DNA methylation in relation to gestational age and brain dysmaturation in preterm infants.

Preterm birth is associated with dysconnectivity of structural brain networks and is a leading cause of neurocognitive impairment in childhood. Variation in DNA methylation is associated with early exposure to extrauterine life but there has been little research exploring its relationship with brain development. Using genome-wide DNA methylation data from the saliva of 258 neonates, we investigated the impact of gestational age on the methylome and performed functional analysis to identify enriched gene sets from probes that contributed to differentially methylated probes or regions. We tested the hypothesis that variation in DNA methylation could underpin the association between low gestational age at birth and atypical brain development by linking differentially methylated probes with measures of white matter connectivity derived from diffusion MRI metrics: peak width skeletonized mean diffusivity, peak width skeletonized fractional anisotropy and peak width skeletonized neurite density index. Gestational age at birth was associated with widespread differential methylation at term equivalent age, with genome-wide significant associations observed for 8870 CpG probes (P < 3.6 × 10-8) and 1767 differentially methylated regions. Functional analysis identified 14 enriched gene ontology terms pertaining to cell-cell contacts and cell-extracellular matrix contacts. Principal component analysis of probes with genome-wide significance revealed a first principal component that explained 23.5% of the variance in DNA methylation, and this was negatively associated with gestational age at birth. The first principal component was associated with peak width of skeletonized mean diffusivity (β = 0.349, P = 8.37 × 10-10) and peak width skeletonized neurite density index (β = 0.364, P = 4.15 × 10-5), but not with peak width skeletonized fraction anisotropy (β = -0.035, P = 0.510); these relationships mirrored the imaging metrics' associations with gestational age at birth. Low gestational age at birth has a profound and widely distributed effect on the neonatal saliva methylome that is apparent at term equivalent age. Enriched gene ontology terms related to cell-cell contacts reveal pathways that could mediate the effect of early life environmental exposures on development. Finally, associations between differential DNA methylation and image markers of white matter tract microstructure suggest that variation in DNA methylation may provide a link between preterm birth and the dysconnectivity of developing brain networks that characterizes atypical brain development in preterm infants.

Introduction

Preterm birth, defined as birth at <37 weeks of gestation, affects around 11% of births worldwide[1] and is a leading cause of neurodevelopmental and cognitive problems that extend across the life course. These include autism spectrum disorder, social difficulties, language impairment, attention-deficit hyperactivity disorder, reduced intelligence quotient, educational underachievement and psychiatric diagnoses.[2-9]

The neural phenotypes that underlie long-term functional impairment include diffuse white matter injury and subsequent dysmaturation of white matter and grey matter neuroaxonal structures, collectively termed the ‘encephalopathy of prematurity’.[10] A consequence of the encephalopathy is generalized dysconnectivity of developing structural networks, which can be inferred from diffusion MRI (dMRI) and neurite orientation dispersion and density imaging (NODDI) during the neonatal period.[11-15] Specifically, normal maturation is characterized by a reduction in mean diffusivity (MD) and increases in both fractional anisotropy (FA) and neurite density index (NDI) in white matter; but MD is increased, and FA and NDI are decreased in preterm infants at term equivalent age, compared with control infants born at term.[16] These changes reflect an increase in water content and a decrease in white matter organization in preterm infants. The peak width of skeletonized mean diffusivity (PSMD) is a method for histogram-based calculation of MD distribution across the entire white matter skeleton; it provides a measure of generalized white matter microstructure, is robust to scanner variation and is predictive of cognition in later life.[17-19] In previous work, we extended the histogram model to neonatal data and included other dMRI and NODDI metrics. We found that PSMD and peak width skeletonized neurite density index (PSNDI) are altered in preterm infants at term equivalent age and that histogram-based measures have utility for investigating upstream determinants of dysmaturation such as systemic inflammation.[20,21]

The mechanisms that link the environmental stress of preterm birth with atypical brain development are uncertain. Variation in DNA methylation (DNAm) is a possible mechanism; DNAm is involved in the regulation of gene expression and cell fate during foetal brain development.[22] Alterations in DNAm contribute to the pathogenesis of neurodevelopmental disorders [Rett syndrome,[23] Immunodeficiency, Centromeric region instability, Facial anomalies (ICE) syndrome[24] and Angleman and Prader–Willi syndromes].[25] There is growing evidence that differential DNAm can mediate the effect of environmental pressures on brain structure and function across the life course.[26,27] The neonatal methylome is sensitive to prenatal factors such as maternal smoking,[28] maternal body mass index,[29] as well as birth weight.[30] It is altered in association with co-morbidities of preterm birth,[26,31,32] and there is some evidence for legacy differences in the methylome two decades after preterm birth.[33] A meta-analysis investigating DNAm from umbilical cord blood identified widespread differential methylation associated with GA at birth (across the range 27−42 weeks) affecting 2% of sites and involving both hypo- and hypermethylation, as measured on the Illumina 450k array.[34] However, due to differences in the cellular composition of samples, epigenetic signatures observed in different tissues are likely to be distinct.[35] The main cellular component of saliva, buccal epithelium, may be more representative of the brain than umbilical cord blood because of ectodermal origin.[36-38]

Here, our first aim was to determine whether low gestational age at birth was associated with variation in the salivary methylome at term equivalent age and to characterize the biological pathways implicated in the DNAm response to preterm birth. Our second aim was to investigate whether the DNAm signal of gestational age explains variance in measures of white matter microstructure at term equivalent age. We tested the hypotheses that low gestational age at birth is associated with widespread differential methylation apparent at the end of neonatal intensive care and that DNAm contributes to variance in peak width skeletonized metrics of white matter connectivity.

Materials and methods

Participants

All participants were born at the Royal Infirmary of Edinburgh, UK. Preterm infants, defined as GA of <33 weeks of gestation and term infants, defined as GA > 37 weeks based on first-trimester ultrasound scan dating, were recruited to the Theirworld Edinburgh Birth Cohort. This is a longitudinal study designed to investigate the effect of preterm birth on brain development.[39] Exclusion criteria were major congenital malformations, chromosomal abnormalities, congenital infection, overt parenchymal lesions (cystic periventricular leukomalacia, haemorrhagic parenchymal infarction) or posthaemorrhagic ventricular dilatation. Ethical approval has been obtained from the National Research Ethics Service, South East Scotland Research Ethics Committee (11/55/0061, 13/SS/0143 and 16/SS/0154). Informed consent was obtained from a person with parental responsibility for each participant. The study was conducted according to the principles of the Declaration of Helsinki. DNAm data were available from 258 neonates, 214 of whom also had successful dMRI acquisition.

DNA extraction and methylation measurement

Saliva obtained at term equivalent age, on the same day as MRI acquisition, was collected in Oragene OG-575 Assisted Collection kits, by DNA Genotek, and DNA extracted using prepIT.L2P reagent (DNA Genotek, ON, Canada). DNA was bisulphite converted and methylation levels were measured using Illumina HumanMethylationEPIC BeadChip (Illumina, San Diego, CA, USA) at the Edinburgh Clinical Research Facility (Edinburgh, UK). The arrays were imaged on the Illumina iScan or HiScan platform, and genotypes were called automatically using GenomeStudio Analysis software version 2011.1 (Illumina). DNAm was processed in two batches.

DNA methylation preprocessing

Raw intensity (.idat) files were read into the R environment (version 3.4.4) using minfi. wateRmelon and minfi were used for preprocessing, quality control and normalization.[40,41] The pfilter function in wateRmelon was used to exclude samples with 1% of sites with a detection P-value > 0.05; sites with beadcount <3 in 5% of samples and sites with 1% of samples with detection P-value > 0.05. Cross hybridizing probes and probes targeting single nucleotide polymorphisms with overall minor allele frequency ≥0.05 were also removed.[42] Control probes were also removed. Samples were removed if there was a mismatch between predicted sex (minfi) and recorded sex (n = 3). Data were danet normalized which includes background correction and dye bias correction.[41] Saliva contains different cells types, including buccal epithelial cells and leucocytes. Epithelial cell proportions were estimated with epigenetic dissection of intra-sample heterogeneity with the reduced partial correlation method implemented in the R package EpiDISH.[43] Probes located on sex chromosomes were removed before analysis. Data from one of each twin pair were removed randomly (n = 20).

MRI acquisition

MRI was obtained at the same appointment as saliva sample collection for DNAm analysis. Structural and dMRI were performed on 93 infants using a MAGNETOM Verio 3T clinical MRI scanner (Siemens Healthcare GmbH, Erlangen, Germany) and 12-channel phased-array head coil, which were used to acquire dMRI using a protocol consisting of 11 T2- and 64 diffusion-weighted (b = 750 s/mm2) single-shot spin-echo echo-planar imaging (EPI) volumes acquired with 2 mm isotropic voxels (echo time (TE) = 106 ms and repetition time (TR) = 7300 ms).

One hundred and twenty-one infants were scanned using a MAGNETOM Prisma 3T clinical MRI scanner (Siemens Healthcare GmbH) and 16-channel phased-array paediatric head and neck coil. dMRI was acquired in two separate acquisitions: the first consisted of eight baseline volumes [b = 0 s/mm2 (b0)] and 64 volumes with b = 750 s/mm2; the second consisted of 8 b0, 3 volumes with b = 200 s/mm2, 6 volumes with b = 500 s/mm2 and 64 volumes with b = 2500 s/mm2. An optimal angular coverage for the sampling scheme was applied.[44] In addition, an acquisition of three b0 volumes with an inverse phase encoding direction was performed. All dMRI volumes were acquired using single-shot spin-echo EPI with 2-fold simultaneous multislice and 2-fold in-plane parallel imaging acceleration and 2 mm isotropic voxels; all three diffusion acquisitions had the same parameters (TR/TE 3500/78.0 ms). Images affected by motion artefact were re-acquired multiple times as required; dMRI acquisitions were repeated if the signal loss was seen in three or more volumes.

Infants were fed and wrapped and allowed to sleep naturally in the scanner without sedation. Pulse oximetry, electrocardiography and temperature were monitored. Flexible earplugs and neonatal earmuffs (MiniMuffs, Natus) were used for acoustic protection. All scans were supervised by a doctor or nurse trained in neonatal resuscitation. Structural images were reported by an experienced paediatric radiologist (A.J.Q.)

dMRI preprocessing

Diffusion images that were acquired on the MAGNETOM Verio scanner were denoised using a Marchenko–Pastur-principal component analysis (PCA)-based algorithm[45-47]; eddy current distortion and head movement were corrected using outlier replacement;[48-50] bias field inhomogeneity correction was performed by calculating the bias field of the mean b0 volume and applying the correction to all the volumes.[51] FA and MD were calculated from the dMRI data.

The two dMRI acquisitions from the MAGNETOM Prisma scanner were first concatenated and then denoised using a Marchenko–Pastur-PCA-based algorithm[45-47]; eddy current, head movement and EPI geometric distortions were corrected using outlier replacement and slice-to-volume registration;[48-50,52] bias field inhomogeneity correction was performed by calculating the bias field of the mean b0 volume and applying the correction to all the other volumes.[51] From the dMRI data, we calculated the three eigenvalues and eigenvectors of the water diffusion tensor, and NODDI (Bingham distribution) parametric maps using cuDIMOT [intracellular volume fraction (NDI) and the overall orientation dispersion index (ODITOT)].[12,13,53] FA and MD were calculated using single-shell data to match the Verio scanner.

The peak width of skeletonized water diffusion parameters

All the subjects were registered to the Edinburgh Neonatal Atlas (ENA50) using DTI-TK.[20] The diffusion tensor derived maps of each subject (FA and MD) were calculated after registration; NDI was then propagated to the template space using the previously calculated transformations. The data were skeletonized using the ENA50 skeleton and then multiplied by a custom mask. Finally, the peak width of the histogram of values computed within the skeletonized maps was calculated as the difference between the 95th and 5th percentiles.[17,20]

Statistical analysis

Epigenome-wide association analyses

Unless otherwise stated, analysis was completed in R version 3.4.4.[54] An overview of the analysis pipeline is shown in Supplementary Fig. 1. Surrogate variable analysis (SVA) of the data matrix was carried out, to adjust for potential confounders such as batch, using the sva function in the sva package in R.[55,56] A fully adjusted model was specified before SVA to retain signal explained by biological variables of interest: CpG ∼ gestational age at birth + age at scan + birthweight Z-score + maternal smoking + sex + epithelial cells. SVA identified 17 significant surrogate variables (SVs) which were subsequently included in the analysis.

An epigenome-wide association study (EWAS) was performed using the limma package in R.[57] Beta values of each of 776 025 CpG sites were regressed (as dependent variables) on gestational age (GA) at birth using linear regression. Covariates were added to adjust for sex, birthweight Z-score, age at sample collection, maternal smoking, estimated epithelial cell proportions and 17 surrogate variables. A significance threshold of 3.6 × 10−8 was selected, which represents genome-wide significance.[58]

Differentially methylated region analysis

Differentially methylated regions (DMRs) were assessed using the dmrff function in the dmrff package in R.[59] Here, a differentially methylated region is a region containing two or more sites separated by ≤500 bp with EWAS analysis P ≤ 0.05 and methylation changes in a consistent direction. Following dmrff’s subregion selection step, DMRs with Bonferroni-adjusted P ≤ 0.05 were significant.

Gene set testing

Gene set enrichment analysis was carried out using the gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) databases, and using the gometh function in missMethyl package, which controls for multiple probe bias.[60] We performed an analysis that included those sites that reached genome-wide significance in EWAS and a second analysis that also incorporated those sites that contributed to differentially methylated regions.

Principal component analysis

Principal component analysis (PCA) was conducted on CpG probes that reached genome-wide significance, using the prcomp function in R. CpGs were precorrected for the effects of biological covariates and surrogate variables via linear regression. The scree plot was visually inspected to select a principal component (eigenvalue > 1) to be carried forward for subsequent analysis.

Linear regression between DNAm and peak width skeletonized metrics

Pearson’s correlation coefficient was used to assess the relationship between the first PC identified from PCA and gestational age at birth. This PC was used in linear regression models, as an independent variable, to test the associations between peak width of skeletonised mean diffusivity (PSMD), peak width of skeletonised neurite density index (PSNDI) and peak width of skeletonised fractional anisotropy (PSFA) with DNAm, conducted in R version 4.0.1.[54] In models testing PSMD and PSFA, MRI scanner was included as a binary covariate as MRI data from two scanners were included. PSNDI was only available from data acquired on one scanner, so it was not necessary to include the scanner as a covariate. We report standardized regression coefficients and P-values.

Data availability

The atlas with templates can be found at https://git.ecdf.ed.ac.uk/jbrl/ena and the code necessary to calculate histogram-based metrics is at https://git.ecdf.ed.ac.uk/jbrl/psmd. Requests for original image data will be considered through the BRAINS governance process: www.brainsimagebank.ac.uk.[61] DNAm data available upon request from Theirworld Edinburgh Birth Cohort, University of Edinburgh (james.boardman@ed.ac.uk or https://www.tebc.ed.ac.uk/2019/12/data-access-and-collaboration/). The DNAm and metadata are not publicly available due to them containing information that could compromise participant consent.

Results

Cohort

DNAm data were collected from 311 neonates. Twenty-nine did not meet DNAm preprocessing QC criteria and were excluded. One participant with a congenital abnormality was removed, as were three participants whose sex predicted from DNAm data did not match their recorded sex. This group of 278 neonates included 20 sets of twins. After random removal of one twin from each set, there was no evidence of imbalance for birthweight between the twins that were removed and the twins that remained in the sample (t = −0.157, P = 0.88) or for sex (χ2 = 0.417, P = 0.52).

The study group consisted of 258 neonates: 155 participants were preterm (gestational age range 23.28–34.84) and 103 were infants born at full term (gestational age range 36.42–42.14), see Table 1 for participant characteristics and Supplementary Fig. 2 for participant flow. Among the preterm infants, 38 (25%) had bronchopulmonary dysplasia (defined as need for supplementary oxygen ≥36 weeks gestational age), 9(6%) developed necrotizing enterocolitis requiring medical or surgical treatment and 31 (20%) had an episode of postnatal sepsis defined as either blood culture positivity with a pathogenic organism, or physician decision to treat for ≥5 days in the context of the growth of coagulase-negative staphylococcus from blood or a negative culture. Of the 258 participants with DNAm data, 214 also had MRI data.

Table 1

Participant characteristics

	Preterm infants (n = 155)	Term infants (n = 103)	P-value[a]
Gestational age at birth/weeks (range)	28.84 (23.28–34.84)	39.7 (36.42–42.14)	<0.05
Gestational age at scan/weeks (range)	40.56 (37.70–45.14)	42.27 (39.84–47.14)	<0.05
Birth weight/g (range)	1177 (500–2100)	3482 (2346–4670)	<0.05
Birth weight Z-score (range)	−0.19 (range −3.13–1.58)	0.43 (range −2.30–2.96)	<0.05
Sex: female (%)	75 (48)	44 (43)	0.37
Maternal folate supplementation in pregnancy (%)	136 (88)	86 (83)	0.33
Maternal age (years)	31.1 (17–44)	33.7 (19–48)	<0.05
Maternal tobacco smoker in pregnancy (%)	29 (19)	2 (2)	<0.05
Maternal diabetes (%)	10 (6)	6 (6)	0.84
Pregnancy-induced hypertension (%)	22 (14)	7 (7)	0.07
Highest maternal qualification			<0.05
None	4 (3)	0 (0)
High school	43 (28)	7 (7)
College/university	106 (68)	95 (92)

Student’s t-test was used to analyse continuous variables and χ2 to analyse proportions. There were three missing for maternal qualifications (two preterm, one term).

Widespread differential saliva DNAm in association with gestational age at birth

We conducted an epigenome-wide association study whereby CpG methylation at 776 025 sites was regressed on gestational age at birth, controlling for birthweight Z-score, infant sex, age at sample collection, maternal smoking, estimated epithelial cell proportion and surrogate variables. The genomic inflation factor was 1.72 (Supplementary Fig. 3). Differential methylation in relation to gestational age at birth was identified at 8870 CpG sites at genome-wide significance (P < 3.6 × 10−8), Fig. 1. Of these, 4250 (47.9%) sites demonstrated a positive association with gestational age, while 4620 (52.1%) were negatively associated. Following Bonferroni adjustment, 1767 DMRs corresponding to 4664 CpG sites were significant at P < 0.05. Of these, 11 had 10 or more CpG sites contributing to the DMR. The largest DMR mapped to a genomic region that encodes two genes: NNAT and BLCAP. The 29 probes mapped to this region were all located within islands and positively associated with gestational age at birth, indicating that longer gestation corresponds to hypermethylation. Of the 10 most significant differentially methylated probes (DMPs), three probes were localized to the IRX4 gene, one probe to the GAL3ST4 gene, and one to LOXL4 (Table 2; Supplementary Fig. 4). The probes with the largest absolute magnitude effect size (top five hypermethylated and hypomethylated) were mapped to the following genes: IRX2, SMIM2, INTS1, HEATR2, ZBP1 and UBXN11 (Table 3).

Figure 1

Manhattan plot for the significance [−log10 (. The solid horizontal line shows the genome-wide significance level and red dots above this line represent probes that are significant at this threshold (P < 3.6 × 10−8).

Table 2

Most significant probes associated with gestational age at birth

Probe	Chromosome	P-value	Gene	Direction of effect	Coefficient[a]	Standard error	Relation to island
cg03558436	5	1.26 × 10−44	−	+	1.02 × 10−2	5.82 × 10−4	Open Sea
cg04466438	9	1.13 × 10−42	−	+	7.55 × 10−3	4.47 × 10−4	Open Sea
cg23701943	10	1.11 × 10−41	LOXL4	+	1.04 × 10−2	6.29 × 10−4	Open Sea
cg18172877	5	2.31 × 10−39	IRX4	−	−6.12 × 10−3	3.85 × 10−4	Island
cg04180086	5	3.22 × 10−39	IRX4	−	−7.36 × 10−3	4.64 × 10−4	Island
cg22645539	7	1.22 × 10−38	GAL3ST4	−	−8.15 × 10−3	5.20 × 10−4	North Shelf
cg17774559	5	2.09 × 10−38	IRX4	−	−6.35 × 10−3	4.06 × 10−4	Island
cg17582074	4	4.12 × 10−38	−	+	5.25 × 10−3	3.38 × 10−4	Open Sea
cg08915267	13	6.01 × 10−38	−	−	−5.15 × 10−3	3.32 × 10−4	North Shelf
cg04441405	5	6.60 × 10−38	−	−	−1.10 × 10−2	7.13 × 10−4	Island

Coefficient corresponds to methylation change per week of gestation.

Table 3

Probes with the largest absolute magnitude of association with gestational age at birth

Probe	Chromosome	P-value	Gene	Direction of effect	Coefficient[a]	Standard error	Relation to island
cg10402321	1	3.11 × 10−36	UBXN11	−	−1.14 × 10-2	7.60 × 10−4	Open Sea
cg04441405	5	6.60 × 10−38	−	−	−1.10 × 10−2	7.13 × 10−4	Island
cg07167946	5	1.94 × 10−32	IRX4	−	−9.85 × 10−3	7.13 × 10−4	Island
cg07803375	7	3.6 × 10−22	HEATR2	−	−9.08 × 10−3	8.50 × 10−4	North Shelf
cg14670058	13	9.24 × 10−23	SMIM2	−	−9.07 × 10−3	8.35 × 10−4	Open Sea
cg16051275	6	7.53 × 10−36	−	+	1.23 × 10−2	8.30 × 10−4	Open Sea
cg11460314	20	4.21 × 10−20	ZBP1	+	1.24 × 10−2	1.24 × 10−3	Open Sea
cg04118102	17	1.00 × 10−30	−	+	1.31 × 10−2	9.86 × 10−4	South Shelf
cg17368297	16	1.55 × 10−25	−	+	1.40 × 10−2	1.20 × 10−3	Open Sea
cg14576951	7	6.73 × 10−30	INTS1	+	1.44 × 10−2	1.10 × 10−3	Island

Coefficient corresponds to methylation change per week of gestation.

Pathways implicated in functional testing

Based on the 8870 sites that met the genome-wide significance threshold (P < 3.6 × 10−8), no KEGG terms remained significant following false discovery rate (FDR) correction for multiple comparisons. Two GO terms were enriched following FDR correction: one for anchoring (GO:0070161; q = 0.0062) and one for adherens junction (GO:0005912; q = 0.0062). In an analysis that incorporated all 11 752 distinct CpGs from both EWAS and DMR analysis, 14 GO terms were enriched (Table 4).

Table 4

Gene ontology terms that were significantly enriched in an analysis of all probes that contributed to DMPs and DMRs

Gene ontology	Term	FDR Q value	Number of probes associated with the gene ontology/total number of probes in the ontology	Type	Description
GO:0005912	Adherens junction	0.00002317	217/545	Cellular component	A cell–cell junction composed of the epithelial cadherin–catenin complex.
GO:0005925	Focal adhesion	0.00907565	154/404	Cellular component	A cell–substrate junction that anchors the cell to the extracellular matrix and that forms a point of termination of actin filaments.
GO:0007155	Cell adhesion	0.01728820	428/1413	Biological process	The attachment of a cell, to another cell or to the extracellular matrix, via cell adhesion molecules.
GO:0007167	Enzyme-linked receptor protein signalling pathway	0.01243116	322/1043	Biological process	A series of molecular signals initiated by the binding of an extracellular ligand to a receptor on the target cell plasma membrane, where the receptor possesses catalytic activity or is closely associated with an enzyme such as a protein kinase, and ending with regulation of a downstream cellular process, e.g. transcription
GO:0007169	Transmembrane receptor protein tyrosine kinase signalling pathway	0.00907565	237/719	Biological process	A series of molecular signals, initiated by the binding of an extracellular ligand to a tyrosine kinase receptor on the target cell plasma membrane, ending with regulation of a downstream cellular process.
GO:0022610	Biological adhesion	0.01266030	431/1420	Biological process	The attachment of a cell to a substrate, another cell, including intracellular attachment between membrane regions.
GO:0030029	Actin filament-based process	0.03635020	257/756	Biological process	Any cellular process that depends upon, or alters, the actin cytoskeleton (comprising actin filaments and their associated proteins).
GO:0030036	Actin cytoskeleton organization	0.04472556	229/664	Biological process	The assembly, arrangement of constituent parts or disassembly of cytoskeletal structures comprising actin filaments and their associated proteins.
GO:0030054	Cell junction	0.00907565	422/1296	Cellular component	Forms a specialized region of connection between two or more cells, or between a cell and the extracellular matrix, or between two membrane-bound components of a cell, such as flagella.
GO:0030055	Cell–substrate junction	0.00907565	155/411	Cellular component	A cell junction between a cell and the extracellular matrix.
GO:0034330	Cell junction organization	0.04159681	119/290	Biological process	The assembly, arrangement of constituent parts, or disassembly of a cell junction. A cell junction is a specialized region of connection between two cells or between a cell and the extracellular matrix
GO:0045296	Cadherin binding	0.00207297	130/331	Molecular function	Interacting selectively and non-covalently with cadherin, a Type I membrane protein involved in cell adhesion.
GO:0050839	Cell adhesion molecule binding	0.00207297	186/499	Molecular function	Interacting selectively and non-covalently with a cell adhesion molecule.
GO:0070161	Anchoring junction	0.00002317	222/561	Cellular component	A cell junction that mechanically attaches a cell, and its cytoskeleton, to neighbouring cells or the extracellular matrix.

Terms in bold were enriched in an analysis of 8870 genome-wide significant DMPs.

Gestational age at birth is associated with metrics of white matter microstructure in neonates

DNAm and PSMD and PSFA were both available for 214 infants (Supplementary Table 1), and DNAm and PSNDI were available for the 121 infants from Phase 2 (Supplementary Table 2). Gestational age at birth was significantly associated with PSMD (β = −0.602, P < 2 × 10−16) and PSNDI (β = −0.594, P = 2.17 × 10−9) but not with PSFA (β = −0.005 P = 0.933) (Fig. 2A and B; Table 5).

Figure 2

Scatter plots with regression lines and 95% confidence intervals showing the relationships between gestational age at birth (weeks) and DNAm with PSMD and PSNDI, where peak width is the difference between the 95th and 5th centile of histogram values across the white matter skeleton. The associations between gestational age (weeks) and PSMD and PSNDI are shown in (A) and (B), respectively. The relationships between DNAm PC1 and PSMD and PSNDI are shown in (C) and (D), respectively. PS metrics are residualized for gestational age at scan; PSMD is additionally residualized for scanner variable.

Table 5

Associations between global metrics of white matter microstructure, DNAm first principal component (left) and gestational age (right)

PS metric	Metric ∼ PC1 DNAm + age at scan + scanner variable[a]		Metric ∼ gestational age at birth + age at scan + scanner variable[a]
	β	P	β	P
PSFA	−0.035	0.510	−0.005	0.933
PSMD	0.349	8.37 × 10−10	−0.602	<2 × 10−16
PSNDI	0.364	4.15 × 10−5	−0.594	2.17 × 10−9

Scanner variable was included in the model for PSFA and PSMD but not PSNDI because NODDI imaging was carried out for a subset using a single MRI scanner (n = 121). Bold print signifies significant associations.

Differential DNAm is associated with white matter microstructure

The first unrotated PC (PC1) derived from the 8870 genome-wide significant CpGs accounted for 23.5% of the variance, and the second PC accounted for 2.5% (Fig. 3A). There was no evidence of batch effects in the PCs (Supplementary Fig. 5). PC1 was significantly correlated with gestational age at birth (r = −0.622; P < 2.2 × 10−16) (Fig. 3B). PC1 was also positively associated with PSNDI (β = 0.364, P = 4.15 × 10−5), and in models adjusted for scanner it was positively associated with PSMD (β = 0.349, P = 8.37 × 10−10) but not PSFA (β = −0.035, P = 0.510) (Table 5). All models were adjusted for gestational age at scan.

Figure 3

Variation in DNAm probes selected by EWAS captured by principal component analysis, and the relationship with gestational age (weeks). A. A scree plot showing percentage of variance accounted for by the first 10 components, with a sharp elbow after the first PC, accounting for 23.5% of variance. B. A scatter plot, with regression line, showing the relationship between gestational age at birth (weeks) and PC1 (r = −0.622) with 95% confidence intervals.

Discussion

By studying a unique database of DNA linked to brain imaging in a population of preterm and term infants, we have identified extensive differential methylation in association with gestational age at birth and revealed an association between the principal axis of methylation and brain dysmaturation within the same sample. Differentially methylated regions and probes were distributed widely across the genome, indicating that gestation duration has a global effect on DNAm. Gene enrichment analysis of changes associated with gestational age identified gene sets pertaining to cell contacts and cytoskeleton. A single principal component that explained 23.5% of the variance in differential DNAm linked to preterm birth was closely associated with markers of generalized dysconnectivity across the white matter skeleton.

The data are consistent with studies that have reported associations between length of gestation and genome-wide variation in DNAm within foetal brain[22] and umbilical cord blood;[34,62] and widely distributed variation is reported in association with postmenstrual age at the sampling of preterm infants (a proxy for gestational age at birth).[31] The signature we identified in saliva sampled at term equivalent age included 233 probes that were previously shown to be differentially methylated in association with gestational age in a meta-analysis of umbilical cord blood samples that reported 8899 gestation-dependent CpGs.[34] The limited overlap could be explained by differences in the cellular composition of assessed tissues, different array types used to measure DNAm, or due to the time of sampling. Cord blood is collected at birth and so methylation changes at this time reflect foetal maturity and/or prenatal experience, whereas the methylation signature at term equivalent age reflects the allostatic load of early postnatal experiences, as well as the prenatal environment. We chose to sample at term equivalent age because postnatal co-morbidities of preterm birth and NICU care practices such as painful stress exposures alter DNAm profiles, and because cumulative DNAm variations over this time period may link exposure to behavioural outcome in preterm infants.[32,63,64]

Functional analyses of DMPs identified two enriched GO terms, for adherens and anchoring junctions. When distinct probes that contributed to both DMPs and DMRs were combined, GO analysis identified an additional 14 terms related to cell–cell adhesion, cell adhesions with the extracellular matrix and signalling from the extracellular membrane; 12 of these were also identified in the meta-analysis of gestational age effects on DNAm obtained at birth from umbilical cord blood.[34] The most significant DMR mapped to a site encoding two genes: NNAT and BLCAP. NNAT encodes the neural fate initiator neuronatin, the expression of which decreases throughout development;[65] there was a positive association with increasing gestational age at birth. Hypomethylation of NNAT is associated with a corresponding increase in the expression of neuronatin.[66]BLCAP encodes the bladder cancer-associated protein. This is a tumour suppressor that induces apoptosis, with high expression in brain and B lymphocyte. The candidacy of this locus, encoding NNAT and BLCAP, as a region of interest whose expression may be modified by perinatal exposures is supported by previous EWAS of gestational age, maternal body mass index, maternal smoking and schizophrenia.[22,29,62,67-70] In addition, BLCAP was found to have reduced methylation in placental samples from women who had preeclampsia.[71]

Probes that demonstrated the largest magnitude of effect in association with gestational age mapped to genes previously associated with gestational age or maternal risk factors in EWAS.[70] Hypermethylated probes were found in genes including ZBP1, which was identified in EWAS investigating gestational age and hypertensive disorders of pregnancy;[62,69,72]INTS1, which has been identified in EWAS of gestational age, hypertensive disorders of pregnancy, maternal body mass index, birthweight and breastfeeding duration.[22,30,62,67,72,73] Hypomethylated probes were found in genes including UBXN11, which was identified in studies of gestational age;[22,62] and IRX4. Three of the 10 most significant DMPs mapped to the IRX4 gene, all of which displayed a negative association with gestational age at birth. IRX4 is associated with cardiac development in vertebrates, including humans.[74] Its homologues have been implicated in retinal axon guidance in zebrafish, and neural patterning in Xenopus,[75,76] and it has been identified in previous EWAS of hypertensive disorders of pregnancy[72] and prenatal maternal stress.[77]

The novel pathways and genes implicated by EWAS studies of gestational age could provide a strong empirical basis for the selection of genes in targeted analyses in association with neuroimaging.[27] For example, one of the genes identified in our EWAS has been implicated in neurodevelopmental disorders; biallelic mutations in INTS1 have been associated with a rare neurodevelopmental syndrome characterized by intellectual disability.[78-80]

We used metrics of generalized white matter connectivity to assess relationships between DNAm and brain development because generalized dysconnectivity of structural networks is a hallmark of preterm brain dysmaturation.[10,14,81] PSNDI and PSMD were strongly associated with the first principal component of gestational age-dependent variation in DNAm but PSFA was not. This suggests that variations in DNAm could contribute to the higher variability in water content and intra-axonal volume that characterize preterm brain dysmaturation.[20] We have previously reported that differential DNAm is associated with FA of the genu of the corpus callosum and tract shape of the right corticospinal tract;[26] it is most likely that we did not observe an association between gestational age or DNAm with PSFA because this metric is subject to histogram shift,[20] meaning that although there are groupwise differences in FA values across the skeleton, the spread of values is the same. Associations between gestational age at birth and both DNAm and image markers of dysconnectivity, and between DNAm and image features, suggest that differential DNAm contributes, in part, to the relationship between gestational age and brain network dysconnectivity in preterm infants. However, some of the DNAm signatures may be related to underlying causes of prematurity that operate in foetal life, such as infection or preeclampsia. We cannot rule out that such changes are involved in the aetiology of preterm birth, which might preclude them from mediating specific associations with brain development.

The strengths of this study are that we studied a population of preterm and term infants across the gestational age range of 23–42 weeks, who were uniquely phenotyped with DNAm and dMRI. We sampled after the period of NICU care to capture the allostatic load of preterm birth. We measured DNAm from neonatal saliva samples, which has consistency with brain DNAm patterns and is non-invasive. The Illumina EPIC platform provided extensive coverage of the methylome (850 000 sites) and we controlled for cell composition. Finally, we used an image phenotype that is robust to scanner variation.[17] There are some limitations. First, control for cell composition was based on estimation of cell proportions rather than measurement, so we cannot rule out the possibility that some of the signal identified was related to variation in cell composition. Second, mediation analysis to assess causation was not possible because the association between the DNAm PC1 and peak width skeletonized metrics might result from the DNAm PC being derived from CpG sites that are associated with gestational age, so variance attributable to the mediating variable cannot be assumed. This could be addressed by out-of-sample validation, which will require other neonatal cohorts with both saliva methylome and dMRI data. Cohorts, with such data, that recruit neonates from across the gestational age continuum would also provide valuable replication of analyses described here. There was also some evidence of inflation of test statistics based on the genomic inflation factor. However, the genomic inflation factor is thought to provide an overestimate of inflation and corrections based on it may be overly conservative.[82] The value of the genomic inflation factor was also similar to those previously reported in neonatal cohorts.[34] In addition, visual inspection of our results via Manhattan plots suggests that our finding, of widespread differences in DNAm in relation to gestational age at birth, is in line with previous studies that have investigated this in cord blood and foetal brain tissue.[22,34] Finally, preterm birth itself can reflect maternal health or circumstance; in this cohort, mothers of preterm infants were more likely to have pregnancy-induced hypertension and the control mothers were more likely to have college or university qualifications. It is possible that these, or other maternal factors associated with preterm birth such as infection, chronic disease, poor nutrition and alcohol, tobacco or drug use, could influence the neonatal methylome at term equivalent age. A large prospective study with longitudinal sampling from birth would be required to test the hypothesis that maternal factors contribute to changes in the methylome apparent after neonatal intensive care.

In conclusion, gestational age at birth has a profound and widely distributed impact on the neonatal saliva methylome at term equivalent age, which reflects the allostatic load of preterm birth itself and postnatal exposures during neonatal intensive care. GO terms related to cell–cell contacts were enriched, indicating that cell contacts and organization are implicated in the phenotype. Associations between DNAm and PSMD and PSNDI suggest that variations in DNAm could contribute to white matter dysconnectivity commonly seen in preterm infants, and this analysis identified several genes and gene regions that could provide further insight into the molecular mechanisms by which early exposure to extrauterine life influences neurodevelopment.

Click here for additional data file.