Prenatal medication exposure and epigenetic outcomes: a systematic literature review and recommendations for prenatal pharmacoepigenetic studies

ABSTRACT

When used during pregnancy, analgesics and psychotropics pass the placenta to enter the foetal circulation and may induce epigenetic modifications. Where such modifications occur and whether they disrupt normal foetal developme nt, are currently unanswered questions. This field of prenatal pharmacoepigenetics has received increasing attention, with several studies reporting associations between in utero medication exposure and offspring epigenetic outcomes. Nevertheless, no recent systematic review of the literature is available. Therefore, the objectives of this review were to (i) provide an overview of the literature on the association of prenatal exposure to psychotropics a nd analgesics with epigenetic outcomes, and (ii) suggest recommendations for future studies within prenatal pharmacoepigenetics. We performed systematic literature searches in five databases. The eligible studies assessed human prenatal exposure to psychotropics or analgesics, with epigenetic analyses of offspring tissue as an outcome. We identified 18 eligible studies including 4,419 neonates exposed to either antidepressants, antiepileptic drugs, paracetamol, acetylsalicylic acid, or methadone. The epigenetic outcome in all studies was DNA methylation in cord blood, placental tissue or buccal cells. Although most studies found significant differences in DNA methylation upon medication exposure, almost no differences were persistent across studies for similar medications and sequencing methods. The reviewed studies were challenging to compare due to poor transparency in reporting, and heterogeneous methodology, design, genome coverage, and statistical modelling. We propose 10 recommendations for future prenatal pharmacoepigenetic studies considering both epidemiological and epigenetic perspectives. These recommendations may improve the quality, comparability, and clinical relevance of such studies. PROSPERO registration ID: CRD42020166675.

KEYWORDS:

• Epigenetics
• DNA methylation
• prenatal exposure
• psychotropics
• analgesics
• literature review
• EWAS
• epigenetic epidemiology
• epidemiology
• recommendations

Background

Every day, pregnant women use medications for which the scientific evidence on foetal safety is limited or inconclusive. As most medications pass both the placental and blood-brain barriers during gestation [1–4], common medications such as analgesics and psychotropics may exhibit pharmaceutical effects in the foetus and potentially disrupt normal foetal development. This reasoning is based on the Developmental Origins of Health and Disease (DOHaD) hypothesis, which is a conceptual framework linking prenatal environmental exposures to health and disease in later life [5–9]. Indeed, many studies have reported a variety of adverse developmental outcomes as sociated with in utero medication exposure, including developmental delays and abnormalities (comprehensively reviewed in [10–18] and the textbook by Schaefer et al. [19]).

T he mechanisms by which prenatal exposure to medications impacts foetal development remain largely unknown. One suggested mechanism is the direct or indirect influence of epigenetic modifications in the developing foetus [9,20]. Epigenetics encompasses regulatory mechanisms that can impact genome stability and gene transcription, such as histone modifications and DNA methylation (DNAm) of cytosine-phosphate-guanine sites (CpGs). Such modifications are reversible and can be influenced by both genetics and environmental factors, such as medications [21], making epigenetic changes plausible mediators of the prenatal environmental impact on developmental outcomes [9,20].

The research on epigenetic modifications in neonates exposed to medications in utero, hereafter referred to as prenatal pharmacoepigenetics, has gained increasing attention in recent years. Although the literature on prenatal pharmacoepigenetics is growing, only one systematic review summarizing the findings on medications that potentially interfere with foetal development is available [22]. However, this review only included studies on antidepressants. Therefore, the primary aim of the current review is to provide an overview of the literature on the association of prenatal exposure to psychotropics and analgesics with epigenetic outcomes. In addition, by evaluating the eligible studies from both epidemiological and epigenetic perspectives, this review also aims to provide recommendations for future prenatal pharmacoepigenetic research to improve the overall quality, comparability, and clinical relevance of prenatal pharmacoepigenetic association studies.

Methods

Search strategy

Literature searches were performed in the MEDLINE, EMBASE, PsycINFO, Scopus, and Web of Science databases. The searches were first completed on 19 January 2020, and any new studies meeting the eligibility criteria, published before 1 September 2020, were included in the final review. In addition, the reference lists of the eligible articles and references of 35 relevant reviews were screened to ensure complete coverage of the literature. Prior to performing the literature searches, a detailed search strategy and vocabulary were developed with support from experienced librarians in medicine, pharmacy, and psychology. We included studies investigating (i) prenatal exposure to (ii) psychotropics and analgesics with (iii) an epigenetic outcome. The search terms for these three criteria are listed in Supplementary Table S1. Supplementary Table S2 provides an example of a search in EMBASE. The review is reported in adherence to the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines [23], and the protocol and search strategy are available in the PROSPERO database (registration ID: CRD42020166675) [24,25].

Inclusion criteria

The studies included in this review were selected based on the participants, intervention/exposure, comparison group, outcome, and study design (PICOS) criteria [26]. Participants were defined as children (<18 y old) prenatally exposed to psychotropics or analgesics for which epigenetic data were available. Anatomical Therapeutic Chemical (ATC) codes were used to identify medication groups in accordance with the World Health Organization ATC index [27]. The exposure was defined as use of antidepressants (ATC code: N06A), psycholeptics (N05), antiepileptic drugs (AEDs; N03), analgesics (N02), or non-steroidal anti-inflammatory drugs (NSAIDs; M01A) during pregnancy. We specifically selected analgesics and psychotropics, based on the expertise of our research group, biological plausibility, and the emerging number of pharmacoepigenetic studies on analgesics and psychotropics. The comparison group included children of mothers who did not use the medication of interest during pregnancy. The outcome was epigenetic measurements in tissue samples from exposed and unexposed children (<18 y old). If the study also included data on immediate or long-term developmental outcomes in the children, we reported these as well. Studies investigating the same data sets were all eligible if they reported on different exposures and/or outcomes. Only original articles with the study designs case-control, cohort, or randomized controlled trial were included. No limitations were applied regarding the time of publication, but only articles in English or Scandinavian languages were eligible.

Data extraction

After searching and retrieving the results from the databases, any duplicates were removed in EndNote X8.2 and the remaining records uploaded to the online systematic review data management platform Rayyan [28]. Two reviewers (KG and EWO) independently screened the titles and abstracts, excluding studies that did not meet the inclusion criteria. If the eligibility of a paper was unclear based on the title and abstract, it was included for the next round of screening. In the second screening, the full-text versions of all papers were read and the final exclusion of papers performed. Any disagreement between the two reviewers was resolved by a third reviewer (HMEN).

Results

Outcomes of the screening and selection process

The initial searches yielded a total of 2,159 records: 488 records in MEDLINE, 880 records in EMBASE, 88 records in PsycINFO, 194 records in Scopus, and 509 records in Web of Science (Figure 1). A total of 871 duplicated records were removed, leaving 1,288 unique articles to screen the titles and abstracts in Rayyan [28]. After the first screening, 1,262 papers were excluded due to being non-original studies (n = 605) or failing to meet the defined PICOS criteria (n = 657). After reading the complete texts of the 26 records remaining from the first round, we excluded 11 records due to wrong exposure according to our criteria (could not differentiate medication exposure across groups; n = 5), wrong population according to our criteria (participants were too old upon exposure or sampling; n = 3), or wrong comparison groups according to our criteria (did not include a non-medicated group; n = 3). By screening the reference lists of the 15 remaining records and 35 relevant reviews, we identified 1 additional article. Two additional studies meeting the eligibility criteria were published during the revision process (before 1 September 2020) and were also included. Consequently, a total of 18 records were included in the final review.

Figure 1. Flow chart of article screening and selection based on the template from PRISMA [23]. ‘Second search’ refers to eligible studies published during the manuscript revision process.

Overview of the eligible studies

