Effect of Prenatal Opioid Exposure on the Human Placental Methylome

Prenatal exposure to addictive drugs can lead to placental epigenetic modifications, but a methylome-wide evaluation of placental DNA methylation changes after prenatal opioid exposure has not yet been performed. Placental tissue samples were collected at delivery from 19 opioid-exposed and 20 unexposed control full-term pregnancies. Placental DNA methylomes were profiled using the Illumina Infinium HumanMethylationEPIC BeadChip. Differentially methylated CpG sites associated with opioid exposure were identified with a linear model using the ‘limma’ R package. To identify differentially methylated regions (DMRs) spanning multiple CpG sites, the ‘DMRcate’ R package was used. The functions of genes mapped by differentially methylated CpG sites and DMRs were further annotated using Enrichr. Differentially methylated CpGs (n = 684, unadjusted p < 0.005 and |∆β| ≥ 0.05) were mapped to 258 genes (including PLD1, MGAM, and ALCS2). Differentially methylated regions (n = 199) were located in 174 genes (including KCNMA1). Enrichment analysis of the top differentially methylated CpG sites and regions indicated disrupted epigenetic regulation of genes involved in synaptic structure, chemical synaptic transmission, and nervous system development. Our findings imply that placental epigenetic changes due to prenatal opioid exposure could result in placental dysfunction, leading to abnormal fetal brain development and the symptoms of opioid withdrawal in neonates.


Introduction
Opioid use disorder (OUD) in pregnancy has increased four-fold over the past decade in the United States (US) and now impacts 1-2% of all pregnancies in the US [1,2]. The primary neonatal outcome, Neonatal Opioid Withdrawal Syndrome (NOWS), typically presents with signs and symptoms of opioid withdrawal after in utero exposure, with often prolonged courses of opioids and lengthened neonatal hospitalizations [3,4]. Prenatal exposure to opioids is known to be associated with adverse pregnancy outcomes, including fetal growth restriction, risk for preterm birth, and higher rates of neonatal intensive care unit (NICU) admission [5,6]. In addition, there is an association between prenatal opioid exposure and risk for altered neuronal function and behavior throughout childhood as well as risk for neurodevelopmental impairment [7][8][9].
The placenta is the critical physiological interface connecting opioid use in pregnant mothers to prenatal exposure of the fetus. In limited cohort studies, maternal OUD has been associated with structural changes within the placental villi, including delayed villous maturation that is associated with stillbirths [10]. Several studies have additionally shown associations between sequestration of toxic substances that inhibit nutrient transport and placental lesions indicating functional insufficiency that includes poor invasion, vascular changes, and vasoconstriction [11,12]. The placenta is also a key target for epigenetic modifications induced by exogenous substances and is the master regulator of the fetal environment [13]. As a uniquely sensitive organ to environmental influences, the placenta is an established site of epigenetic modifications with subsequent differences in gene expression and tissue differentiation [14,15]. Previous studies have shown an association between placental DNA methylation in key genes after exposure to maternal stress and differences in infant neurobehavioral profiles in the first few months after birth [16,17]. These neurobehavioral profiles, as measured by the NICU Network Neurobehavioral Scale (NNNS), are known to be associated with later neurodevelopmental outcomes [16,17]. Though previous studies have demonstrated placental epigenetic modifications induced by exposure to other addictive substances or maternal psychosocial stress, the effect of prenatal opioid exposure in humans on genome-wide epigenomic modifications in the placenta has yet to be explored [17][18][19][20].
The link between maternal OUD and epigenetic modification in other tissues, such as maternal blood and maternal/infant saliva, is well established [21,22]. Opioid exposure in non-pregnant adults is known to lead to changes in DNA methylation in the mu-opioid receptor gene (OPRM1; the primary molecular target for the analgesic and addictive properties of opioids), and it has been associated with DNA methylation changes at a genome-wide level [23,24]. In recent years, researchers have examined changes in DNA methylation in OPRM1 in infants with NOWS as predictive of NOWS severity [25,26]. Specifically, associations between increased cytosine:guanine (CpG) dinucleotide methylation within the OPRM1 promoter region and NOWS severity have been observed in maternal and infant saliva samples, with increased methylation at select CpG sites associated with more severe NOWS with higher pharmacologic treatment rates, or higher rates of requiring a multi-drug regimen for NOWS treatment [25,26]. Our previous pilot study probed methylation in placental tissues at six CpG sites within the OPRM1 promoter. Methylation at these sites was not found to be associated with NOWS severity, nor was it associated with opioid exposure [27].
It is well known that there is a high density of hypomethylated CpG sites (or CpG islands) in the promoter region of many genes (including OPRM1), and the methylation of promoter CpG sites can inhibit gene transcription. Moreover, altered DNA methylation in the gene body (usually hypermethylated to avoid generating an alternative transcription start site) and the enhancer region (can be up to 1 Mbp away from the gene) can also influence the transcription of genes. Therefore, the present study took a broader and non-biased approach to detect methylome-wide alterations in placental DNA methylation following prenatal exposure to opioids. To our knowledge, genome-wide DNA methylation analyses have not been reported in placental tissues from opioid-dependent women. Identification of genes that contain differentially methylated CpG sites in opioid-exposed placental tissues can facilitate the narrowing of genes associated with NOWS risk and developmental outcomes. Given the associations between epigenetic modifications and differences in infant neurobehavioral profiles following prenatal exposure to drugs of abuse, as well as the findings from our previous studies showing a relationship between NOWS and epigenetic modifications in OPRM1, we sought to compare genome-wide DNA methylation in placental tissues of opioid-exposed versus opioid-unexposed pregnancies.

