DNA methylation in cocaine use disorder–An epigenome-wide approach in the human prefrontal cortex

Background Cocaine use disorder (CUD) is characterized by a loss of control over cocaine intake and is associated with structural, functional, and molecular alterations in the human brain. At the molecular level, epigenetic alterations are hypothesized to contribute to the higher-level functional and structural brain changes observed in CUD. Most evidence of cocaine-associated epigenetic changes comes from animal studies while only a few studies have been performed using human tissue. Methods We investigated epigenome-wide DNA methylation (DNAm) signatures of CUD in human post-mortem brain tissue of Brodmann area 9 (BA9). A total of N = 42 BA9 brain samples were obtained from N = 21 individuals with CUD and N = 21 individuals without a CUD diagnosis. We performed an epigenome-wide association study (EWAS) and analyzed CUD-associated differentially methylated regions (DMRs). To assess the functional role of CUD-associated differential methylation, we performed Gene Ontology (GO) enrichment analyses and characterized co-methylation networks using a weighted correlation network analysis. We further investigated epigenetic age in CUD using epigenetic clocks for the assessment of biological age. Results While no cytosine-phosphate-guanine (CpG) site was associated with CUD at epigenome-wide significance in BA9, we detected a total of 20 CUD-associated DMRs. After annotation of DMRs to genes, we identified Neuropeptide FF Receptor 2 (NPFFR2) and Kalirin RhoGEF Kinase (KALRN) for which a previous role in the behavioral response to cocaine in rodents is known. Three of the four identified CUD-associated co-methylation modules were functionally related to neurotransmission and neuroplasticity. Protein-protein interaction (PPI) networks derived from module hub genes revealed several addiction-related genes as highly connected nodes such as Calcium Voltage-Gated Channel Subunit Alpha1 C (CACNA1C), Nuclear Receptor Subfamily 3 Group C Member 1 (NR3C1), and Jun Proto-Oncogene, AP-1 Transcription Factor Subunit (JUN). In BA9, we observed a trend toward epigenetic age acceleration (EAA) in individuals with CUD remaining stable even after adjustment for covariates. Conclusion Results from our study highlight that CUD is associated with epigenome-wide differences in DNAm levels in BA9 particularly related to synaptic signaling and neuroplasticity. This supports findings from previous studies that report on the strong impact of cocaine on neurocircuits in the human prefrontal cortex (PFC). Further studies are needed to follow up on the role of epigenetic alterations in CUD focusing on the integration of epigenetic signatures with transcriptomic and proteomic data.