All of the eligible articles were based on data from single birth cohorts, except for one study validating results in an independent cohort [29], and one study being a randomized controlled trial [30]. Of the eligible articles, nine were epigenome-wide association studies (EWASs; median sample size 241 neonates [interquartile range; IQR: 284]), and eleven were candidate gene studies (median sample size 115 neonates [IQR: 168]). Hence, two studies, including 46 and 58 neonates, combined epigenome-wide and candidate gene approaches [31,32]. The medications included were the psychotropics antidepressants (12 studies; median sample size 201 neonates [IQR: 354]), and AEDs (2 studies; 18 and 201 neonates), and the analgesics paracetamol (2 studies; 281 and 384 neonates), acetylsalicylic acid (1 study; 358 neonates), and methadone (1 study; 53 neonates). The epigenetic outcome investigated in all papers was DNAm, in cord blood (13 studies; median sample size 201 neonates [IQR: 341]), placental tissue (5 studies; median sample size 236 neonates [IQR: 38]), and/or buccal cells (2 studies; 236 and 53 neonates). The neonatal tissues were sampled within 72 h after birth in all studies, except by Cardenas et al. (2019), who also collected blood from children aged 3–5 y and 7–11 y [29]. All studies adjusted for potential covariates and/or confounders in their statistical analyses or by design, but the number of variables under consideration differed greatly (Supplementary Table S3). The covariates most frequently accounted for were maternal age (n = 16), smoking during pregnancy (n = 13), infant sex (n = 12), gestational age (n = 10), and folate use in pregnancy (n = 10).

In addition to the epigenetic outcomes, several studies reported phenotypic outcomes in children, specifically poor foetal growth (n = 1) [33], birth weight (n = 2) [31,34], severity of neonatal abstinence syndrome (n = 1) [34], ADHD (n = 1) [35], stress reactivity (n = 2) [36,37], and soothability (n = 1) [38] (Supplementary Table S4). One study performed a mediation analysis of medication exposure, epigenetic modification, and neonatal phenotypic outcome [37]. This study assessed whether the DNAm of a CpG in the placental NR3C2 gene acted as a mediator of the effect of maternal depressive symptoms on cortisol reactivity in 12-month-old infants [37]. The effects of maternal depression on cortisol levels were decomposed into direct effects and DNAm-mediated indirect effects, finding that, although the indirect effect of DNAm was positive, it did not overcome the larger negative direct effect of depressive symptoms on infant cortisol levels [37]. However, the analysis demonstrated an increased DNAm at the NR3C2 CpG upon in utero antidepressant exposure, suggesting that maternal antidepressant use during pregnancy enhances the indirect effect of NR3C2 DNAm on the infant stress response [37].

All EWASs used the Illumina platform [39] to assess DNAm with the MethylationEPIC (n = 2), HumanMethylation 450 (n = 3), or HumanMethylation 27 (n = 4) bead chips. To assess the association between CpG DNAm and medication exposure, the majority of the EWASs used linear regression models (n = 6). In most of the EWASs, a result was considered significant if the false discovery rate (FDR) adjusted p-value was <0.05 (n = 8), except for one study that used an FDR adjusted p < 0.1 [32].

In the candidate gene studies, several methods were used to investigate DNAm: the Illumina platform (n = 1) [33], the SEQUENOM MassARRAY EpiTYPER platform (n = 3), and the PyroMark system (n = 7). The studies reported the methylation percentages (n = 6), mean methylation percentages of triplicates (n = 4) or the β value from the Illumina microarray (n = 1) [33]. Various statistical tests were applied to assess differential DNAm. For these tests, three studies used FDR adjusted p-values (p < 0.25 in one study; p < 0.05 in two studies), two studies used the Bonferroni-corrected p-value, and the last six studies applied an unadjusted p < 0.05. Further details on the studies are available in Tables 1 and 2. For excellent discussions of statistical approaches in epigenetic studies, we recommend the recently published reviews by Teschendorff and Relton [40], van Rooij et al. [41], and Mansell et al. [42].

Table 1. Overview of the studies included in the literature review

	Reference	Country; setting (study period)	Sample sizen(groups)	Research aim(s)	Exposure(s) (data source)	Pharmacoepigenetic outcome(s)
Epigenome-wide association studies	Yeung et al. (2020) [30]	USA; the EAGeR (Effects of Aspirin in Gestation and Reproduction) randomized trial (2006–2012)	n = 358 (acetylsalicylic acid, n = 185; no acetylsalicylic acid, n = 173)	Investigate the impact of maternal use of low-dose acetylsalicylic acid prior to and during pregnancy on cord blood DNAm	Acetylsalicylic acid (randomly assigned prior to conception; 81 mg/d until gestational week 36 [49])	• Differential DNAm at 1 CpG associated with prenatal acetylsalicylic acid exposure (cg2002882; 3,500 bp upstream of the POU4F1 promoter)
	Addo et al. (2019) [48,100]	USA; Extremely Low Gestational Age Newborns (2002–2004)	n = 281 (paracetamol, n = 165; no paracetamol, n = 116)	Examine DNAm of CpGs in placentae of paracetamol-exposed and non-exposed neonates	Paracetamol (≥ 1 during pregnancy; self-reported Tylenol use)	• Differential DNAm at 24 CpGs associated with prenatal paracetamol exposure
	Cardenas et al. (2019) [29]	USA; Project Viva (1999–2002) The Netherlands; Generation R Study (2002–2006)	nProV = 479 (antidepressant^a, n = 14; depression, n = 33; anxiety, n = 40) nGenR = 999 (antidepressant, n = 14; depression, n = 31; anxiety, n = 56)	Identify DNAm differences in neonates associated with exposure to maternal anxiety, depression, or antidepressant use in pregnancy	Antidepressants (≥ 1 prescription in pregnancy; medical record in Project Viva, self-reported and prescription-validated in Generation R) Maternal depression (EPDS at mid-pregnancy in Project Viva; BSI at 20 weeks in Generation R) Maternal anxiety (PRAS at mid-pregnancy in Project Viva; BSI at 20 weeks in Generation R)	• Differential DNAm at 130 CpGs in Project Viva in neonates prenatally exposed to antidepressants compared to non-exposed neonates, 5 confirmed in Generation R (1 under Bonferroni significance; reduced DNAm on cg22159528 ZNF575 of exposed children) • No DMRs associated with prenatal antidepressant exposure
	Gervin et al. (2017) [35]	Norway; the Norwegian Mother, Father and Child Cohort Study (1999–2008)	n = 384 (no paracetamol & no ADHD, n = 96; paracetamol & no ADHD, n = 96; no paracetamol & ADHD, n = 96; short-term paracetamol & ADHD, n = 77; long-term paracetamol & ADHD, n = 19)	Investigate if differential DNAm is associated with prenatal paracetamol exposure and ADHD development	Paracetamol (long-term [≥ 20 d] and short-term [6–19 d] exposure; self-reported in questionnaires during pregnancy) ADHD (diagnosis in the Norwegian Patient Registry)	• In children with ADHD, prenatal long-term exposure to paracetamol was associated with differential DNAm compared to children without ADHD not exposed to paracetamol (6211 CpGs), children with ADHD and not exposed to paracetamol (193 CpGs), and short-term paracetamol-exposed children with ADHD (2089 CpGs)
	Emes et al. (2013) [46]	UK; EFFECT-M study	n = 18 (epilepsy & AEDs^b, n = 9; no epilepsy & no AEDs, n = 9)	Examine association between prenatal AED exposure and DNAm, and if AEDs affect the foetal DNAm by lowering the maternal folate level	AEDs (self-reported and validated by medical record)	• DNAm difference in 662 CpGs (652 different genes) in AED-exposed compared to non-exposed neonates • No difference in global DNAm levels between AED-exposed and non-exposed children
	Smith et al. (2012) [47]	USA; the Emory Women’s Mental Health Program	n = 201 (AEDs^c & epilepsy/psychiatric disorder^d, n = 53; no AEDs or epilepsy/psychiatric disorder, n = 148)	Examine the impact of prenatal AED exposure on DNAm patterns in neonates	AEDs (self-reported every 4–6 weeks in pregnancy and validated by concentrations of AEDs in maternal blood) Psychiatric disorder and/or epilepsy (questionnaire on medical and psychiatric history and SCID at intake, SCID and seizure history at 4–6-week intervals)	• Prenatal AED exposure associated with decreased global DNAm in cord blood, not in placenta • Longer prenatal AED exposure associated with decreased DNAm in 14 cord blood CpGs • In placental tissue, 3 of the 14 cord blood CpGs also exhibited decreased DNAm (PGC, ZNF384, and C15orf2) • DNAm patterns neither specific to AED type nor more extreme differences for polydrug treatment
	Schroeder et al. (2012) [43]	USA; the Emory Women’s Mental Health Program	n = 201 (Several different comparisons^†: • Current MDD, n = 118; no current MDD, n = 83 • Antidepressants^e, n = 151; no antidepressants, n = 50)	Investigate the association of maternal psychiatric disorder, symptoms and severity of depression, and medication treatment in pregnancy with neonatal DNAm patterns	Psychotropics (medication evaluation upon visits every 4–6 weeks during pregnancy) Maternal mood disorder diagnosis (life-time history and MDEs; SCID every 4–6 weeks during pregnancy, SCID Mood Module assessed MDEs) Depressive symptoms (depression severity and clinically significant depressive symptoms; HRSD17 and BDI)	• Prenatal antidepressant exposure associated with DNAm in 2 CpGs (1.9% decrease in TNFRSF21 and 3% increase in CHRNA4), independent of antidepressant type and duration • No association between prenatal hypnotic, antiemetic, benzodiazepine, or atypical antipsychotic exposure and differential DNAm, independent of duration of exposure • Prenatal exposure to atypical antipsychotics associated with DNAm at 1 CpG