Setting
Boston Medical Center (BMC) is the largest urban safety-net hospital in New England, with a specialized prenatal clinic for women with OUD. Treatment options for pregnant women with OUD include methadone and buprenorphine. BMC practices a rooming-in model of care where infants room-in with their mothers in the postpartum room until maternal discharge, and then the infants are transferred to the pediatric inpatient unit for continued NOWS monitoring and treatment, where their mothers can continue to room-in. The Eat, Sleep, and Console (ESC) method is utilized for NOWS assessments, with an "as needed" methadone pharmacologic treatment protocol.

Subjects
This study was approved by the Boston University Medical Campus Institutional Review Board. Placentas were collected from human subjects between July 2019 and July 2020. Eligibility criteria for inclusion in the opioid cohort included confirmed maternal OUD with opioid exposure for at least 30 days prior to delivery based on chart review and toxicology screening results, singleton pregnancies, gestational age (GA) of ≥36 weeks, and delivery at BMC. For the control group, eligibility criteria for placental collection included delivery at BMC, gestational age (GA) ≥ 36 weeks, and absence of a known substance use disorder per the electronic medical record (EMR) problem list, admission note, and review of toxicology screen lab results. Controls were matched based on month of delivery with the opioid group. Informed consent was waived for this study due to the collection of discarded placental tissue only for these subjects, along with limited de-identified basic demographic data from the EMR. The EMRs for these subjects were accessed once at the time of placental collection, with no identifiers or master code collected per IRB guidelines. As shown in Table 1, 40 pregnant women (20 cases and 20 controls) were recruited for this study. They were all non-Hispanic. Except for two black women (1 case and 1 control), all others were white women (19 cases and 19 controls).

Phenotype Data Collection
Limited data points were collected at the time of placental collection, including maternal race and ethnicity, maternal age, smoking status, gestational age at delivery, infant birth weight, and infant sex. All data were hand abstracted and input into an electronic database. Data were checked for accuracy and missingness, with any discrepancies addressed prior to data analysis.