Introduction
Cocaine is one of the most frequently consumed psychostimulants in the world with around 21 million individuals having used cocaine in 2020 (1). During the last decade, a stable increase in cocaine consumption and cocaine-associated hospitalization rates has been observed (1,2). While a substantial fraction of individuals retain control over cocaine intake, up to 21% of regular cocaine users develop cocaine dependence characterized by repeated cycles of binge intoxication, withdrawal, and compulsive cocaine-seeking (3). Chronic use of cocaine, especially crack cocaine, increases the risk for the development of cardiovascular, cerebrovascular, and respiratory diseases (4). Also, it is associated with substance abuse of other drugs such as alcohol, cannabis, and opioids (5). Thus, chronic cocaine consumption as observed in cocaine use disorder (CUD) enhances the overall disease burden, at the individual level significantly reduces life quality, and thereby contributes to an increased risk for premature death. Investigating the pathophysiology of CUD enables the characterization of the underlying disease processes and at the same time provides a basis for the development of novel therapeutic approaches.
Substance use disorders (SUDs) including CUD are understood as brain disorders with profound alterations at the structural, functional, and molecular levels in the human brain (6)(7)(8). Druginduced effects on epigenetics are hypothesized to contribute to structural and functional changes in the brain via the establishment of altered transcriptional programs [for review see (9)(10)(11)]. Thus, profiling of epigenetic alterations contributes to a better understanding of molecular mechanisms underlying the neuroplastic changes associated with CUD. DNA methylation (DNAm) represents the most prominent epigenetic mechanism involved in gene expression regulation and DNAm changes in the brain following acute and chronic cocaine exposure are well described (12). So far, the effects of cocaine on DNAm were most extensively studied in the rodent brain, which has led to the identification of cocaineassociated differentially methylated genes. Interestingly, several of these genes were previously related to the neurobiology of addiction, such as Bdnf (13) and Pp1c (14) as well as the immediate early genes c-Fos (15) and fosB (14). Cocaine was also shown to affect methylation and expression levels of DNA methyltransferases (Dnmt1, Dnmt3a) (15)(16)(17)(18) and Tet demethylases (19) thereby directly altering DNAm homeostasis.
In contrast, much less is known about cocaine-associated DNAm changes in humans. Due to an easy and minimally invasive sampling procedure, epigenome-wide association studies (EWASs) of different cocaine use phenotypes have been performed in human blood samples (20, 21). These studies revealed an enrichment of cocaine-associated differential methylation within biological pathways involved in synaptic transmission, neurogenesis, cell division, and chromatin organization. However, as DNAm levels are known to be tissue-specific, the predictive value of differential methylation in blood for methylation signatures of CUD in the brain might be limited. Thus, profiling DNAm levels in brain tissue is of critical importance. So far, two epigenome-wide characterizations in CUD have been performed in human post-mortem brain tissue from the caudate nucleus (CN) (22) and the nucleus accumbens (NAc) (23). By analyzing post-mortem tissue of the CN from N = 25 CUD cases and N = 25 controls using a reduced representation bisulfite sequencing approach, a total of 173 differentially methylated regions (DMRs) were found that were significantly associated with CUD status. The IRX2 gene was among the top findings of CUD-associated differential methylation in the CN. Using a similar approach in the NAc, the key finding was a hypermethylated region within the tyrosine hydroxylase (TH) gene encoding an enzyme involved in dopamine metabolism. Interestingly, the identified DMR in TH was shared between the CN and the NAc. Both studies provide some first mechanistic insights into the biological relevance of CUD-associated DNAm changes in the human brain.
While epigenome-wide alterations associated with CUD have been studied in the human striatum, they have not been reported for cortical brain regions so far. The prefrontal cortex (PFC) has a key role in the addiction cycle due to its relation to dysregulated cognitive processes in SUDs such as inhibitory control and response to drugrelated cues (24). In CUD, resting-state connectivity alterations have been found in frontostriatal circuits representing functional connections between the PFC and the striatum (25). As DNAm alterations in the PFC have been described for the Brodmann area 9 (BA9) subregion in other SUDs such as alcohol use disorder (AUD) (26,27) and opioid use disorder (OUD) (28, 29), BA9 depicts a Frontiers in Psychiatry 02 frontiersin.org promising brain region for the epigenome-wide profiling of DNAm in CUD. Based on epigenome-wide DNAm data, epigenetic clocks provide the opportunity to estimate biological age, a measure that increases naturally with chronological age, but is strongly influenced by disease processes such as chronic inflammation (30). Two commonly used methodologies are Horvath's (31) and Levine's et al. (30) (PhenoAge) epigenetic clocks. Epigenetic age acceleration (EAA), defined as the difference between the estimated epigenetic and chronological age, is of particular interest in the context of diseases, with its interpretation as faster or slower aging considering chronological age as the reference. Existing literature reports on a trend toward positive EAA in the PFC of individuals with SUDs such as AUD and OUD compared to a control population, especially for Levine's clock (28,32). This implies that the SUD brain appears to be epigenetically older compared to a brain from an individual without SUD at the same chronological age.
In the present study of CUD in human post-mortem brain tissue of BA9, we aimed to expand the epigenome-wide characterization of CUD to the PFC as another key brain region in addiction. Next to performing an EWAS of CUD in BA9, we aimed to identify functional enrichment of CUD-associated differential methylation within biological pathways and utilized epigenetic clocks to investigate the relationship between chronological and biological age in CUD.

Post-mortem human brain samples
Post-mortem human brain tissue of BA9 from N = 42 male individuals was received from the Douglas Bell Canada Brain Bank (DBCBB, Montréal, QC, Canada). Brain samples were obtained as fresh-frozen tissue and were stored at −80 • C until processing. All sample donors were of European ancestry. At the brain bank, phenotypic information was routinely obtained using a psychological autopsy including next-of-kin interviews complemented by the coroner's notes and the medical record of the deceased donor. Psychiatric diagnoses related to substance use and comorbidities were based on the Diagnostic and Statistical Manual of Mental Disorders (DSM)-IV classification system, however, we used the more recent nomenclature of DSM-5, i.e., SUD instead of substance dependence throughout this study. We included subjects with donor age > 18 and information on CUD status (yes = CUD, no = no CUD). Exclusion criteria were severe neurodevelopmental and neuropsychiatric disorders, SUDs other than CUD or AUD, and psychiatric diagnoses except for major depressive disorder diagnoses (MDD, depressive disorder NOS). Matching was performed by age and comorbidities resulting in N = 21 CUD case/control pairs.