Combined epigenome-wide and candidate gene studies	Gurnot et al. (2015) [31]	Canada; University of British Columbia	nepigenome-wide = 23 (depressed & SRIs^f, n = 11; depressed & no SRIs, n = 12) ncandidate gene = 44 (depressed & SRIs^g, n = 19; depressed & no SRIs, n = 25)	Examine if prenatal SRI exposure and/or maternal mood is associated with DNAm across the genome and in CYP2E1; investigate if DNAm is also associated with birth outcomes	SRIs (cord blood/maternal blood drug ratios at birth; maternal whole blood and neonatal cord blood concentrations of SRIs, drawn directly after birth) Maternal mood (HAM-D at 26 and 36 weeks of pregnancy, mean score)	• Prenatal SRI exposure associated with DNAm in 3 CpGs (CYP2E1, EVA1, and SLMAP; EWAS) • In CYP2E1, DNAm highly negatively correlated with prenatal maternal depressive mood only if prenatally exposed to SRIs concurrently • Pyrosequencing of CYP2E1 yielded DNAm values that correlated with the microarray findings
	Non et al. (2014) [32]	USA; Harvard Epigenetic Birth Cohort (2007–2009)	n = 58 (SSRIs^h, n = 22; depression/anxiety & no SSRIs, n = 13; no depression/anxiety or SSRIs, n = 23)	Examine differences in DNAm patterns across the genome in neonates prenatally exposed or non-exposed to SSRIs and/or maternal depression/anxiety	SSRIs (medical charts) Maternal depression/anxiety (explicitly noted by obstetrician in labour and delivery forms)	• No association between prenatal SSRI exposure and differential DNAm (EWAS) • No regional clusters of CpGs (1 kb) associated with prenatal SSRI exposure (EWAS) • Pyrosequencing: 6 CpGs in Col7a1 exhibited lower DNAm at all 6 CpGs and a lower mean DNAm in neonates prenatally exposed to SSRIs • Pyrosequencing: prenatal SSRI exposure associated with 1 CpG in NFKB2, 1 CpG in SLC6A4, 1 CpG in FKBP5, and 1 CpG in DNMT3a
Candidate gene studies	Galbally et al. (2020) [37]	Australia; Mercy Pregnancy and Emotional Wellbeing Study (2012–2015)	n = 236 (antidepressants, n = 43; non-medicated & depression/dysthymia, n = 24; non-medicated & no depression/dysthymia, n = 169)	Investigate associations between maternal depression during pregnancy and the DNAm of placental and buccal NR3C1 and NR3C2, which may mediate an indirect effect of maternal depression on 12-month infant cortisol reactivity	Antidepressants (self-reported and hospital records, validated by concentration in cord and maternal whole blood, converted to SEDs) Maternal mental health (SCID at ≤ 20 weeks, EPDS/STAI in week 20, 3^rd trimester, and 6 and 12 months after birth) Maternal stress (PRAMS at week 20 and in 3^rd trimester)	• 1 differentially methylated CpG in placental NR3C2 when comparing medicated and non-medicated depression • No differential methylation of NR3C1 or NR3C2 in buccal cells
	Galbally et al. (2018) [51]	Australia; Mercy Pregnancy and Emotional Wellbeing Study (2012–2015)	n = 239 (untreated current MDD, n = 24; antidepressant-treated current MDD, n = 28; antidepressant-treated not meeting MDD diagnostic criteria, n = 15; no current or past MDD, n = 172)	Explore DNAm of OXTR in the placentae of women depressed during pregnancy and in women using antidepressant(s) during pregnancy	Antidepressants (self-reported at recruitment and in 3^rd trimester, hospital records in 3^rd trimester [converted to SEDs], and measurement of whole blood and cord blood concentrations at birth) Maternal depression (SCID at recruitment and EPDS in 3^rd trimester)	• Decreased DNAm in OXTR CpG 8 upon foetal exposure to antidepressants (self-reported) • Increased DNAm of OXTR CpG 8 upon higher cord blood concentrations of antidepressants
	McLaughlin et al. (2017) [34]	UK; Princess Royal Maternity Hospital	n = 53 (methadone-maintained opioid-dependent mothers, n = 21; smoking, “deprived” & opioid-naive mothers, n = 17; non-smoking, “affluent” & opioid-naive mothers, n = 15)	Explore if prenatal opioid exposure is associated with a differential DNAm in opioid-related genes (ABCB1, CYP2D6, and OPRM1)	Methadone (venous blood concentration 24–72 hours after birth to methadone-maintained opioid-dependent mothers) Smoking (case records) Poverty (DepCat score calculated from postal codes provided in case records)	• Increased DNAm in ABCB1, CYP2D6, and OPRM1 in neonates of methadone-maintained opioid-dependent mothers compared to neonates of opioid-naive mothers
	Mansell et al. (2016) [45]	Australia; The Barwon Infant Study (2010–2013)	n = 481 (Various group comparisons: • Depressive symptoms, n = 88; • No depressive symptoms, n = 357 • Anxiety, n = 77; no anxiety, n = 368 • Stress scores, n = 481)	Investigate the association between maternal mental well-being and the DNAm of cord blood NR3C1	Antidepressants (self-reported at 28 weeks) Depression (EPDS ≥ 10 at 28 weeks) Anxiety (EPDS anxiety subscale ≥ 5 at 28 weeks) Stress (PSS score week 28)	• Increased DNAm of NR3C1 CpG 35 associated with prenatal antidepressant exposure, but the association diminished when included as a covariate in the multivariate model of maternal pregnancy well-being and neonatal DNAm of NR3C1
	Gartstein et al. (2016) [38]	Canada; University of British Columbia	n = 115 (SSRI exposure, n = 46; no SSRI exposure, n = 69)	Examine the association between prenatal SSRI exposure and neonatal SLC6A4 DNAm, and the influence on soothability	SSRIs (self-reported prescription retrieval and number of days used) Maternal internalizing symptoms (sum of EPDS, HAM-D, and HAM-A scores in 33–36 weeks of pregnancy)	• Prenatal SSRI exposure positively associated with neonate SLC6A4 CpG DNAm of CpGs 3, 5, 7, and 9 and mean DNAm status of CpGs 9 and 10
	Ciesielski et al. (2015) [33]	USA; Women and Infants Hospital in Providence Rhode Island (2008–2010)	n = 184 (psychiatric diagnosis^i & antidepressants^j, n = 12; psychiatric diagnosis & no antidepressants, n = 13; no psychiatric diagnosis & no antidepressants, n = 159)	Investigate how DNAm in placentae is related to growth restriction observed in mothers with psychiatric illnesses (some of which are treated with antidepressants)	Antidepressants (medical records) Maternal psychiatric disease (depression, anxiety, and/or OCD/panic disorder prior to and/or during pregnancy; medical records)	• No difference in likelihood of unequal median DNAm between the antidepressant-exposed group, and the group with no psychiatric diagnosis and no antidepressants
	Soubry et al. (2011) [52]	USA; the Newborn Epigenetics Study (2005–2008)	nIGF2 = 356 (antidepressants^k, n = 35; depressed, n = 56; non-medicated non-depressed, n = 265) nH19 = 411 (antidepressants, n = 43; depressed, n = 65; non-medicated, non-depressed, n = 303)	Examine the association between prenatal antidepressant and/or maternal depression exposure and DNAm of two DMRs in IGF2 (H19)	Antidepressants (any use; medical charts) Depression during pregnancy (self-reported and interviewer-based questionnaires, diagnosis validated by clinical charts)	• Prenatal antidepressant exposure not associated with DNAm in neonatal IGF2 or H19 DMRs • Higher DNAm in H19 DMR associated with prenatal antidepressant exposure in African-Americans but not in Caucasians
	Devlin et al. (2010) [44]	Canada; cohort part of a study on how psychotropic medication exposure impact neonatal health	n = 82 (SRIs, n = 33; no SRIs, n = 49)	Examine the impact of maternal MTHFR C677T genotype on maternal mood, and the association with DNAm in maternal and neonatal SLC6A4 and BDNF	SRIs (data source not stated) Depressed mood (EPDS and HAM-D in the 2^nd [week 26] and 3^rd [week 33] trimesters)	• SRI exposure not associated with neonatal DNAm levels in SLC6A4 or BDNF
	Oberlander et al. (2008) [36]	Canada; cohort part of a study on how psychotropic medication exposure impact neonatal health	n = 82 ^l (depression & SRIs^m, n = 36; depression & no SRIs, n = 13; no depression & no SRIs, n = 33)	Investigate any association between maternal depressed or anxious mood during pregnancy and DNAm of NR3C1 in neonates; examine the association between NR3C1 DNAm and stress reactivity at 3 months	SSRIs (data source not stated) Maternal mood (EPDS, HAM-A and HAM-D in week 26 and 33; EPDS, HAM-A, HAM-D, PSI-SF 3 months after birth)	• Prenatal SRI exposure not associated with neonatal NR3C1 DNAm status • Global neonatal DNAm did not differ between SRI exposed and non-exposed infants