DNA Extraction
Genomic DNA was isolated from the placental tissue per standard protocols in batches every 6 months. Genomic DNA extraction was performed by the Boston University Genomics Core lab using the DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany), following manufacturer instructions. DNA concentration and purity were quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and DNA samples were stored at −20 • C until the Illumina Infinium MethylationEPIC array assay was performed.

Illumina EPIC DNA Methylation Array Assay and Raw Data Processing
Extracted genomic DNA was sent to the Yale Center for Genome Analysis (YCGA) for DNA methylome analysis using the Illumina Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA, USA). The GenomeStudio software (Illumina) was used to generate β values for each CpG site. Data quality control (QC) and normalization, as well as statistical analysis, were performed as in a recent study that investigated peripheral blood DNA methylomic changes in European-American women with OUD using the Illumina Infinium MethylationEPIC BeadChip [23]. Briefly, QC and normalization were performed using the 'minfi' R package (v.1.36.0) [28]. Samples were excluded if more than 1% of the probes had a detection p-value > 0.01, and probes were excluded if more than half of samples had a detection p-value > 0.01. Moreover, cross-reactive probes and probes overlapping genetic variants at targeted CpG sites or single base extension sites, as well as probes with genetic variants overlapping the body of the probe, were excluded [29]. Additionally, probes mapped to X and Y chromosomes were excluded prior to analysis. The 'ComBat' method in the 'sva' R package (v.3.38.0) [30] was applied to correct for batch effects. Stratified quantile normalization of beta and M values was conducted using the 'preprocessQuantile' function in the 'minfi' R package (v.1.36.0). This method stratified probes by genomic region (CpG islands, shores, and shelves) to account for differences in region-specific methylation distributions. Density plots were generated to evaluate the distribution of beta (β) values before and after quantile normalization.

Statistical Analysis
All statistical analyses were performed within R v.4.0.3 (www.r-project.org, accessed on 1 March 2022). Baseline demographics for the opioid and control cohorts were summarized. M-values of normalized probe intensity were used for differential methylation analysis, and beta values were used for all data visualization. Differentially methylated probes (DMPs) associated with opioid exposure were identified with a linear model using the 'limma' R package (v.3.46.0) [31]. The model included coefficients to account for variance attributable to confounding factors (infant sex, infant birth weight, and batch). In limma, the design matrix was set up as:~opioid status + infant sex + infant birth weight + batch. The effect of the 'opioid status' coefficient was extracted to determine DMPs. We employed the recommended significance threshold for DNA methylation association analyses of p < 9.4 × 10 −8 as determined from permutation analysis [15]. Methylation of the top differentially methylated CpG sites were visualized via volcano plots [32], heatmaps with hierarchical clustering [33], and site-specific plots of normalized beta values by opioid-exposure group generated with R scripts. To identify differentially methylation regions (DMRs) spanning multiple CpG sites, we utilized the 'DMRcate' R package (v.2.4.1) [34].

Subjects and Demographic Data
Placental samples (20 control and 20 opioid-exposed) were collected and processed for DNA methylation analysis. One sample from the opioid-exposed group was omitted from the analysis following QC inspection due to a poor signal intensity. Twenty control and 19 opioid-exposed samples were included in the analysis, and demographic information for these samples is shown in Table 1. Opioid-exposed subjects were significantly younger than controls (p = 0.010). There were no significant differences in gestational age (GA) at delivery (p = 0.220); however, neonates from opioid-exposed moms showed a near-significant trend for lower birthweights compared to those of controls (p = 0.050). Notably, five subjects (three opioid-exposed and two controls) had missing data for maternal age, gestational age at delivery, and delivery type.

