Study population

Between October 2009 and December 2010, samples from 170 highly educated pregnant women at 6 to 12 weeks of gestation were recruited within the Prenatal Early Life Stress (PELS) study in the University Hospitals Leuven, Belgium. A detailed description of data and sample collection can be found elsewhere [14].

Briefly, each trimester, the pregnant women completed several questionnaires regarding mental health and general well-being. From 80 of these 170 participating mothers, we were able to collect a cord blood sample at birth, of which DNA was subsequently isolated. Of these 80 samples, 45 DNA samples were selected for genome-wide DNA methylation analysis based on maternal psychological measurements, while all 80 DNA samples were subsequently used for the candidate gene confirmation study. A schematic overview of the study design is represented in Additional file 1: Figure S1.

Psychological measurements

At each trimester of pregnancy, future mothers were asked to complete questionnaires assessing depression, anxiety, and mother-fetus relationship. The revised pregnancy-related anxiety questionnaire (PRAQ) was used to assess worries and anxiety specifically regarding the participating mother’s pregnancy [16]. The validated revised version, which we used for the current study, contains three subscales: (1) fear of integrity of the fetus, (2) concerns about one’s own appearance, and (3) fear of the delivery itself [17]. We previously reported a high internal reliability for the total scale (Cronbach’s α ≥ 0.95) and good internal reliability for the subscales (Cronbach’s α ≥ 0.73) of the revised questionnaire [14]. To select the samples for EWAS, we focused on the fear of integrity subscale of PRAQ, which appeared as most associated with DNA methylation in our candidate gene studies [14, 15]. This subscale provides a score ranging from zero to seven and measures maternal anxiety regarding the health of the unborn infant and fear of integrity of the newborn.

Child HPA axis measurements

In order to measure infant HPA axis stress response as well as baseline activity, cortisol was measured in response to a vaccination at 2 and 4 months old and during the day at 12 months old. In both vaccination studies and the day profile study, saliva samples for cortisol measurements were taken using BD Visitec eye sponges (Becton, Dickson and Company, Waltham, USA). To determine cortisol levels in saliva a High Sensitivity Salivary Cortisol Enzyme Immunoassay Kit (ELISA kit, Salimetrics, Europe) was used. This assay captures the full range of salivary cortisol levels (0.003 to 3.0 μg/dL) while using only 25 μL of saliva per test and it was designed by Salimetrics to be resilient to the effects of interference caused by collection techniques that affect pH.

Cortisol was measured in response to a physical stressor at 2 and 4 months of age. All visits were scheduled in the afternoon to control for circadian rhythm. Mothers were asked not to feed their infants from 1 hour before the test until after the test since this might compromise assay performance by lowering sample pH and influencing bacterial growth. The routine vaccination consisted of two intramuscular injections: one injection for Pneumococcus and one combined injection for diphtheria, whooping cough, tetanus, poliomyelitis, and Haemophilus influenza type B. Saliva for cortisol measurements was collected at arrival to the examination room and at 15, 30, 45, and 60 min following the injection. The time between arrival and vaccination included the time when the parent undressed the infant, when both waited for the doctor to arrive, and when the doctor was taking the health history and the physical examination. The use of this vaccination stressor has a number of advantages: (1) it is routinely performed on all infants and therefore is convenient; (2) it is standardized, and so equal treatment to all infants can be assumed (especially if all the data are collected at the same facility); (3) the onset of the application of the aversive stimulation can be precisely determined; and (4) the inoculation is aversive enough to evoke behavioral distress [18]. The maximum value following vaccination was noted and the area under the curve with respect to increase (AUCi), initial response (value at 15 min minus baseline), and recovery slope (value at 60 min minus 15 min) were calculated. Data for several of these variables were skewed (Shapiro-Wilk p value < 0.05) and were logit transformed to correct skewness. If measurements for one-time point were missing, AUCi, initial response and recovery slope were not calculated for this individual.

When infants were 12 months of age, mothers were instructed to collect saliva samples of their infant at awakening, 30 min, 4, and 12 h after awakening. Samples were sent back to our lab in a prepaid, addressed envelope, and processed to measure cortisol levels. The cortisol awakening response (CAR) was calculated and a repeated t test was carried out to verify a significant cortisol increase between awakening and 30 min.

Illumina Infinium HumanMethylation450 BeadChip Assay