^†Complete list of comparisons is available in Supplementary Table S4.

^a1 amitriptyline; 1 bupropion; 2 citalopram hydrobromide; 2 desipramine; 3 fluoxetine; 2 paroxetine; 5 sertraline (12/14 women used SSRIs only).

^b4 carbamazepine; 3 lamotrigine; 2 polytherapy (1 carbamazepine/valproic acid; 1 lamotrigine/valproate).

^c36 lamotrigine; 3 valproate; 3 levetiracetam; 2 carbamazepine; 1 topiramate; 1 phenytoin; 1 gabapentin; 6 polytherapy (not specified).

^d26 epilepsy; 27 psychiatric disorder (20 bipolar disorder; 5 MDD; 2 anxiety).

^eDefine two classes: class 1 of SRIs (SSRIs; SNRIs; TCAs), and class 2 of bupropion.

^f2 paroxetine; 2 fluoxetine; 2 sertraline; 1 citalopram; 4 venlafaxine.

^g2 paroxetine; 3 fluoxetine; 2 sertraline; 2 citalopram; 10 venlafaxine.

^h11 sertraline; 6 fluoxetine; 4 citalopram; 2 paroxetine.

ⁱDepression, anxiety, obsessive-compulsive disorder, or panic disorder.

^j18 SSRIs (13 sertraline); 5 atypical antidepressants.

^kSSRIs (72%); SNRIs; TCAs; SARIs; bupropion.

^lNumbers not fixed; some mothers in untreated, non-depressed group became depressed during the study, whereas others started receiving pharmacological treatment.

^m18 paroxetine; 6 fluoxetine; 5 sertraline; 2 venlafaxine; 5 citalopram.

ADHD: attention-deficit/hyperactivity disorder; AED: antiepileptic drug; BDI: Beck Depressive Inventory; BSI: Brief Symptom Inventory; DepCat: measure of socioeconomic status in Scotland (affluent–deprived); DMR: differentially methylated region; DNAm: DNA methylation; EPDS: Edinburgh Postnatal Depression Scale; GO: Gene Ontology; HAM-A: Hamilton Rating Scale for Anxiety; HAM-D: Hamilton Rating Scale for Depression; HRSD17: 17-item Hamilton Rating Scale for Depression; IBQ: Infant Behaviour Questionnaire; MDD: major depressive disorder; MDE: major depressive episode; OCD: obsessive-compulsive disorder; PRAS: Pregnancy-related Anxiety Scale; PSI-SF: Parenting Stress Index – Short Form; PSS: Perceived Stress Scale; SARI: serotonin antagonist and reuptake inhibitor; SCID: Structured Clinical Interview for DSM-IV; SED: sertraline-equivalent dosage; SNRI: serotonin and noradrenaline reuptake inhibitor; (S)SRI: (selective) serotonin reuptake inhibitor; STAI: State-trait Anxiety Inventory; TCA: tricyclic antidepressant.

Table 2. Overview of the methodology and statistical analysis

	Reference	Sample tissue	Individual sites, regional or global level	Method	Epigenetic outcome	Statistical tests	Statistical thresholds	Model covariates^†	Cell count correction	Authors’ inference about exposure–outcome associations^*
Epigenome-wide association studies	Yeung et al. (2020) [30]	Cord blood	Site-by-site	Infinium MethylationEPIC BeadChip (Illumina)	Subset quantile normalized β values	Linear mixed-effects models	FDR cut-off < 0.05	6	Reference data for cord blood [101]	Association
	Addo et al. (2019) [48,100]	Placenta	Site-by-site	Infinium MethylationEPIC BeadChip (Illumina)	M values for statistical procedures, β values for data presentation	Robust linear regression models	FDR cut-off < 0.05 (2 CpGs reached Bonferroni significance at p < 6.3 × 10^−8)	14 (and PCA)	Houseman [102]	Association
	Cardenas et al. (2019) [29]	Cord blood	Site-by-site and regional	Infinium Human-Methylation450 BeadChip (Illumina)	M values for statistical procedures, β values for data presentation	Robust linear regression models	FDR cut-off < 0.05 (16 CpGs reached Bonferroni significance at P < 1.34 × 10^−7 in Project Viva, 1 site of which was confirmed in the replication cohort)	10 (and PCA)	Houseman [102]	Association
	Gervin et al. (2017) [35]	Cord blood	Site-by-site and regional	Infinium Human-Methylation450 BeadChip (Illumina)	M values for statistical procedures, β values for data presentation	Linear regression models	FDR cut-off < 0.05	6 (and SVA)	Houseman [102] and reference data for cord blood [103]	Association
	Emes et al. (2013) [46]	Cord blood	Site-by-site and global (LINE-1)	Infinium Human-Methylation27 BeadChip (Illumina)	β values (log2-transformed upon t testing)	Hierarchical clustering of the β values of all sites, then t-tests of the mean scores across the identified clusters to infer differentially methylated sites	FDR cut-off < 0.05	0	No	Association
	Smith et al. (2012) [47]	Cord blood and placenta	Site-by-site and global (average measure across all investigated CpGs)	Infinium Human-Methylation27 BeadChip (Illumina)	M values for statistical procedures, β values for data presentation	Linear mixed effects models	FDR cut-off < 0.05 Follow-up in placental tissue: one-sided p-values	4	No	Association
	Schroeder et al. (2012) [43]	Cord blood	Site-by-site	Infinium Human-Methylation27 BeadChip (Illumina)	M values for statistical procedures, β values for data presentation	Linear mixed effects models	FDR cut-off < 0.05	4	No	Association