DNA extraction and methylation data analysis
Extraction of DNA from N = 42 BA9 post-mortem brain samples was performed using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). DNA was stored at −20 • C and epigenome-wide DNAm was determined using the Illumina Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA, USA). Samples were randomized based on CUD case/control status and comorbidities such as depressive disorders and AUD prior to the array processing steps.
Quality control (QC) was performed on raw data using a modified version of the CPACOR-pipeline (33) as described in Zillich et al. (27). All statistical analyses were performed using the statistical computing environment R v4.2.1 (34). In brief, sample exclusion criteria were missing rate > 0.10 and discordance between phenotypic and predicted sex based on DNAm. cytosine-phosphateguanine (CpG) sites were excluded if (1) the call rate was below 0.95, (2) they are located on sex chromosomes, and (3) their genomic position is related to a genetic variant with MAF > 0. 10.
Raw intensities were quantile-normalized followed by a conversion to beta-values of methylation. Beta-values were logittransformed to M-values as recommended by Du et al. (35). To assess DNAm signatures of CUD in post-mortem brain tissue, we applied a linear regression model using M-values as the dependent variable and CUD status as the predictor. Donor age, post-mortem interval (PMI), brain pH, diagnosis of a depressive disorder, and AUD status were included as covariates in the regression model. We performed a variance partition analysis using the R package variancePartition v1. 26.0 (36) to confirm the covariates (see also Supplementary Figure 1). As variance explanation by cause of death was only minimal and would have made the statistical model significantly more complex, we did not include it as a covariate. Based on the methylation data, we estimated the fraction of neuronal and non-neuronal cells with the Houseman (37) approach using the dorsolateral PFC (dlPFC) dataset (38) as reference. We included the fraction of neuronal cells as a covariate to correct for cell-type heterogeneity in the samples. To adjust for technical effects from microarray processing, we derived principal components from the internal control probes of the EPIC array and included the first 10 principal components (cpPC1-10) into the statistical model. This resulted in the following linear model: DNAm ∼ CUD status + age + PMI + pH +depressive disorder status + AUD status We performed Benjamini-Hochberg false discovery rate (FDR) correction (39) to adjust for multiple testing. Annotation of CpG sites to genes and genomic background was based on the Illumina manifest as downloaded on 10th August 2018: http://webdata.illumina.com.s3-website-us-east-1.amazonaws. com/downloads/productfiles/methylationEPIC/infiniummethylationepic-v-1-0-b4-manifest-file-csv.zip.
To test for enrichment of CUD-associated CpG sites within categories of gene regulatory features (TSS1500, TSS200, 5 -UTR, 1st Exon, ExonBnd, Body, 3 -UTR, and intergenic regions) and the CpG island (CGI) background (N-Shelf, N-Shore, Island, S-Shore, S-Shelf, and open sea regions), we performed chi-squared tests against the distribution of features among all probes on the EPIC array. The results of the enrichment analysis were adjusted for multiple testing using Bonferroni correction.

Differentially methylated regions
Following the assessment of CUD-associated methylation differences at single CpG sites, we investigated DMRs associated with CUD status. We used the dmrff algorithm v1.0.0 (40) for the identification of DMRs which was recently shown to outperform several other DMR identification tools when applied to methylation array data (41). We applied dmrff on the EWAS summary statistics (N = 654,448 CpG sites). The maximum distance between individual CpG sites (maxgap) was set to 500 bp and an association p-value threshold (p. cutoff) of 0.01 was used. The implemented Bonferroni correction was used to adjust for multiple testing during the identification of DMRs.

Replication of CpG-sites identified in the striatum and in whole blood
We further tested whether DNAm signatures, which had previously been associated with CUD in whole blood (20), the CN (22), and the NAc (23) were also associated with CUD in the present sample of BA9. For DMRs overlapping between the two studies in the striatum, we tested CpG sites in close proximity (500 bp upand downstream) for their association with CUD in our sample. The comparison with CUD-associated CpG sites derived from whole blood was restricted to the subset of CpG sites that was available in our analysis in BA9.