Differential Methylation in Opioid-Exposed Placentas
Post-filtering, a total of 710,952 CpG sites (from 865,918 CpG sites in total) were used for differential methylation analysis. Density plots of raw and quantile-normalized beta values in opioid-exposed and control groups are shown in Figure 1. There were no individual CpG probes/sites that passed the FDR threshold (p < 9.4 × 10 −8 ) for significance, with a minimum unadjusted p-value of 7.9 × 10 −6 for cg05771324, which is mapped to the Phospholipase D1 gene (PLD1). There were 1,958 CpG sites with an unadjusted p-value ≤ 0.005. Of these sites, 684 (~35%) had |∆β| ≥ 0.05 where ∆β = β (opioid-exposed) − β control . A volcano plot of ∆β values by −Log 10 p-value (unadjusted) is shown in Figure 2, and the 684 sites with both p-value ≤ 0.005 and |∆β| ≥ 0.05 are plotted in red. Of these sites, 397 (~58%) were hypomethylated, and 287 (~42%) were hypermethylated in opioidexposed samples versus controls. A heatmap displaying normalized beta values indicative of methylation levels at all 684 sites is shown in Figure 3. CpG sites (rows) were clustered by hierarchical sorting, and those with similar trends in methylation were grouped regardless of opioid exposure.
We identified 199 opioid exposure-associated DMRs that had a minimum smoothed FDR < 0.05. Of these, 94 (~47%) were hypomethylated, and 105 (~53%) were hypermethylated in opioid-exposed samples relative to controls. The top 15 DMRs (FDR < 0.05) as- Normalized beta values are shown by groups (controls vs. opioid-exposed) for the top 6 annotated CpG sites associated with opioid exposure. CpG site names are above each plot, and chromosome numbers, associated genes, genomic regions, and unadjusted p-value are below each plot.
We identified 199 opioid exposure-associated DMRs that had a minimum smoothed FDR < 0.05. Of these, 94 (~47%) were hypomethylated, and 105 (~53%) were hypermethylated in opioid-exposed samples relative to controls. The top 15 DMRs (FDR < 0.05) associated with opioid exposure are shown in Table 2. One DMR on Chromosome 2, spanning 16 CpGs (1,322 bp), was mapped to the first exon of the Boule Homolog, RNA Binding Protein gene (BOLL; Figure 5A), which is involved in meiosis and gamete development [41]. In opioid-exposed samples, methylation levels were 9.0% higher on average at CpG sites within this region relative to control samples (min smoothed FDR = 8.5 × 10 −9 ). Another DMR (1,120 bp) on Chromosome 10 spanning 10 CpGs was mapped to the first exon of the Potassium Calcium-Activated Channel Subfamily M Alpha 1 gene (KCNMA1; Figure 5B), which has been associated with substance use disorders [42,43]. Opioid exposed samples had methylation levels that were an average of 11.8% lower than control samples across this region (min smoothed FDR = 2.6 × 10 −5 ).  In opioid-exposed samples, methylation levels were 9.02% higher on average at CpG sites within this region relative to control samples (min smoothed FDR = 8.5 × 10 −9 ). (B) A DMR (1,120 bp) on chromosome 10 spanning 10 CpG probes mapped to the first exon of the Potassium Calcium-Activated Channel Subfamily M Alpha 1 gene (KCNMA1). Opioid exposed samples had methylation levels that were an average of 11.8% lower than control samples across this region (min smoothed FDR = 2.6 × 10 −5 ). Mean normalized beta values across the DMR are shown for each group (control vs. opioid exposed) at the bottom of each panel. In opioid-exposed samples, methylation levels were 9.02% higher on average at CpG sites within this region relative to control samples (min smoothed FDR = 8.5 × 10 −9 ). (B) A DMR (1120 bp) on chromosome 10 spanning 10 CpG probes mapped to the first exon of the Potassium Calcium-Activated Channel Subfamily M Alpha 1 gene (KCNMA1). Opioid exposed samples had methylation levels that were an average of 11.8% lower than control samples across this region (min smoothed FDR = 2.6 × 10 −5 ). Mean normalized beta values across the DMR are shown for each group (control vs. opioid exposed) at the bottom of each panel.