Combined epigenome-wide and candidate gene studies	Gurnot et al. (2015) [31]	Cord blood	Genome-wide Site-by-site Gene-specific CYP2E1 both individual sites (16 CpGs) and global methylation	Genome-wide Infinium Human-Methylation27 BeadChip array (Illumina) Gene-specific PyroMark Q96 ID pyrosequencer (Qiagen)	Genome-wide β values Gene-specific Methylation percentage	Genome-wide Non-parametric Wilcoxon tests and linear regression models Gene-specific Linear regression models	FDR cut-off < 0.05	0	No	The authors propose an epigenetic mediation mechanism based on three observations, suggesting that the DNAm of CYP2E1 is an epigenetic mechanism to protect the unborn child from the adverse effects of SRIs
	Non et al. (2014) [32]	Cord blood	Genome-wide Site-by-site and regional Gene-specific Site-by-site; 10 candidate genes	Genome-wide Infinium Human-Methylation450 BeadChip (Illumina) Gene-specific Pyromark Q24 Pyrosequencer (Qiagen)	Genome-wide β values Gene-specific Methylation percentage	Multivariate robust standard error regression models	Genome-wide FDR cut-off < 0.1 Gene-specific Bonferroni correction accounting for the number of probes tested in each gene	4	No; the authors do not expect shifts in cell populations to influence much, as neither of the identified genes are important in inflammation or immune system functioning	Association
Candidate gene studies	Galbally et al. (2020) [37]	Placenta and buccal cells (24–72 hours after birth)	Site-by-site; NR3C1 and NR3C2 (13 placental CpGs and 11 buccal CpGs for both genes)	SEQUENOM MassARRAY EpiTYPER platform	Mean methylation percentage of triplicate samples	Univariate ANOVAs	FDR cut-off < 0.25	0	No, but state this as a limitation	Hypothesize that the DNAm of placental NR3C2 CpG 24 mediates the indirect effect of maternal depression on cortisol reactivity at 12 months. In this model, antidepressants may modify the effect of depression on DNAm. To strengthen the hypothesis, the authors performed a mediation analysis based on Hayes et al. [104]. From this analysis, the authors concluded that CpG 24 methylation in the placental NR3C2 reduces the association between maternal depression and infant cortisol reactivity
	Galbally et al. (2018) [51]	Placenta	Site-by-site; OXTR (16 CpGs)	SEQUENOM MassARRAY EpiTYPER platform	Mean methylation percentage of triplicate samples	One-way ANOVAs Significant antidepressant–DNAm associations explored with scatter plots and multiple regression models	p < 0.05	1	No	Association
	McLaughlin et al. (2017) [34]	Buccal cells (24–72 hours after birth)	Regional; ABCB1, CYP2D6, and OPRM1 (mean of all CpG DNAm values within a gene to compare across samples)	Pyromark Q24 Pyrosequencer (Qiagen)	Methylation percentage	One-way ANOVAs	p < 0.05	0	No	Association
	Mansell et al. (2016) [45]	Cord blood (mononuclear cells only)	Site-by-site; NR3C1 promoter (21 CpGs)	SEQUENOM MassARRAY EpiTYPER platform	Log-transformed methylation percentage (mean of triplicate arrays; log base not specified) for regression modelling; mean methylation with no transformation for data presentation	Student’s t-test, ANOVA or pairwise correlation tests (as appropriate) to investigate the associations between exposures, DNAm, and covariates. Multivariate linear regression models to investigate the association of maternal mental well-being and DNAm of NR3C1	Bonferroni corrected p < 0.00079 (accounting for three maternal well-being measures and 21 CpG units) Unadjusted p < 0.05	10	Determine monocyte and lymphocyte frequencies by FACS, used this as a covariate in the model	Association
	Gartstein et al. (2016) [38]	Cord blood	Site-by-site; SLC6A4 promoter (10 CpGs)	PyroMark MD System (Biotage, Qiagen)	Methylation percentages	Identified the most important regions of methylation over the 10 CpGs of SLC6A4 by PCA and literature searches. Using these CpGs, performed hierarchical multiple regression models including interaction terms	p < 0.05^b	2	No	Association
	Ciesielski et al. (2015) [33]	Placenta	Site-by-site; 15 genes (in 27 CpGs)	Infinium Human-Methylation27 BeadChip array (Illumina)	Adjusted^a β values	Logistic regression to determine if the methylation of CpG sites were associated with low birth weight. Fisher’s exact tests and follow-up exact binomial tests were applied to compare the significant CpGs from the regression across the included groups	FDR cut-off < 0.05	1	No	Findings suggest that decreased DNAm of one placental LEPR CpG may be part of a biological mechanism linking ongoing maternal psychiatric disease and poor foetal growth. The authors base this proposal on the observation of associations across their variables and state that they rule out other explanations of this apparent DNAm pattern. Notably, the authors also emphasize that the observed associations do not by default suggest causation and that additional validation and characterization studies are needed
	Soubry et al. (2011) [52]	Cord blood	Regional level; IGF2 DMR (3 CpGs) and H19 DMR (4 CpGs)	Pyromark MD System (Biotage, Qiagen)	Mean methylation percentage from triplicate assays	t-tests between mean methylation over all samples of the DMRs and the covariates Multiple regression models	p < 0.05^b	7 (and SVA)	Validated that IGF2 and H19 DMRs were equal across different cell fractions from cord blood [105]	Association
	Devlin et al. (2010) [44]	Cord blood	Site-by-site; SLC6A (10 CpGs) and BDNF (12 CpGs)	PyroMark MD System (Biotage, Qiagen)	Methylation percentage	MANCOVA models	p < 0.05	0	No, but the authors recognize this limitation	Association
	Oberlander et al. (2008) [36]	Cord blood (mononuclear cells only)	Site-by-site; NR3C1 (13 CpGs) and global (LINE-1)	PyroMark MD System (Biotage, Qiagen)	Methylation percentage	Multiple regression models	p < 0.05	0	No	Findings suggest an association between increased 3^rd trimester maternal depressed mood and higher HPA stress responsiveness in the 3-month-old infant. This association may potentially be mediated epigenetically through DNAm of human NR3C1. However, the authors acknowledge that to infer any functional consequences of neonatal cord blood NR3C1 CpG3 DNAm, they will need more direct evidence, and any mediational relationship cannot be deduced from their study

^†Includes all covariates adjusted for in the statistical models. Most studies also employed other strategies to adjust for or assess suspected covariates; these are presented in the Supplementary Table S3.

*Full explanation of the authors’ reasoning behind a causal relationship is in Supplementary Table S5.

^aWhat is meant by ‘adjusted’ is not specified in the article.

^bAssumed p < 0.05 based on which results are considered significant, but the p was not stated clearly in the article.

(M)AN(C)OVA: (multivariate) analysis of (co)variance; CpG: 5ʹ–Cytosine–phosphate–guanine–3ʹ site; DMR: differentially methylated regions; DNAm: DNA methylation; FACS: fluorescence-activated cell sorting; FDR: false discovery rate; HPA: hypothalamic-pituitary-adrenal; PCA: principal component analysis; (S)SRI: (selective) serotonin reuptake inhibitor; SVA: surrogate variable analysis.

Results of prenatal medication exposure and neonatal DNA methylation