Gene-Ontology enrichment analysis
To investigate the potential biological relevance of CUD-associated differential methylation, we analyzed the overrepresentation of differentially methylated CpG sites within biological pathways. Under the assumption of hypomethylation leading to increased gene expression and vice versa, we applied missMethyl v1.30.0 (42) separately to hypermethylated and hypomethylated CpG sites after selection for a CUD association p-value < 0.001. The set of N = 654,448 CpG sites remaining after QC was used as the reference set during the enrichment analysis. Further, we performed an enrichment analysis on DMRs resulting from dmrff using the "goregion" function in missMethyl. Correction for multiple testing in results from missMethyl was performed using FDR.

WGCNA
To further characterize the biological and functional relevance of CUD-associated methylation signatures in the brain, we performed a weighted correlation network analysis using WGCNA v1.71 (43). WGCNA identifies co-methylation modules and provides a correlation estimate between modules and traits such as CUD status. We prioritized CpG sites based on their annotation to promoter/transcription start sites (TSS 200 and TSS1500) under the assumption of a cis-regulatory role of their methylation levels on gene expression. Quantile normalized beta values of the resulting 117,876 CpG sites were used as input data for WGCNA. Soft power threshold picking was performed according to the R signed 2 > 0.90 criterion resulting in a power threshold of five for the methylation data in BA9. Parameters for constructing the signed network in a block-wise approach were minModuleSize = 30, maxBlockSize = 36,000, and mergeCutHeight = 0.25. For the CUD-associated co-methylation modules, enrichment analysis was performed in missMethyl using the N = 117,876 input CpG sites as background. Module hub genes were derived by annotation of CpG sites to genes and were defined as the top 0.5% of annotated genes in a co-methylation module ranked by the product of module membership (correlation of module eigengene with methylation profile) and gene significance (correlation between CpG site and trait). Gene significance was calculated for each CpG site. Protein-protein interaction (PPI) networks of hub genes from the CUD-associated co-methylation modules were generated in Cytoscape v3.9.1 (44) using the stringApp v1.7.0 (45). Networks of the "full STRING network" type were generated by applying the STRING protein query option on the hub gene list. An interaction score threshold of 0.7 was used to include only high-confidence interactions. PPI networks were visualized using the "Prefuse Force Directed Layout" option while unconnected nodes ("Singletons") and hidden edges were removed from the visualization. Information on stringdb identifiers corresponding to the nodes of the final PPI network for each module is provided in Supplementary Tables 7A-D.

Epigenetic age estimates
Using quantile-normalized beta-values of the BA9 methylation data, we estimated EAA using epigenetic clocks as implemented in the R package methylclock v1.2.1 (46). We applied Horvath's (31) (Horvath) and Levines's (30) (PhenoAge) epigenetic clocks to assess epigenetic age in the context of CUD directly in the brain. We first investigated the relationship between epigenetic age and chronological age using correlation analysis. Next, we estimated EAA (ageAccel), defined as the residuals of regressing chronological age on epigenetic age. We performed t-tests to assess the group difference between individuals with and without CUD and linear regression models to evaluate the association with CUD status while adjusting for the covariates neuronal cell fraction, pH, PMI, depressive disorder, and AUD status.

Ethics approval statement
Post-mortem brain tissue sampling at the DBCBB was performed according to their established ethical standards including written informed consent from the next-of-kin for each subject. Our study design was approved by the Ethics Committee II of the University of Heidelberg, Medical Faculty Mannheim, Germany, under the register number 2021-681.

Results
Demographic data of sample donors are summarized in Table 1, while a detailed overview of toxicology and cause of death is provided in Supplementary Table 1. After QC of methylation data, all N = 42 subjects, and a total of 654,448 CpG sites remained for further analysis.

Epigenome-wide association study of CUD
None of the CpG sites passed epigenome-wide significance (p < 7.64e-08). The strongest association with CUD in BA9 was detected for the CpG site cg04659218 (p = 7.98e-07, FDR = 0.52) on chromosome 11. A Manhattan plot summarizing the results is displayed in Figure 1A.