Cord blood DNA samples were processed by the Barts and the London Genome Centre (London, UK), where bisulfite treatment was carried out using the Zymo EZ DNA methylation kit (Zymo Research) prior to running the Infinium HumanMethylation450 BeadChip (HM450) arrays (Illumina Inc., California, USA). In total, 48 samples were randomized over four chips based on PRAQ score, maternal smoking, maternal alcohol use, gender of the newborn, parity and gravidity, gestational age at birth, and maternal age. One sample was replicated on each chip to account for batch effects [19], resulting in a sample size of 45 unique samples analyzed on the HM450 array.

Data processing and quality control

Raw data generated by the iScan Illumina array were imported using GenomeStudio software (Illumina, Inc.) and the subsequent quality control and normalization were implemented using the wateRmelon package in R (available from the Bioconductor repository http://www.bioconductor.org) [20]. Data clean-up included removal of samples and CpG probes with insufficient data quality. In more detail, samples were removed if more than 5% of its sites had a poor detection p value (> 0.01), and probes with a detection p value of more than 0.01 in more than 1% of samples or a bead count of less than 3 in more than 5% of samples were removed from the analysis. Furthermore, cross-hybridizing probes, probes containing a common single nucleotide polymorphism (SNP) in the sequence or within 10 bp of the sequence, and probes on the X and Y chromosomes were removed [21]. The final analysis included 414,733 probes, and all samples passed the stringent quality control. A schematic overview of the analysis from data preprocessing until validation can be found in the Supplementary information (Additional file 1: Figure S1).

Genome-wide DNA methylation analysis

Differentially methylated regions (DMRs) were first explored using the R package DMRcate [22]. The DMRcate statistical model was corrected for gender, gestational age of the newborn at birth, batch (chip Sentrix ID), and cell type compositions as calculated with the Houseman algorithm [23]. Using the top 500 uniquely annotated DMRs from DMRcate, we performed a pathway analysis. Enrichment for Gene Ontology (GO) classes was conducted using the over-representation analysis (ORA) method in ErmineJ with a minimum gene set of 5, and a maximum of 100 [24].

DMR verification and validation using Sequenom EpiTYPER

In order to verify the top DMRs from the DMRcate analysis, we used a second region-specific method: the Python module comb-p [25]. Comb-p employs probe locations and p values from the differentially methylated CpG probe (DMP) analysis as input to identify differentially methylated regions. We, therefore, first calculated DMPs using a linear regression model with pregnancy anxiety as a categorical measurement (high versus low) corrected for gender, gestational age of the newborn, batch (chip Sentrix ID) effects, and cell type compositions, similar to the DMRcate model. A list of DMPs was further used for comb-p DMR analysis. Both comb-p and DMRcate have distinct underlying statistical methods to identify DMRs. DMRcate runs a model using a kernel smoothing function on logit-transformed beta values (“M values”), independently from the single CpG site significance. Comb-p, on the contrary, identifies DMRs based on the output of a probe-by-probe analysis, i.e., using the DMP analysis p values and locations of each CpG probe. Both methods calculate DMRs only accounting for their location in the genome, i.e., non-annotated. We combined their results to rank DMRs and further use this ranking to select candidate genes. To statistically test the overlap of unique, protein-coding genes between DMRcate and comb-p lists, the R package GeneOverlap (Bioconductor) was used [26].

Based on the HM450-identified DMR, an amplicon of 446 base pairs was designed using EpiDesigner. A schematic overview and the genomic sequence of the amplicon located at chr6:29594930–29595375 (hg19; forward primer: GGTTAGGGGGTTAGGTTTGTTAGTT; reverse primer: ACTCCCTCAAAAAATCAATATCTCC) can be found in Additional file 1: Figure S2. In silico analysis, using the MassArray package in Bioconductor [27] preceded laboratory experiments in order to assess (un)successful CpG units in advance [28]. EpiTYPER DNA methylation analysis was performed for a total of 80 cord blood samples (including 45 samples used for the HM450 analysis and 35 additional samples) of the PELS cohort as described [15]. In short, cord blood DNA (1 μg) was bisulfite treated using the MethylDetector kit with long cycling protocol (Active Motif, Carlsbad, CA, USA) [29]. Subsequently, amplicons were amplified in triplicate for each sample, followed by T-cleavage reactions and detection by mass spectrometry using the Sequenom MassARRAY (San Diego, CA, USA) protocol. The signal was subsequently translated by the EpiTYPER software resulting in DNA methylation percentages for each CpG unit. Values were excluded if the standard deviation between triplicates was more than 10%, and CpG units were excluded when having a success rate of less than 70%. A summary of CpG sites analyzed by EpiTYPER and their corresponding HM450 CpG probe identifiers is shown in the Supplementary data (Additional file 2). For this part of the study, all 80 samples of our PELS cohort were used and, therefore, PRAQ could be analyzed as a continuous variable. The association of DNA methylation at each CpG unit with PRAQ was analyzed using Pearson’s correlation tests. SPSS Statistics, version 23 (IBM Corp., Armonk, NY, USA) was used for these statistical analyses.

Demographic data of the study cohort selected for EWAS

Cord blood DNA samples for the genome-wide methylation analysis were selected from the 80 mother-infant dyads of our PELS cohort. We aimed for the upper and lower quartiles of PRAQ-integrity scores and included 23 mothers with the highest and 22 mothers with the lowest scores.

Global DNA methylation is not influenced by pregnancy anxiety

Global DNA methylation difference between the low and the high prenatal anxiety group was assessed. We compared DNA methylation of all long interspersed nuclear element (LINE) repeats by analyzing preprocessed and normalized beta values at CpG probes in both LINE-1 and LINE-2 repeats extracted from the HM450 data as described previously [30]. This analysis involved 15,612 CpG probes at LINE repeats localized throughout the genome. The average DNA methylation percentage for these repetitive elements was 68.58 versus 68.61% for the high and low prenatal anxiety group, respectively.

Differentially methylated regions associated with pregnancy anxiety

Because a region of nearby CpG sites that are differentially methylated often provides stronger evidence for effective transcriptional regulation, we investigated the presence of DMRs using DMRcate. This method starts from preprocessed, normalized methylation values to analyze the HM450 data. DMRcate identified 901 DMRs with 3 or more probes that were differentially methylated between newborns exposed to high versus low pregnancy anxiety (top 50 DMRs presented in Additional file 3).

Epigenome-wide pathway analysis was performed in order to determine specific pathways, which are over-represented. For this purpose, we used ErmineJ on the top 500 uniquely annotated DMRs from DMRcate, enabling us to additionally assess multifunctionality of the resulting GO terms. The top enriched pathways are shown in Additional file 4, showing medium to high multifunctionality scores. Three of the top five GO terms relate to brain development, i.e., “pallium development” (p = 0.001), “hippocampus development” (p = 0.004), and “telencephalon development” (p = 0.004). Other top GO terms seem to be related to general cellular function and tight junctions with multifunctionality scores ranging from 0.002 to 0.932 (Additional file 4).

GABBR1 DNA methylation is positively associated with pregnancy anxiety in male newborns

We aimed at identifying a candidate gene to be further investigated in a larger sample set that includes the full PELS cohort (n = 80). For this purpose, we wanted to strengthen the DMR identification analysis by applying to our genome-wide data a second statistical method to identify DMRs, namely comb-p. The latter, in contrast to DMRcate, identifies DMRs based on the output of a probe-by-probe analysis, using the location and p values of DMPs. To this end, we first tested each CpG probe for differential methylation between neonates exposed to high or low pregnancy anxiety, while correcting for batch effects, gestational age, gender of the newborn, and cord blood cell type estimates (Additional file 5). In total, 17,180 CpG sites were differentially methylated with an unadjusted p value of < 0.05 (referred to as “DMP”) between the low and high pregnancy anxiety group. However, following the false discovery rate (FDR) correction for multiple testing none of the individual probes remained significant.

Subsequently, the comb-p tool was applied, returning ten DMRs (Table 2). Individual probes within the DMRs identified by comb-p are highlighted in a Manhattan plot (Fig. 1). The top ten DMRs according to the comb-p method are presented in Table 2 with indication of the rank order of the equivalent genomic region in the DMRcate results. There was a significant overlap between genes annotated to the DMRcate and comb-p lists, as assessed by GeneOverlap (p = 3.4E-7). The combined ranking of both comb-p and DMRcate DMR analyses identified homeobox C4 (HOXC4), palmitoyl-protein thioesterase 2 (PPT2) and GABA-B receptor subunit 1 gene (GABBR1) as top three DMRs of interest with a maximum beta fold change (FC) of 0.044 (minimum comb-p p = 1.8E-06, DMRcate p = 5.36E-14), 0.050 (minimum comb-p p = 4.5E-05, DMRcate p = 5.43E-18), and 0.053 (minimum comb-p p = 0.00015, DMRcate p = 1.08E-08), respectively. Contrary to HOXC4 and PPT2, GABBR1 DMR is located within a CpG island near active regulatory elements (Additional file 1: Figure S2) and was therefore selected for our validation study with Sequenom EpiTYPER. DNA methylation levels for each HM450 CpG probe in the GABBR1 DMR identified by DMRcate, including percent methylation difference between both groups (range 0.24 to 4) and effect sizes (range 0.170 to 0.825) can be found in Additional file 6.

Rank comb-p	Rank DMRcate	RefGene name	Location (GRCh37/hg19 assembly)	Minimum p value	Number of probes	Stouffer-Liptak-Kechris p value	Šidák corrected	CpG island
1	624	TP53INP1	chr8:95962083–95962464	0	6	0	0	yes
2	NA	PRSS50	chr3:46759334–46759699	0	8	0	0	yes
3	3	HOXC4	chr12:54446278–54446577	0	6	0	0	no
4	169	PKP3	chr11:396685–397078	0	3	0	0	yes
5	NA	ZNF764	chr16:30572738–30573014	0	5	0	0	yes
6	2	PPT2	chr6:32120954–32121421	0	14	0	0	no
7	21	GABBR1	chr6:29595001–29595316	0	7	2.70E-07	0.07	yes
8	358	PITPNM3	chr17:6358362–6358600	0	4	2.88E-07	0.1	yes
9	NA	ESRRG	chr1:217306589–-217306764	0.01	3	1.01E-06	0.39	no
10	NA	HS3ST2	chr16:22959593–22959820	0	2	0	1	yes

Both uncorrected (Stouffer-Liptak-Kechris) p values and p values corrected for multiple testing (Sidak) are shown, as well as the rank order of the DMR obtained using the DMRcate method

An overview of the GABBR1 gene and the region analyzed by Sequenom EpiTYPER platform can be found in Additional file 1: Figure S2. A summary of CpG units analyzed or excluded for technical reasons is shown in Additional file 2 and correlations between EpiTYPER and HM450 probes are shown in Additional file 7. The 17 CpG sites present in the GABBR1 DMR corresponded to 14 CpG units, of which 4 could not be analyzed due to technical reasons inherent to the MassArray and EpiTYPER methodology (e.g., low mass, fragment overlap) and one had a low success rate (CpG_4). Correlations between pregnancy anxiety and DNA methylation of the GABBR1 DMR were not significant in the total sample set (mean GABBR1 DMR methylation: Pearson R = 0.100, p = 0.379). However, as DNA methylation at several GABBR1 CpG units was associated with the pregnancy anxiety x gender interaction term (CpG_8 p = 0.045, CpG_9.10 p = 0.014, CpG_14.15 p = 0.029, and a trending effect at CpG_3 p = 0.06), data were stratified by gender. When analyzing the EpiTYPER data for this GABBR1 separately for male and female newborns, we observed a positive correlation between fetal GABBR1 DNA methylation and pregnancy anxiety in male newborns on average across the amplicon (Pearson R = 0.332, p = 0.039), at CpG_6, CpG_8, CpG_9.10, CpG_14.15, and a trend was found (p < 0.1) at CpG_1.2, CpG_16, and CpG_17 (Table 3). Only single CpG units CpG_8 (Pearson R = 0.517, p = 0.002) and CpG_14.15 (Pearson R = 0.462, p = 0.003) remained significant following Bonferroni correction for multiple tests (at α = 0.0056) (Fig. 2).

		GABBR1 DNA methylation
		CpG_1.2	CpG_3	CpG_6	CpG_8	CpG_9.10	CpG_11	CpG_14.15	CpG_16	CpG_17	Average GABBR1
Boys (n = 39)	Pearson R	0.302	0.184	0.37	0.517a	0.37	0.156	0.462a	0.295	0.299	0.33
	p value	0.061	0.306	0.043	0.002	0.023	0.395	0.003	0.095	0.086	0.039
Girls (n = 41)	Pearson R	− 0.029	− 0.288	0.032	0.035	− 0.192	− 0.128	− 0.040	− 0.078	− 0.106	− 0.049
	p value	0.860	0.093	0.858	0.844	0.278	0.451	0.824	0.655	0.531	0.761

^aSignificant after correction for tests at nine CpG sites (α = 0.0056)

GABBR1 DNA methylation significantly correlates with HPA axis function

At 12 months, salivary cortisol levels of infants were measured to calculate cortisol awakening responses. However, performing a paired t test indicated that there was no significant increase in cortisol between awakening and 30 min (mean difference = − 0.009, p = 0.647). This is also apparent from Additional file 5: Figure S3B and S3C. Since there was no significant awakening response found, association analyses with GABBR1 DNA methylation levels were not carried out.