The most examined medication group in prenatal pharmacoepigenetics was antidepressants, investigated in 12 studies of 3,320 neonates (2 EWASs, 8 candidate gene studies, and 2 studies combining an epigenome-wide and a candidate gene approach). In the most recent EWAS, Cardenas et al. [29] discovered 130 differentially methylated CpGs in cord blood samples collected from neonates exposed to antidepressants in utero. One of these sites that mapped to ZNF575 was replicated in an independent cohort [29]. Schroeder et al. [43] found that the exposed neonates had two differentially methylated CpGs in TNFRSF21 and CHRNA4. However, the authors disregarded these findings as false positives considering the small effect sizes (DNAm changes of 1–3%) [43]. In the EWAS conducted by Gurnot et al. [31], three CpGs were differentially methylated in neonates prenatally exposed to serotonin reuptake inhibitors (SRIs; CYP2E1, EVA1, and SLMAP). However, in the EWAS by Non et al. [32], no CpGs were significantly different in neonates exposed to selective serotonin reuptake inhibitors (SSRIs) in utero.

The candidate gene studies investigated CpGs in a total of 32 different genes (Supplementary Table S6). Most of the included genes were chosen based on their suggested association with psychiatric disorders (e.g., the serotonin transporter gene SLC6A4) [32,38,44], stress reactivity (e.g., the glucocorticoid and mineralocorticoid receptor genes NR3C1 and NR3C2) [32,36,37,45], or adverse early life events (e.g., the brain-derived neurotrophic factor gene BDNF) [32,44]. In the studies combining epigenome-wide and candidate gene approaches [31,32], the candidate gene investigation was used to verify the epigenome-wide results. However, except for the verification of CYP2E1 DNAm by Gurnot et al. [31], neither of the significant genes in either of the candidate gene studies were also significant in the EWASs of antidepressants.

Four genes involved in neurotransmitter receptor or transporter activity (NR3C1, SLC6A4, and FKBP5) or neuronal differentiation (BDNF) were investigated across several studies (Table 3). The DNAm of neither NR3C1 nor BDNF was associated with prenatal exposure to antidepressants in any of the studies investigating these genes [32,36,37,45]. For SLC6A4, the results were contradictory. Although Gartstein et al. [38] found an increase in DNAm at six CpGs in cord blood upon prenatal SSRI exposure, Non et al. [32] reported a decrease in DNAm at one CpG in cord blood upon prenatal SSRI exposure when examined by pyrosequencing but not in the epigenome-wide approach. Finally, Devlin et al. [44] found no association between in utero exposure to SSRIs or serotonin and noradrenaline reuptake inhibitors (SNRIs), and DNAm of SLC6A4 in cord blood. A CpG in FKBP5, which encodes a co-regulator of the glucocorticoid receptor, was negatively associated with in utero SSRI exposure in cord blood [32], but not in the placenta [33]. In summary, the results of studies on prenatal antidepressant exposure and DNAm are largely inconsistent.

Table 3. Overview of genes examined in more than one candidate gene study on antidepressants

	NR3C1 Glucocorticoid receptor	SLC6A4 Serotonin transporter	BDNF Brain-derived neurotrophic factor	FKBP5 Regulator of the glucocorticoid receptor
Galbally et al. (2020) [37]	N.S.
Mansell et al. (2016) [45]	N.S.
Gartstein et al. (2016) [38]		+
Ciesielski et al. (2015) [33]				N.S.
Non et al. (2014) [32]	N.S.	–	N.S.	–
Devlin et al. (2010) [44]		N.S.	N.S.
Oberlander et al. (2008) [36]	N.S.
TOTAL STUDIES	4	3	2	2

‘+’: significantly increased DNAm level in medication-exposed group; ‘–’: significantly decreased DNAm level in medication-exposed group; ‘N.S.’: no significant difference between the medication-exposed and non-exposed groups.

Prenatal AED exposure was investigated in two EWASs [46,47], which reported discrepant results. Emes et al. [46] found no global DNAm differences in the cord blood of neonates exposed to AEDs in utero, whereas Smith et al. [47] observed decreased global DNAm in the cord blood of neonates prenatally exposed to AEDs and no global DNAm differences in placental tissue. Furthermore, Emes et al. [46] reported differential DNAm at 662 CpGs when comparing cord blood from neonates that were exposed and not exposed to AEDs in utero, whereas Smith et al. [47] found 14 CpGs with significantly reduced DNAm in the same tissue, including three CpGs that were also significant in placentae.

Two EWASs examined the association between in utero paracetamol exposure and DNAm in placentae [48] or cord blood [35]. Addo et al. [48] reported 24 differentially methylated CpGs in placental tissue when comparing exposed and unexposed pregnancies. Using a different study design, Gervin et al. [35] compared DNAm in long-term paracetamol-exposed children with ADHD to short-term-exposed children with ADHD (2,089 differentially methylated CpGs), unexposed children with ADHD (192 differentially methylated CpGs), and unexposed children without ADHD (6,211 differentially methylated CpGs). Although the studies report vastly different numbers of significant CpGs, both Gervin et al. [35] and Addo et al. [48] concluded that prenatal paracetamol exposure may be associated with DNAm in cord blood from susceptible individuals or placentae, respectively.

Yeung et al. [30] investigated the association between prenatal acetylsalicylic acid exposure and DNAm in cord blood. In the randomized controlled trial, women were randomly assigned to receive 81 mg of acetylsalicylic acid or placebo every day until conception (within six menstrual cycles) and during pregnancy [49]. The DNAm of one CpG (3,500 base pairs upstream of the POU4F1 promoter) in cord blood was significantly associated with prenatal exposure to acetylsalicylic acid [30]. However, Yeung et al. concluded that the association of prenatal acetylsalicylic acid exposure with DNAm in cord blood is negligible, as only one CpG with a minor effect size (1% increase in CpG DNAm) was discovered in their association study [30].

In the EWAS on prenatal methadone exposure, McLaughlin et al. [34] reported a significant increase in buccal cell DNAm of ABCB1, CYP2D6, and OPRM1 in neonates of mothers who were methadone-maintained during pregnancy. The authors argued that their results demonstrated that opioids interact with epigenetic mechanisms, and that the altered DNAm of the opioid metabolism-related genes may have a functional significance that needs further investigation [34].

Discussion

In this review, we have systematically summarized the literature investigating associations between prenatal medication exposure and epigenetic differences in neonates. We included a total of 18 studies on DNAm, examining in utero exposure to antidepressants, AEDs, paracetamol, acetylsalicylic acid, or methadone. We found substantial inconsistency across studies, including heterogeneity in methodology, materials, design, genome coverage, and statistical modelling, making the interpretation of findings and cross-study comparisons challenging. The novelty of the field combining epidemiological and pharmacoepigenetic methods may partly explain this heterogeneity due to a lack of consensus on how to perform analyses and report findings. Therefore, we discuss the results of the reviewed studies with respect to both epidemiological and epigenetic considerations, and suggest 10 recommendations for future studies in prenatal pharmacoepigenetics, as summarized in Box 1.