Differentially methylated regions associated with CUD status
A total of 20 DMRs were significantly associated with CUD status after multiple testing correction. We detected more hypermethylated (N = 17, 85%) than hypomethylated (N = 3, 15%) DMRs. The DMR consisting of six hypermethylated CpG sites annotated to the genes LOC101927196 and Fibrous Sheath Interacting Protein 2 (FSIP2) displayed the strongest association with CUD status in BA9 (p = 1.74e-16). Results of the DMR analysis including summary statistics for all 20 identified regions are visualized in Table 2. 3.3. Replication of CpG sites identified in whole blood and striatal post-mortem brain tissue Of the 186 CpG sites associated with cocaine and crack cocaine use in whole blood, 138 were available in our dataset after QC. Six CpG sites annotated to the genes PDSS2, EHMT2, C2, PDE12, ARF4, ARF4-AS1, RMDN2, RMDN2-AS1, and CDKN1C were nominally associated in the present analysis (see also Supplementary Table 3A) but did not remain significant after multiple testing correction. The DMRs previously identified as being conserved between the CN and NAc annotated to the genes SPEG, TH, PDXDC1, LDHD, and SYT5 were not replicated in our sample of BA9. Similarly, CpG sites located in a genomic window 500 bp up-and downstream of these overlapping striatal DMRs were not significantly associated with CUD in BA9 (Supplementary Table 3B). However, when we compared the DMRs in BA9 with DMRs detected in the CN and the NAc individually, we found one hypermethylated CUDassociated DMR annotated to the Aldehyde Dehydrogenase 3 Family Member B1 (ALDH3B1) gene on chromosome 11 that was conserved between BA9 and the CN.

Functional annotation and enrichment analysis of CUD-associated CpG sites
To analyze enrichment within genomic features (Figure 2A) and the CGI background (Figure 2B), we prioritized our findings based on CpG sites associated with CUD at nominal significance (p < 0.05). This resulted in a subset of N = 29,176 CpG sites of which N = 14,737 were hypermethylated and N = 14,439 were hypomethylated. We detected a balanced ratio between hyper-and hypomethylation in this subset ( Figure 2C). For the CUD-associated CpG sites, a statistically significant enrichment within 1st Exon, N-Shore, and CGI regions was detected, whereas a significant depletion from intergenic and open sea regions was found.
The most significantly enriched Gene Ontology (GO) term was response to exogenous dsRNA (p = 3.71e-04) for hypermethylated CpG sites and regulation of calcium ion transmembrane transporter activity (p = 2.39e-04) for the subset of hypomethylated CpG sites. No enrichment remained statistically significant after multiple testing correction (Supplementary Table 4A). Descriptively, genes harboring hypermethylated CpG sites were overrepresented within biological pathways involved in cation and especially mechanosensitive calcium ion channel activity, cyclic adenosine monophosphate (cAMP)-dependent protein kinase activity (PKA), and immunological signaling mainly related to the cyclic GMP-AMP synthase (CGAS)-mediated response to nucleic acids ( Figure 2D). In contrast, the gene set derived from hypomethylated CpG sites was enriched within pathways related to calcium transmembrane transport, negative regulation of muscle relaxation, and cAMP binding ( Figure 2E). Also, positive regulation of inositol-1,4,5trisphosphate-sensitive calcium release channel activity was among the enriched GO terms detected for hypomethylated CpG sites. The strongest overrepresentation resulting from the enrichment analysis using DMRs was observed for "thiosulfate-thiol sulfurtransferase activity" (Supplementary Table 4B).