Functional Enrichment
The 684 DMPs with both p ≤ 0.005 and |∆β| ≥ 0.05 were mapped to 258 annotated genes that were then used as input for enrichment analysis. The top 15 Gene Ontology (GO) enrichment terms (all unadjusted p ≤ 6.4 × 10 −4 and all adjusted p ≤ 0.054) are shown in Figure 6A Figure 6B. None of the KEGG enrichment terms reached the adjusted significance threshold.

Functional Enrichment
The 684 DMPs with both p ≤ 0.005 and |∆β| ≥ 0.05 were mapped to 258 annotated genes that were then used as input for enrichment analysis. The top 15 Gene Ontology (GO) enrichment terms (all unadjusted p ≤ 6.4 × 10 −4 and all adjusted p ≤ 0.054) are shown in Figure 6A. Ten of the top fifteen GO enrichment terms met the adjusted threshold for significance, including regulation of synapse assembly (GO:0051963; adjusted p = 1.  Figure 6B. None of the KEGG enrichment terms reached the adjusted significance threshold.   We additionally performed enrichment analysis for the DMRs. A total of 199 DMRs (minimum smoothed FDR < 0.05) were mapped to 174 genes that were used as input. Of the top 15 GO enrichment terms and the top 9 KEGG enrichment terms (Figure 7), only one reached the adjusted threshold for significance: the GO term 'integral component of plasma membrane' (GO:0005887; adjusted p = 3.7 × 10 −4 ; 31 genes). This term was also implicated in the DMP enrichment analysis. None of the KEGG enrichment terms reached the adjusted significance threshold. We additionally performed enrichment analysis for the DMRs. A total of 199 DMRs (minimum smoothed FDR < 0.05) were mapped to 174 genes that were used as input. Of the top 15 GO enrichment terms and the top 9 KEGG enrichment terms (Figure 7), only one reached the adjusted threshold for significance: the GO term 'integral component of plasma membrane' (GO:0005887; adjusted p = 3.7 × 10 −4 ; 31 genes). This term was also implicated in the DMP enrichment analysis. None of the KEGG enrichment terms reached the adjusted significance threshold.

Discussion
This is the first genome-wide study of placental DNA methylation alterations that were associated with opioid exposure during pregnancy. Although no CpG sites reached the adjusted significance threshold (p < 9.4 × 10 −8 ), a subset of 684 CpG sites and affiliated genes with methylation influenced by in utero opioid exposure were identified using nominal thresholds (p ≤ 0.005 and |∆β| ≥ 0.05). Several of the top DMPs were annotated to genes of potential interest, including PLD1, a gene implicated in a variety of processes, including receptor endocytosis [44] and cytoskeletal organization [45]. The other isoform of this gene, PLD2, is required for µ-opioid receptor endocytosis and opioid desensitization [46]. The effects of aberrant PLD1 regulation in the placenta are not clear, but as this gene is involved in a variety of processes for cellular homeostasis, this alteration could have a significant impact on placental function with subsequent effects on fetal development. We additionally identified a DMP in MGAM, a membrane digestive enzyme involved in starch digestion and linked to chronic diarrhea in children [39]. MGAM has been most thoroughly studied in the mucosa [47] and intestinal tissue [48], but considering the prevalence of feeding problems and low birth weight in opioid-exposed neonates [6,49], it is possible that dysregulation of this gene impacts metabolic processes in the placenta, impacting nutrient transfer to the developing fetus. One of the 684 DMPs used in the enrichment analysis mapped to the first exon of OPRM1 (cg15085086; p = 0.004). In opioid-exposed samples, this site had 5.54% less methylation than controls. As we previously established an association between OPRM1 and NOWS outcomes [25], it is likely that a larger study would be more successful in the identification of additional CpG sites annotated to OPRM1 in association with opioid exposure.
We additionally identified several interesting genes annotated to DMRs that were associated with opioid exposure, including KCNMA1 ( Figure 5B). The KCNMA1 subunit is part of a large-conductance potassium channel that affects neuronal excitability and has been implicated in opioid analgesia [50], neuropathic pain [51], and kappa opioid receptor-mediated cardiomyocyte response to ischemia [52]. It has also been implicated in other substance use disorders (SUDs), including alcohol use disorder (AUD) [53,54] and methamphetamine use disorder (MUD) [55]. Additionally, KCNMA1 has been found to have decreased gene expression in endometrial stromal fibroblasts during decidualization, supporting its function as a regulator of embryonic development [56]. The exact role of KCNMA1 in the placenta, particularly in the context of opioid exposure, remains unclear. Future studies can better ascertain whether this is a viable target gene in the treatment and study of NOWS.
The results of our enrichment analyses for the top DMPs and DMRs indicate a disruption in gene pathways related to synaptic assembly and organization, chemical signal transduction, neuron differentiation, and nervous system development. There is substantial evidence from rodent models that prenatal opioid exposure influences neuron differentiation, synaptogenesis, and postsynaptic density protein expression in the central nervous system (CNS) [57][58][59]. While it is not feasible to analyze fetal or infant brain samples, the implication of these pathways in our placental enrichment analyses may be somewhat predictive of epigenetic adaptations present in the CNS [60].
This study has several limitations. First, it is limited by the small sample size that was underpowered to detect differences in DNA methylation at individual CpG sites that passed the strict FDR threshold for significance. Moreover, given that these were samples from subjects who did not sign informed consent, we were limited by few clinical variables (as well as missing data) to correlate with the samples. Future studies should include more detailed clinical and drug exposure information through a consented study-specific questionnaire. Second, accounting for additional variables in the linear regression, including maternal age, gestational age, and smoking status, could significantly improve the ability to identify significant associations between opioid exposure and placental DNA methylation. Third, there was also variability in opioid exposure type, which primarily included methadone (~53%), buprenorphine (~16%), and fentanyl (26%). Methadone and buprenorphine both have much longer elimination half-lives (~30 h) [61,62] versus that of fentanyl (3-7 h, i.v.) [63]. Future studies of larger sample sizes would better serve to determine the effects of individual opioid medications on placental DNA methylation. Fourth, beyond the epigenetic modifications induced by exposure to opioids themselves, there are likely additional implications related to the presentation of spontaneous opioid withdrawal during pregnancy. Future studies could better disentangle the effects of opioid agonism versus opioid withdrawal on placental DNA methylation and subsequent NOWS outcomes. Fifth, we need to further explore the function of those identified hypermethylated or hypomethylated CpG sites to clarify whether altered methylation at these CpG sites could influence their gene transcription. Finally, the villous placental tissue samples evaluated for this analysis contain a heterogeneous mix of cell populations, so our overall results do not give exact information on cell-specific DNA methylation changes. Ongoing studies evaluating gene methylation differences within isolated key placental cell populations will provide further insight into how gene expression changes can be correlated with altered placental function and fetal development.

Conclusions
This pilot work is the first step to the further examination of epigenetic variation in placental tissue within opioid-exposed pregnancies, with potential importance for improving neonatal outcomes and informing transgenerational risk. Future studies could include additional assays of epigenetic modifications (such as ATAC-seq and quantification of histone acetylation) for a more comprehensive understanding of how opioids influence epigenetic profiles and chromatin accessibility in the placenta. Furthermore, the inclusion of high-throughput RNA-seq would enable a concomitant analysis of transcriptomic and epigenetic adaptations and help facilitate the identification of genomic regions underlying variation in behavioral features of NOWS. Genes determined to be associated with NOWS severity and/or developmental outcomes could be potential targets for therapeutics [64][65][66]. Informed Consent Statement: Informed consent was waived for this study due to the collection of discarded placental tissues only for these subjects, along with limited de-identified basic demographic data from the electronic medical record (EMR).

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.