Box 1. (1)HYPOTHESIS: candidate gene studies should use a plausible hypothesis to guide the study design Hypotheses should be defined prior to designing a candidate gene study, and be guided by principles of teratology, knowledge of pharmacological mechanisms, and epidemiological and biological observations. Hypothesis-free EWASs are also important as the field of prenatal pharmacoepigenetic studies is still emerging.(2)MEDICATION SELECTION: investigate individual medications rather than medication classes Unless the pharmacological and epigenetic mechanisms of action of medications are expected to be similar across the medication class, medications should be analysed on an individual substance level.(3)STATISTICAL POWER: ensure sufficient sample sizes to detect relevant DNAm differences To detect biologically relevant DNAm associations and to ensure valid interpretation of the results, tools developed for power assessments in epigenetic studies should be used when planning such studies.(4)STUDY DESIGN: include a disease comparison group to disentangle medication from indicationStudies should include a disease comparison group to better differentiate the effects of exposure to medication from the underlying maternal disease. This may reduce the impact of confounding by indication.(5)SYSTEMATIC ERROR: assess selection bias, information bias, and confounding Selection bias should be assessed by comparing characteristics of study samples to the target population. The validity of medication exposure, neonatal phenotype, and other covariates should be reported, and information bias and misclassification addressed. Measured confounders of the exposure–outcome association(s) are to be adjusted for and residual confounding investigated. Importantly, cell type heterogeneity should be considered a confounding factor in epigenetic studies.(6)TISSUE SELECTION: biomarkers and extrapolation of DNAm patterns across tissues If the research aim is not only to report a tissue-independent biomarker, but to extrapolate results to other target tissues, the limitations of such translation should be recognized,, and reduced using software applications or data sets on cross-tissue correlations of modifications.(7)LONGITUDINAL PERSPECTIVE: assess persistence of DNAm patterns throughout childhoodThe follow-up of epigenetic patterns later in childhood is essential to assess the relevance of these changes over time, as they may suggest a long-term impact on the phenotypic outcome. (8)DATA INTEGRATION: integrate epigenetic data with complementary omics dataIntegration of complementary omics data, such as genomic and transcriptomic data, can strengthen functional and causal inferences of the findings. (9) CAUSAL INFERENCE: provides a framework for interpreting exposure-outcome associations Causal inference methods, such as two-step Mendelian randomization, may support the inference of causation from exposure–outcome associations, including how medication may impact phenotypic outcome via DNAm changes. Importantly, the underlying assumptions of causal methods are often untestable and, therefore, such methods should be used carefully.(10)REPLICATION: replicate findings using different methods and independent cohorts Replication both across methods and in independent cohorts is essential to increase the validity of the findings and the generalizability of the results to enhance clinical relevance..

Prenatal pharmacoepigenetic candidate gene studies should have a clearly defined hypothesis guided by teratological principles [50] and pharmacological, epidemiological, and biological knowledge (pt. 1, Box 1). Founding the research question on a well-informed hypothesis is fundamental for a transparent and well-designed prenatal pharmacoepigenetic study. This was mostly done in the candidate gene studies included in this review, which provided a rationale for selecting the genes being studied, such as the gene being related to psychiatric disorders (the serotonin transporter gene SLC6A4) [32,38,44] or stress reactivity (the glucocorticoid and mineralocorticoid receptor genes NR3C1 and NR3C2) [32,36,37,45].

Small molecular and structural differences between drugs are known to cause variations in toxicity and teratogenicity [50]. Although the reviewed studies on analgesics focused on one specific medication [30,34,35,48], the studies on psychotropics investigated the effect of medication classes on neonatal DNAm [29,31–33,36–38,43–47,51,52]. In the two studies on AEDs [46,47], several medications were investigated, which may be too broad considering the various different pharmacological [53,54] and epigenetic mechanisms of action of AEDs [21,55,56]. For example, Smith et al. [47] jointly analysed seven different AEDs among 53 women, but also performed a stratified analysis of carbamazepine monotherapy (36 women). In contrast, Emes et al. [46] jointly analysed valproate, lamotrigine, and carbamazepine among nine women, and did not stratify their analyses on individual medications. Analyses on a medication class level may mask effects or give heterogeneous results that are difficult to interpret. Consequently, prenatal pharmacoepigenetic studies should aim to investigate individual medications rather than medication classes and ensure sufficient study power to do this (pt. 2 & pt. 3, Box 1).

The median sample size of the reviewed studies was 201 (IQR of 289), with sample sizes as low as 18 [46] and 23 [31]. None of the studies reported a power assessment to justify the selected sample size. The power of an EWAS depends on many variables, including the significance level, effect size, sample size, array technology, tissue type, and distribution of DNAm differences [57]. Therefore, power calculations are challenging, but simulation studies for power estimation [42,58,59], as well as power assessment tools [42,57], may support the investigation of power in epigenetic studies. In epigenetic epidemiology, the effect sizes are expected to be small, ranging from 0.05 to 0.1 and upwards [60,61], as evidenced in studies on prenatal smoking exposure, with effect sizes commonly ranging from 0.02 to 0.1 [62,63]. We recommend that future prenatal pharmacoepigenetic studies perform and report power assessments in order to ensure sufficient power to detect genuine epigenetic differences between comparison groups (pt. 3, Box 1).

The indication for medication use is an important potential confounder in prenatal pharmacoepigenetic studies, as the observed outcome may be associated with the underlying maternal illness and not the medication used to treat it [64,65]. Among the 14 studies on psychotropics in this review, only eight studies included such a comparison group [29,32,33,36,37,44,51,52]. Notably, seven of these studies found an association between the underlying maternal depression and DNAm in the neonate [29,32,33,36,37,44,51], emphasizing the importance of including this comparison group in future prenatal pharmacoepigenetic studies (pt. 4, Box 1).

When defining medication-exposed comparison groups, more than half of the reviewed studies relied partly or entirely on self-reported medication use during pregnancy [29,35,37,38,43,45–48,51]. This measure does not necessarily reflect the actual medication use [66–68] and is vulnerable to recall bias if reported retrospectively [69]. In five studies, medication exposure was assessed at birth using maternal and/or neonatal blood concentrations of the medication [31,34,37,47,51]. Although informative at birth, this measure does not reflect medication use in earlier stages of pregnancy. Similarly, the eight studies investigating the association between DNAm and neonatal phenotypic outcomes included various outcome definitions without assessing the validity of the measurement [31,33–38,46]. In five studies, data were measured objectively (birth weight [31,33,46] and cortisol levels [36,37]), whereas two studies relied on diagnoses by specialists [34,35], and one study used parent reports on infant temperament [38]. An assessment of medication exposure, neonatal phenotype, and covariates is crucial to avoid misclassification. Therefore, we recommend that future prenatal pharmacoepigenetic studies perform sensitivity analyses to assess the robustness of the findings, taking into account the validity of the measures (pt. 5, Box 1). For example, methods to quantify the impact of exposure and outcome misclassifications, such as probabilistic bias analysis [70], are highly recommended [64].

A wide range of different covariates were considered in the studies (Supplementary Table S3). Selecting an appropriate set of confounders to control for is critical to avoid systematic bias (pt. 5, Box 1). Ten studies selected confounders by assessing associations between covariates, exposure, and outcome [31,33,34,36–38,45,46,51,52]. We suggest to control for covariates that are assumed to be confounders, i.e., covariates that are not part of the causal path, and that are both a cause of the exposure and the outcome. Therefore, the covariates to be accounted for should occur upstream of the prenatal exposure, while mediators (which are part of the causal pathway) should not be accounted for when investigating the total effect of the exposure on the outcome. The specific covariates to be taken into account need to be assessed for each individual study, as the relevance of the covariates depends on several factors (e.g., study design and tissue type), and include both technical covariates related to laboratory procedures and biological covariates. Examples of biological covariates to be evaluated are maternal age, smoking during pregnancy, infant sex, gestational age, and folate use in pregnancy, which were the covariates most commonly accounted for in the studies included in this review (Supplementary Table S3). To this end, future prenatal pharmacoepigenetic studies may also benefit from implementing causal inference tools, such as directed acyclic graphs (DAGs) [71], to identify a sufficient set of confounders for adjustments [64]. Such investigations can be complemented by assessing whether the selected confounders largely capture the model variability, as was performed by surrogate variable analysis in Gervin et al. [35], and by principal component analyses in Addo et al. [48] and Cardenas et al. [29]. For an excellent overview of a general approach to identify relevant confounders in observational studies, please refer to the review by VanderWeele [72].