Weighted correlation network analysis
Due to the known cis-regulatory role of DNAm on gene expression, we extracted TSS1500-and TSS200-associated CpG sites from the methylation matrix to identify transcriptionally relevant comethylation networks. After hierarchical clustering in WGCNA, the sample tree was pruned at the height of 20 resulting in the removal of one sample from the analysis that was considered as an outlier (N = 41 remaining). The sample dendrogram was related to traits resulting in the clustering of CUD cases and non-CUD individuals within subgroups ( Figure 3A). Network construction led to the identification of a total of 69 modules ranging in size between 38 and 35,758 CpG sites. A module-trait correlation heatmap for all resulting modules is shown in Supplementary Figure 3.
Of the 69 modules, four were significantly associated with CUD status: brown (N = 6,309 CpG sites, p = 9.15e-03), brown4 (N = 187 CpG sites, p = 0.04), blue (N = 14,050 CpG sites, p = 0.02), and steelblue (N = 340 CpG sites, p = 0.02). A positive, highly significant correlation between module membership and gene significance for CUD status was detected within these modules ( Figure 3B). Whereas module brown was positively associated with CUD status (r = 0.40), a negative correlation was identified for modules brown4, blue, and steelblue (−0.37 ≤ r ≤ −0.32, see Figure 3C). In addition, for the brown module, significant negative correlations were detected with age (r = −0.5, p = 8.11e-04) and the fraction of neuronal cells (r = −0.91, p = 7.69e-17). The blue module also depicts a cell type-specific network displaying a highly significant positive correlation with neuronal fraction (r = 0.96, p = 1.78e-22). To further assess the biological processes related to CUD-associated modules, we performed GO enrichment analyses (Figure 3D). After multiple testing correction, GO enrichment in modules brown and blue remained statistically significant. Descriptively, module brown was enriched for synaptic signaling and neurotransmitter transport, whereas module brown4 revealed GO terms related to transcription factor binding and neuroblast as well as synapse maturation. The blue module displayed enrichment of large GO terms involved in basic cellular processes and compartments, such as extracellular region, plasma membrane, and cell adhesion. Within module steelblue, nominally significant pathway enrichment was related to ion channel activity, membrane potential homeostasis, and vesicle transport. Full enrichment statistics for the top enriched GO terms in each module are provided in Supplementary Table 5.
To further characterize the molecular processes enriched in the CUD-associated modules, we constructed PPI networks based on module hub genes (listed in Supplementary Table 6). The resulting PPI network plots are displayed in Supplementary Figures 4A-D. Among PPI network hub nodes of modules brown and blue (Supplementary Table 8), we found members of the Transforming Growth Factor Beta (TGF-β) signaling cascade such as TGF Alpha (TGFA), TGFB1, and Endoglin (ENG). Reflecting the pathway enrichment for neuron-related biological processes in modules brown, brown4, and steelblue, we detected multiple hub nodes with an important role in neuronal activity and function within these

Epigenetic age in CUD
A moderate (R Levine 2 = 0.44) to strong (R Horvath 2 = 0.80) correlation was found between estimated epigenetic age and chronological age (Figures 4A, B). Levine's epigenetic clock produced strongly negative PhenoAge estimates ( Figure 4B). When assessing ageAccel we observed a trend toward a more positive ageAccel in CUD cases compared to individuals without CUD, however, not statistically significant (p Horvath = 0.39, p Levine = 0.15, Figures 4C,  D). For both epigenetic clocks, we consistently found positive associations for CUD status while adjusting for covariates in the regression models for ageAccel, although not statistically significant (Supplementary Table 9).

