Study population
For the BBB study, women aged 18 years or over, and less than 30 weeks pregnant were recruited through screening programs at the Northern Hospital and Mercy Hospital for Women, Melbourne, Australia, and via other health professionals and services in the public (e.g., obstetricians, GPs, and PaNDA; a Perinatal Anxiety and Depression helpline) and private sector (e.g., Northpark Private Hospital). The participating institutions were reached through advertisement and encouraged to refer women with suspected clinical depression. Women scoring 13 points or higher on the Edinburgh Postnatal Depression Scale (EPDS), the optimal score for detecting depression during pregnancy [38], were referred to the study for assessment by a psychologist if they consented. They were included in the study if they met DSM-IV criteria for a minor or major depressive disorder or an adjustment disorder with mixed depression and anxiety [39]. Severity of depression and anxiety symptoms was measured with the Beck Depression and Anxiety Inventories [40, 41]. Women with comorbid axis I disorders or medical conditions that were likely to interfere with study participation, risk requiring crisis management, participation in other psychological programs, or significant difficulty with English were excluded [24]. Women included in the study (N = 54) were randomized to receive pregnancy-specific CBT (N = 28) or TAU (N = 26). The CBT program consisted of seven individual sessions and one partner-session. TAU consisted of case-management by a midwife or a general practitioner and referral to other services of agencies as necessary. For ease of interpretation, in the results sections of this paper, the group of children of mothers from the CBT group will be referred to as the “intervention” group, and the group of children of mothers from the TAU group will be referred to as the “control” group. For participation in the current study, starting approximately 5 years after the BBB program had ended, all participants were invited through a letter. If they agreed to participate, an appointment at the Melbourne Brain Institute was planned, and informed consent was signed prior to or on the day of their visit to the clinic. If women were not able to attend the clinic, they were invited to send a buccal sample through the mail. The study was approved by the Human Research Ethics Committees of Austin Health, Melbourne, Australia.
Data collection
A questionnaire on current sociodemographic data and current symptoms of depression and anxiety was sent to each woman’s home address. Baseline demographics, including symptoms of depression and anxiety as well as the child’s birth weight, were taken from the BBB study files. At the Melbourne Brain Centre, a cognitive assessment by means of the Wechsler Preschool and Primary Intelligence Scale (WWPSI-III) [42] was performed on the child, an MRI scan of the child’s brain was conducted, of which results are described elsewhere, and a buccal cell sample from the child was obtained by a researcher who was blinded to the allocation status of the women.
Buccal cell samples
Buccal cells were collected using a dedicated swab (OraCollect 100, DNA Genotek Inc., Ontario, Canada). Children were instructed not to eat or drink 30 min prior to taking the swab. Women who were not able to visit the Melbourne Brain Centre were instructed how to apply the swab on their child, and asked to send the sample via mail. The swabs were stored at room temperature at the Parent-Infant Research Institute and transported to the Murdoch Children’s Research Institute (Melbourne, Australia) for DNA extraction within 2 weeks after collection.
DNA extraction and genome-wide methylation detection
DNA extraction of all samples was performed using the NucleoBond CB20 DNA extraction kit. Purification of DNA was assessed using Nanodrop Spectrophotometry. Bisulfite conversion was performed using the EZ-96 DNA methylation kit (ZYMO Research Corporation) according to the manufacturer’s instructions. DNA methylation profiling was performed at the Australian Genome Research Facility, on bisulfite converted DNA using the Illumina Infinium Methylation EPIC BeadChip Array (HM850) (Illumina), which measures CpG methylation at > 850,000 genomic sites.
Candidate gene approach
We extracted 729 probes spanning 16 a priori selected genes for linear regression analysis. Candidate genes were those that had previously been assessed in relation to prenatal exposure to maternal stress, depression, and/or anxiety in earlier studies [20]. Genes of interest were genes encoding brain-derived neurotrophic factor (BDNF; 91 probes), corticotrophin releasing hormone (CRH; 21 probes), corticotrophin-releasing factor-binding protein (CRHBP; 25 probes), corticotrophin-releasing hormone receptors 1 and 2 (CRHR1; 41 probes, CRHR2; 40 probes), FK506 binding protein (FKBP5; 49 probes), a long noncoding RNA (H19; 57 probes), hydroxysteroid 11-beta dehydrogenase 1 and 2 (HSD11B1; 25 probes, HSD11B2; 23 probes), insulin-like growth factor (IGF2; 15 probes), maternally expressed 3 (MEG3; 87 probes), mesoderm-specific transcript homolog protein (MEST; 63 probes), the glucocorticoid receptor (NR3C1; 89 probes), the mineralocorticoid receptor (NR3C2; 50 probes), the oxytocin receptor (OXTR; 22 probes), and the serotonin transporter (SLC6A4; 31 probes) [20]. Additionally, considering the especially strong evidence for this gene, we separately analyzed the probes of the promoter region of the glucocorticoid receptor gene (NR3C1 promoter-associated probes; 34 probes) for differential methylation.
Statistical analysis
DNA methylation was defined as a continuous variable varying from 0 (completely unmethylated) to 1 (completely methylated). Methylation data were processed in R using the minfi package. Normalization of the data was performed using the SWAN method [43]. Probes on X and Y chromosomes, probes that were associated with SNPs with a minor allele frequency > 1%, and cross-reactive probes [44] were removed from the dataset. This resulted in data for 770,668 probes available for subsequent analysis.
Sources of variation
Main contributors to the variation in the methylation data were identified by principal component analysis (PCA). We included the following variables in the analysis to assess associations with PC’s: participant ID, chip ID, HM850 array chip position, allocation, sex, child age, birth weight, maternal age, gestational age, current income, baseline depression symptoms, baseline anxiety symptoms, current depression symptoms, and current anxiety symptoms. Results of the PCA showed that the first five principal components contributed most to the variation in the methylation data, and all variables associated with any of these PC’s were added as covariate in all analyses (Fig. 3a). The heatmap demonstrated that allocation was associated with the third principal component. Birth weight, child age, sex, and HM850 array chip position were associated with the first four principal components and they were included in the analyses as covariates. None of the variables included in our model was significantly associated with the fifth principal component, and this PC was therefore not included in our model as covariate (Fig. 3b). Unsupervised analysis by multidimensional scaling was conducted in order to examine sources of variation within the dataset. Beta values (methylation level) at all HM850 probes for all samples were used to produce multidimensional scaling (MDS) plots, with samples colored according to intervention (turquoise)/control (orange) status, showing the relatedness of samples over the first two principal components of variation (Fig. 4a). Coloring by intervention/control revealed no distinct separation by allocation. Additional MDS plots of samples over other principal components also failed to show a distinct separation between the two groups (Figs. 4b c).
Differential methylation according to allocation
Linear regression analysis was used to identify associations between the intervention status and epigenome-wide DNA methylation. We took into account variation associated with the covariates birth weight, HM850 array chip position, child sex and age, to account for PC1, PC2, PC3, and PC4, as identified by PCA. The Benjamini-Hochberg False-Discovery-Rate method [45] was used to correct for multiple testing. However, none of the analyses yielded significant differentially methylated probes between the intervention and control group after correcting for multiple testing. In an explorative analysis, we extracted differentially methylated probes between the intervention and control group at a nominal significance level set at p < 0.01, prior to correcting for multiple testing. We assessed differences in mean DNA methylation of all significant probes between the intervention and control group using an unpaired Mann-Whitney-Wilcoxon test. We additionally compared mean beta differences of 16 candidate genes, and the promoter region of the NR3C1 gene between the intervention and control group using an unpaired Mann-Whitney-Wilcoxon test.
Differential methylation according to baseline depression or anxiety symptom score
As additional explorative analyses, two separate linear regression models were also used to investigate associations between baseline depression (BDI–II score) and baseline anxiety (BAI- score) with methylation profiles in the children. For ease of interpretation, the sample was divided into two groups in both analyses. The rationale behind this approach was to explore widespread methylation variation between women with severe symptoms compared to those with mild symptoms using clinically relevant cut-offs, rather than investigating the direction of correlations between increasing depression and anxiety scores on all probes separately. Baseline depression was converted to a dichotomous variable using clinically relevant Beck questionnaire cut-offs. Women with BDI-II ≥ 29 were classified as “highly depressed” (n = 13), whereas those with a score below 29 were classified as “mildly depressed” (n = 9) [46]. This procedure was repeated for baseline anxiety (BAI-score). The cut-off for clinically relevant anxiety is set at 16, and therefore we classified women with BAI ≥ 16 as “highly anxious” (n = 8), and women with BAI below 16 as “mildly anxious” (n = 14) [47]. One woman had missing data on baseline depression and anxiety and was excluded from the analysis. We took into account allocation status, birth weight, HM850 array chip position, child sex, and age as covariates, as identified by PCA. Differentially methylated probes at a nominal significance level set at p < 0.01, prior to correction for multiple testing, were extracted. We compared differences in mean DNA methylation in groups of children of women with high baseline symptoms and low baseline symptoms using an unpaired Mann-Whitney-Wilcoxon test, both for depression and anxiety. We additionally compared mean beta differences of 16 candidate genes, and the promoter region of the NR3C1 gene between groups of children of women with high baseline symptoms and low baseline symptoms using an unpaired Mann-Whitney-Wilcoxon test, both for depression and anxiety.