The majority of the reviewed studies were based on cord blood, which consists of cells exhibiting cell type-specific DNAm patterns [40]. Therefore, prenatal pharmacoepigenetic studies should consider whether DNAm differences associated with medication exposure reflect variation in constituent cell types, which are known to mediate or confound the exposure associations [40,73]. When investigating cell-type proportions as a mediator, in order to assess the direct effect of medication exposure on associated DNAm differences, it may be necessary to adjust for estimated or measured cell-type composition (see, e.g., Liu et al. [74] and Gervin et al. [75]). However, if the total effect of medication exposure on DNAm is more interesting, for instance when searching for potential biomarkers of a phenotypic outcome, the cell-type composition should not be accounted for, as it may remove relevant DNAm–phenotypic outcome associations (see, e.g., Ollikainen et al. [76]). Cell type composition may act as a confounder when assessing the extent to which DNAm mediate the effect of drug exposure on a phenotypic outcome, and should in such instances be accounted for [73]. In summary, variation in cell-type composition confers an important covariate in epigenetic studies, and should be appropriately evaluated [40,60,73,77]. Surprisingly, only six of the reviewed studies considered cell-type composition in their analyses [29,30,35,45,48,52], emphasizing the need for increased awareness among prenatal pharmacoepigenetic researchers to evaluate cell types in future studies (pt. 5, Box 1). There are several different methods to determine and account for the cell-type composition in tissue samples, and these are extensively described in the excellent recent review by Teschendorff and Relton [40].

Among the reviewed records, eight studies hypothesized an association between DNAm and neonate phenotypes [31,33–38,46], including five studies concerning brain-related phenotypic outcomes [34–38]. However, DNAm in peripheral surrogate tissues does not necessarily resemble DNAm in the target tissue [78–80], which challenges the accuracy of the extrapolation of the findings. Although this limitation was acknowledged in most studies [29–32,35–38,43–45,47,48,51], only one study attempted to reduce the constraint by including a correlation analysis of select CpGs across adult whole blood and brain tissue [29]. Importantly, investigation of peripheral tissues is still considered valuable, since biomarkers of maternal disease or child developmental outcomes do not need to be from the relevant tissue (i.e., do not need to be tissue-specific). However, to realize the ultimate aim of prenatal pharmacoepigenetics of gaining direct mechanistic insights into how medication exposure impacts the foetus with potential phenotypic consequences, future studies should validate tissue extrapolation by, for example, investigating cross-tissue correlations in available databases [61,81–85] (pt. 6, Box 1). Yet, current databases are mostly available on adult tissues, limiting the relevance to prenatal pharmacoepigenetic studies. Researchers have been calling for initiatives to develop biobanks of foetal and child brain specimens, while also taking into account the ethical issues of building such biobanks [86].

Only one of the reviewed studies investigated DNAm patterns longitudinally during childhood, finding that DNAm at a CpG in ZNF575 persisted into early childhood [29]. Though investigating the DNAm at birth provides information on the immediate impact of prenatal medication exposure, the follow-up of epigenetic patterns later in childhood is valuable to assess the persistence over time and increase the clinical relevance of the findings (pt. 7, Box 1).

The clinical relevance of prenatal pharmacoepigenetic research may also be strengthened by functional and causal interpretations of the results. Using a multi-omics approach with integration of omics data (e.g., genomics, epigenomics, and transcriptomics data) [87] could substantiate the epigenetic findings (pt. 8, Box 1). However, the reviewed studies only used single omics data (i.e., epigenomics). Although single omics data are potentially useful both as biomarkers and in providing insight into biological pathways, this is limited to correlations or associations often reflecting reactive, rather than causative, processes. We recommend that future studies include additional omics data, as this may enable (i) investigation of the functional consequences of DNAm on gene expression [40], (ii) adjustment for the genetic variation associated with DNAm variation [88–90], and (iii) utilization of genomic methylation quantitative trait loci (mQTLs) to implement causal inference methods, such as two-step Mendelian randomization [40,91] (pt. 9, Box 1). There are several openly accessible resources making omics data available for integration, as thoroughly reviewed by Walton, Relton and Caramaschi [92].

Causal modelling and reasoning are increasingly being applied in genetic epidemiology to strengthen the ability to make causal inferences about associations, but it is still new to the field of pharmacoepigenetics [40,64]. For example, two-step Mendelian randomization [91,93–95] has been used to assess how DNAm can mediate an association between prenatal exposure and phenotypic outcomes in children [96,97]. Notably, only one reviewed study attempted to make causal inferences about an association between prenatal antidepressant exposure, cord blood DNA, and infant stress reactivity in a mediation analysis [37]. We foresee important advances in future prenatal pharmacoepigenetic studies using the causal inference framework (pt. 9, Box 1). Importantly, the causal models rely on assumptions that need to be met for them to be valid [98]. As these assumptions are often untestable [98], careful use of the causal inference framework in pharmacoepigenetic studies is essential.

Lastly, to validate findings, replication using a different technology and in an independent cohort is essential, both to determine the robustness of the associations and to assess the level of technical and biological variation. Notably, only two of the reviewed studies applied more than one method to assess DNAm [31,32], and only one study attempted to validate their results in an independent cohort [29], emphasizing the need for an increased focus on replication in prenatal pharmacoepigenetics (pt. 10, Box 1). Several multi-cohort consortia to enable replication of studies are already in place, such as the Pregnancy and Childhood Epigenetics consortium (PACE) [99].

Conclusion

Investigating the potential effects of pharmacological treatment in pregnancy is essential to establish foetal epigenetic safety, understand the underlying mechanisms, and recognize the clinical consequences for the offspring. However, studies on prenatal medication exposure and epigenetic changes are largely heterogeneous and inconsistent. To improve the quality, comparability, and interpretability of future prenatal pharmacoepigenetic studies, we propose 10 recommendations bridging the fields of prenatal epidemiology and epigenetics. Epidemiological approaches and causal inference frameworks will reduce systematic bias and improve our ability to interpret exposure–outcome associations, including how medications may impact phenotypic outcomes via changes in DNAm. Furthermore, it is essential to consider the persistence of DNAm patterns over time and the potential for cross-tissue extrapolation when assessing the biological relevance of the epigenetic contribution. Importantly, integrating more omics data and implementing two-step Mendelian randomization can strengthen the functional and causal inferences of the findings. In conclusion, a consensus on how to perform and report prenatal pharmacoepigenetic studies will fuel the development of the field and contribute to future high-quality studies of clinical relevance.

List of Abbreviations

ADHD	=	Attention-deficit/hyperactivity disorder
AED	=	Antiepileptic drug
ATC	=	Anatomical Therapeutic Chemical
CpG	=	5ʹ–Cytosine–phosphate–guanine–3ʹ site
DAG	=	Directed acyclic graph
DOHaD	=	Developmental Origins of Health and Disease
DNAm	=	DNA methylation
EWAS	=	Epigenome-wide association study
FDR	=	False discovery rate
IQR	=	Interquartile range
mQTL	=	Methylation quantitative trait locus
NSAID	=	Non-steroidal anti-inflammatory drug
SNRI	=	Serotonin and noradrenaline reuptake inhibitor
SRI	=	Serotonin reuptake inhibitor
SSRI	=	Selective serotonin reuptake inhibito

Author’s contribution

EWO, HMEN, and KG conceived the idea of the systematic literature review. EWO, HMEN, and KG planned the searches. EWO performed the searches. EWO and KG performed the screening and selection of records. EWO extracted the data and drafted the first version of the paper. EWO, HMEN, and KG all revised the paper. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Acknowledgments

We thank Christian Page and Mollie Wood for their critical appraisal of the article. We also thank the librarians at the Library of Medicine and Science of the University of Oslo, particularly Bente Kathrine Rasch and Hilde Strømme, for excellent support in developing the search strategy for the review.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the PharmaTox Strategic Research Initiative, Faculty of Mathematics and Natural Sciences, University of Oslo, Oslo, Norway and the European Research Council (ERC) Starting Grant “DrugsInPregnancy” [grant number: 678033].