Discussion
We performed an epigenome-wide analysis of DNAm differences in BA9 post-mortem human brain tissue from N = 21 CUD cases and N = 21 individuals without CUD. In the BA9 subregion of the human PFC, none of the associations of CpG site methylation levels with CUD status passed the level of epigenome-wide significance. However, we identified 20 DMRs significantly associated with CUD.
Among the strongest single-site associations with CUD status, there were multiple CpG sites annotated to genes involved in transcriptional regulation either via intrinsic transcription factor activity or due to the interaction with RNA POL-II. The transcription factor ERG is involved in the regulation of vasculature integrity (47). For ZBTB4, binding to methylated CpG sites was reported (48) and ZBTB22 is assumed to be involved in transcription regulation via interaction with POL-II (49). In our EWAS, we also detected CUD-associated differential methylation within the largest subunit of POL-II (POL2RA), suggesting a role of CUD-associated methylation signatures in gene expression regulation. Results of the enrichment analysis for CUD-associated CpG sites within genomic features are in line with this: the significant enrichment of CUD-associated differential methylation within features containing transcription factor binding sites and the depletion from intergenic and non-CGIassociated (open sea) regions points toward a regulatory potential of CUD-associated DNAm patterns on gene expression. Results from the epigenome-wide association study (EWAS) of CUD were prioritized based on association p-values resulting in a subset of N = 29,176 CpG sites with p assoc < 0.05 that was further analyzed by annotation to panel (A) gene regulatory features and (B) the CpG island (CGI) background. The proportion of CUD-associated CpG sites (blue) within different groups was compared to the EPIC array background (gray). Differences remaining statistically significant after Bonferroni multiple testing correction are highlighted using asterisks (* = p adj < 0.05, ** = p adj < 0.01, *** = p adj < 0.001). For several of the genes harboring DMRs significantly associated with CUD in BA9, there is evidence for their relevance for neuronal function. The DMR displaying the strongest association with CUD status was related to the LOC101927196 gene. Overexpression of the long non-coding RNA LOC101927196 was shown to reduce cell proliferation and support apoptotic processes in the brain in a rat model of autism (50). For NPFFR2 and KALRN, we found direct links to their importance for the biological response to cocaine. NPFFR2 interferes with the opioid system in the brain and its activation contributes to anxiety-like behavior in mice (51). Intraventricular application of NPFF reduced cocaine-induced conditioned place preference in rats pointing toward an important role of NPFF signaling in cocaine reward processing (52). The KALRN gene encodes several isoforms and Kalirin-7 was shown to be involved in the development and function of dendritic spines (53). After repeated cocaine exposure in a mouse model, expression levels of Kalirin-7 were elevated in the NAc contributing to the formation of new dendritic spines. Kalirin-7 knock-out mice displayed increased cocaine-induced locomotor sensitization and a decreased conditioned place preference (54). Looking at multiple brain regions allows a more complete understanding of the role of DNAm alterations in CUD. In line with results from previous studies on epigenome-wide alterations of CUD in the CN (22) and in the NAc (23), we observed more hypermethylation than hypomethylation at the DMR level in BA9. In our systematic comparison of the DMRs that were consistently associated with CUD in both the CN and NAc, none were detected in our analysis. We detected one DMR annotated to the ALDH3B1 that was conserved between BA9 and the CN. Hypermethylation of this DMR was also consistent between studies. It has to be noted that the EPIC array covers only a finite number of CpG sites which may limit the comparison with results from reduced representation bisulfite sequencing. In addition, we also compared our results with DNAm signatures of CUD in whole blood. While some nominally significant overlap was observed, there was little evidence for a systematic relationship between the CUD-associated CpG sites in BA9 and whole blood, as the overlap was not significant after multiple testing correction. DNAm is highly tissue-specific, which could partly explain the missing overlap. At the same time, the samples were phenotypically different in age and substance use.
In our analysis in BA9, no enrichment of hypo-and hypermethylated CUD-associated CpG sites among GO terms remained statistically significant after multiple testing correction. Further, pathway enrichment was driven by relatively few genes leading to the emergence of multiple biologically related pathways among the results. The functional relevance of CUD-associated differential methylation should thus be evaluated together with the enriched biological processes and hub genes derived from WGCNA co-methylation modules. Three of the four CUD-associated modules were enriched for neuron-specific processes such as synaptic signaling, ion channel activity, and neurotransmitter transport. Next to the association with CUD, modules brown and blue displayed a strong correlation with neuronal fraction, which might also contribute to the enrichment of neuronal processes within these comethylation modules. PPI networks based on module hub genes revealed CACNA1C and NR3C1 as highly connected genes for which a previous relation to cocaine use is known. CACNA1C is a well-described risk gene for psychiatric disorders such as bipolar disorder, MDD, and schizophrenia (55). By encoding a subunit of the voltage-gated Ca v 1.2 calcium channel, it depicts a key regulator of neuronal excitability. Further, the expression of cocaine-induced locomotor sensitization was shown to be dependent on Ca v 1.2 signaling in mice (56). The NR3C1 gene encodes a glucocorticoid receptor for which downregulation in blood was observed in chronic cocaine users (57) and its DNAm levels were shown to predict substance use phenotypes in adolescents (58). Further, brain-specific conditional knock-out of NR3C1 flattened the dose-response curve for cocaine self-administration and reduced behavioral sensitization in an intravenous cocaine self-administration model in mice (59).
Thus, CACNA1C and NR3C1 depict important regulators of the reinforcing effects of cocaine. We further identified JUN as the hub node with the strongest connectivity in module steelblue. By forming heterodimers with FOS family proteins, JUN is involved in the formation of the AP-1 transcription factor complex, a key regulator of neuroplasticity in addiction (60). In a mouse model, decreased cocaine-induced conditioned place preference was found during expression of a dominant-negative mutant of c-Jun in the brain (61). Epigenetic dysregulation of the JUN interaction network might therefore contribute to neurocircuit alterations in CUD.
In our analysis of epigenetic age in the brain, we observed a trend toward a more positive ageAccel in CUD cases compared to individuals without CUD that was more pronounced with Levine's clock compared to Horvath's clock. This is in line with results from previous studies that report on positive ageAccel, particularly for Levine's epigenetic clock, in the PFC of individuals with opioid intoxication, AUD, and OUD (28,32). As cocaine intake is associated with neuroinflammatory processes in the brain (62, 63), ageAccel in CUD could be a consequence of an inflammatory response that is particularly reflected by Levine's (30) epigenetic clock. However, as comorbid diseases are frequently observed in CUD, an increase in ageAccel could at least in part be driven by systemic diseases that are themselves associated with an accelerated epigenetic age. Due to the lack of information on somatic comorbidities, we were not able to address this point in the present study. Further, for the interpretation of epigenetic age in the brain, it must be noted that epigenetic clocks are at least to some extent tissue-specific and might not be applicable to every tissue. We only assessed epigenetic clocks that were either trained on non-blood tissues (Horvath) or at least showed a high correlation of findings across tissues (Levine). The negative values of Levine's PhenoAge we observed in some of the brain samples might represent a post-mortem brain tissue artifact, as Levine's clock was initially trained on blood methylation data. Similarly, in a recent study in post-mortem human brain tissue, negative PhenoAge estimates were observed for Levine's clock (32). Within-subject comparisons using blood and brain tissue from the same person should be performed to investigate such potential tissue artifacts. Further, analyses with larger sample sizes, careful adjustment for comorbidities, and appropriate estimation methods are required to evaluate if accelerated epigenetic aging in the brain is a true phenomenon in CUD.
Several limitations apply to our study on DNAm signatures of CUD in post-mortem human brain tissue. First, our sample size of N = 42 is comparatively low for an epigenome-wide analysis and is associated with limitations in statistical power. Another limitation is the frequent comorbidity of CUD with other SUDs and mood disorders, which is well reflected in our sample. Drug intoxication at death and medication prior to death could have affected methylation levels and therefore might confound the results for CUD. Next to the CUD phenotype, assessment of used substances and comorbidities including their treatment is crucial to adjust for confounding effects on DNAm signatures. Due to the small sample size, we were not able to adjust for drug intoxication at death and medication prior to death. Thus, we cannot rule out a potential confounding of the identified CUD-associated methylation signatures by epigenetic effects of other substances. Third, the case-control design of the EWAS does not allow the assessment of temporal changes in DNAm levels. Thus, within the present study, we are not able to separate differential methylation acquired during the disease course of CUD from an epigenetic predisposition to CUD.
In the present study, we identified CUD-associated DNAm signatures in the BA9 subregion of the human PFC. We detected CUD-associated differential methylation within several transcription factor genes, an enrichment within genomic features of transcription regulation, and an association with gene regulatory networks involved in neurotransmission and neuroplasticity. Specifically, epigenetic alterations related to calcium and cAMP signaling as well as the CACNA1C and JUN PPI networks might contribute to the development of the altered neurocircuits observed in the CUD brain.
To follow up on these hypotheses, further studies are required. These should be based on larger sample sizes with detailed phenotyping to minimize the confounding of results in association studies. Further, integration of epigenome-wide methylation data with other omics data such as transcriptomics or proteomics is necessary to prioritize and substantiate findings from the EWAS. Such multi-omics studies depict promising approaches for deeper profiling of altered biological processes in SUDs as recently shown for AUD (64) and OUD (65). To address the temporal dynamics of molecular changes in the disease course of CUD, a longitudinal study assessing different stages of CUD within the same individuals should be performed. While a longitudinal design is nearly impossible using post-mortem brain tissue, it is a feasible approach for blood samples. The application of such innovative methods could pave the way for the development of precision medicine approaches thereby addressing the urgent need for novel therapeutic options in CUD.

Data availability statement
Methylation data presented in this study have been deposited at the European Genome-phenome Archive (EGA), which is hosted by the EBI and the CRG, under accession number: EGAS00001006826 (https://ega-archive.org/studies/EGAS00001006826).

Ethics statement
The studies involving human participants were reviewed and approved by the Ethics Committee II of the Medical Faculty Mannheim, Heidelberg University, Mannheim, Germany under the register number: 2021-681. Written informed consent was not provided because post-mortem brain tissue was obtained from individuals that died by suicide or sudden death. Following the established ethical standards of the DBCBB, written informed consent was obtained from the next-of-kin for each subject.

Funding
Funding was provided by the German Federal Ministry of Education and Research (BMBF) through the e:Med research program "A systems-medicine approach toward distinct and shared resilience and pathological mechanisms of substance use disorders" (01ZX01909 to RS, MR, SW, and AH). Further, by "Toward Targeted Oxytocin Treatment in Alcohol Addiction (Target-OXY)" (031L0190A to JCF), the ERA-NET program: Psi-Alc (FKZ: 01EW1908), the Hetzler Foundation for Addiction Research (to AH), and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID 402170461-TRR 265 (66) to RS, AH, and MR.