Neuronal-specific methylome and hydroxymethylome analysis reveal significant loci associated with alcohol use disorder

Background: Alcohol use disorder (AUD) is a complex condition associated with adverse health consequences that affect millions of individuals worldwide. Epigenetic modifications, including DNA methylation (5 mC), have been associated with AUD and other alcohol-related traits. Epigenome-wide association studies (EWAS) have identified differentially methylated genes associated with AUD in human peripheral and brain tissue. More recently, epigenetic studies of AUD have also evaluated DNA hydroxymethylation (5 hmC) in the human brain. However, most of the epigenetic work in postmortem brain tissue has examined bulk tissue. In this study, we investigated neuronal-specific 5 mC and 5 hmC alterations at CpG sites associated with AUD in the human orbitofrontal cortex (OFC). Methods: Neuronal nuclei from the OFC were evaluated in 34 human postmortem brain samples (10 AUD, 24 non-AUD). Reduced representation oxidative bisulfite sequencing was used to assess 5 mC and 5 hmC at the genome-wide level. Differential 5 mC and 5 hmC were evaluated using the methylKit R package and significance was set at false discovery rate < 0.05 and differential methylation > 2. Functional enrichment analyses were performed, and gene-level convergence was evaluated in an independent dataset that assessed 5 mC and 5 hmC of AUD in bulk cortical tissue. Results: We identified 417 5 mC and 363 5hmC significant differential CpG sites associated with AUD, with 59% in gene promoters. Some of the identified genes have been previously implicated in alcohol consumption, including SYK, DNMT3A for 5 mC, GAD1, DLX1, DLX2, for 5 hmC and GATA4 in both. Convergence with a previous AUD 5 mC and 5 hmC study was observed for 28 genes. We also identified 5 and 35 differential regions for 5 mC and 5 hmC, respectively. Lastly, GWAS enrichment analysis showed an association with AUD for differential 5 mC genes. Discussion: This study reveals neuronal-specific methylome and hydroxymethylome dysregulation associated with AUD, identifying both previously reported and potentially novel gene associations with AUD. Our findings provide new insights into the epigenomic dysregulation of AUD in the human brain.


Introduction
The detrimental effects of alcohol use disorder (AUD) are substantial, resulting in more than 150,000 deaths globally (GBD, 2019Risk Factors Collaborators, 2020).AUD is characterized by persistent, uncontrollable, and excessive alcohol consumption despite its negative consequences.Although genome-wide association studies (GWAS) have identified genetic risk factors of AUD (Gelernter et al., 2019;Zhou et al., 2020;2021;2023), these only account for a portion of the variation observed.
Epigenetic mechanisms, such as DNA methylation (5mC), have been implicated in AUD and alcohol-related traits in human studies evaluating various tissues, including saliva, blood, and brain (Longley et al., 2021;Clark et al., 2022;Montalvo-Ortiz et al., 2022;Zillich et al., 2022).The 5mC mechanism involves the addition of a methyl group to the carbon 5 position of the nucleotide, which is catalyzed by DNA methyltransferases (DNMTs) (Gibney and Nolan, 2010).DNA hydroxymethylation (5hmC) occurs when this methyl group is removed through oxidation catalyzed by a family of ten-eleven translocase proteins (TET1, TET2, and TET3) during the DNA demethylation process.Recent work from our group and others has shown that 5hmC is functionally distinct from 5mC.This epigenetic mark is associated with transcriptional activation and highly prevalent in the brain (Rompala et al., 2022).Several studies have implicated 5hmC in anxiety-related behaviors (Papale et al., 2017), schizophrenia, bipolar disorder (Guidotti et al., 2014;Madrid et al., 2016), autism (Papale et al., 2015), and Alzheimer's disease (Kuehner et al., 2021).Interestingly, a recent study evaluating 5mC and 5hmC in bulk tissue from the human postmortem brain identified a role for 5hmC in AUD (Clark et al., 2022).
Epigenetic patterns, such as 5mC and 5hmC, are tissue-and celltype specific, and particularly 5hmC is highly enriched in the brain and abundant in neuronal cells, underscoring the need to investigate this epigenetic mark in brain tissue, particularly in neurons (Kriaucionis and Heintz, 2009;Szulwach et al., 2011;Mellén et al., 2012).However, most epigenetic studies have used bulk brain tissue, which can mask cell-type specific biological signals, highlighting the need for a cell-type-specific approach when evaluating the epigenetic landscape of AUD in the human brain.
The orbitofrontal cortex (OFC) has been implicated in decision making and motivated reward-related behavior (Morisot et al., 2019;Lohoff et al., 2021), and recent neuroimaging studies have associated alterations in this brain region with AUD (Shields and Gremel, 2020;Bracht et al., 2021;Atmaca et al., 2023).Individuals diagnosed with AUD exhibit a reduction in the OFC volume, accompanied by a decrease in gray matter, and an impact on dopaminergic pathways (Volkow et al., 2007;Coleman et al., 2011;Le Berre et al., 2014;Nimitvilai et al., 2017;Moorman, 2018;Morisot et al., 2019;Hernandez and Moorman, 2020).Recent 5mC studies from our group and others have revealed a role of epigenetic mechanisms in OFC in the context of substance use disorders (SUDs) (Kozlenkov et al., 2017;Rompala et al., 2022).
In this study, we examined neuronal-specific 5mC and 5hmC profiles in the OFC of AUD (n = 10) and non-AUD (n = 24) groups to identify epigenetically dysregulated genes and evaluate the differences between 5mC and 5hmC marks in the OFC.We also identify the functional pathways enriched by these epigenetically dysregulated genes, evaluate replication in an independent dataset, and assess its relationship with GWAS studies.

Study cohort
Our study cohort comprised 34 postmortem brain samples obtained from the National posttraumatic stress disorder (PTSD) Brain Bank51 (NPBB) (Friedman et al., 2017), a brain tissue repository at the U.S. Department of Veterans Affairs (VA).Consisting of European American and African American men with a mean age of 41 (s.d ± 12) (Friedman et al., 2017).The tissue samples were collected after obtaining informed consent from the next-of-kin and processed as described by (Friedman et al., 2017).The clinical diagnosis followed the antemortem assessment protocol (AAP) and postmortem diagnostic assessment protocol (PAP) based on the DSM-IV criteria (Friedman et al., 2017).The samples were categorized into AUD and non-AUD groups.The AUD group included 10 donors with alcohol use disorder (AUD) history, which refers to those diagnosed with alcohol dependence or Cases (10) Controls ( 24) alcohol abuse.The non-AUD group included 24 donors without an AUD diagnosis.AUD and non-AUD groups were matched by (PTSD), opioid use disorder (OUD), and current smoking.Table 1 presents the demographic and clinical characteristics of the study cohort.

High-throughput bisulfite sequencing and data processing
Sequencing data were obtained by reduced representation oxidative bisulfite-sequencing (RRoxBS), carried out at the Weill Cornell Epigenomics Core (New York, NY).The library preparations for 5mC and 5hmC were made using Mspl digestion for 400 ng of gDNA with the Ovation RRoxBS Methyl-Seq library preparation kit (TrueMethyl oxBS; Tecan, Switzerland).Bisulfite conversion was followed by a single-end 1 × 50 bp sequencing with the Illumina NovaSeq6000 system (mean depth of 42.7 ± 1.5 (µ ± SEM) million reads per library).The hg38 genome reference was used for adapter trimming, alignment information, and mapping efficiency of the sequencing data using an in-house BiSeq pipeline (Garrett-Bakelman et al., 2015).

Differential methylation analysis
The MethylKit R package (Akalin et al., 2012) was used to conduct differential methylation (5mC) and hydroxymethylation (5hmC) analyses at CpGs.The samples were filtered using a read coverage above ×10 and under the 99.9th percentile.Normalization was performed using the median method, where the mean was used to calculate the scaling factor to reduce coverage bias in the statistical analysis.Differential 5mC and differential 5hmC analyses were performed using logistic regression with correction for overdispersion and chi-squared significance testing (Akalin et al., 2012).The covariates included in the model were ancestry, age of death, PTSD, OUD, smoking, and postmortem interval (PMI).The sliding linear model (SLIM) was used to fit the p-values to q-values (Wang et al., 2011).The significance level for differential 5mC and 5hmC was defined as q-value < 0.05 and a greater than 2% difference of 5mC and 5hmC between AUD and non-AUD groups.Genome-wide significant (GWS) sites were used for annotation.We used the Genomation R package (Akalin et al., 2015) to annotate CpGs in promoters, introns, exons, and intragenic regions as well as CpG island, CpG islands shores, CpG shelves, and open sea (Akalin et al., 2015).Gene name annotation was performed using Ensembldb (Rainer et al., 2019).
Differentially methylated regions (DMR) and differentially hydroxymethylated regions (DhMR) with a window of 1,000 bases and at least 10 CpG sites covering the region were also examined for both 5mC and 5hmC epigenetic marks, respectively.Significant DMRs and DhMR were defined based on a q value < 0.05.

Functional enrichment analysis
Genomation R package was used to perform genomic feature annotation (Akalin et al., 2015).Functional enrichment analysis was made using the gene annotation obtained by scan_region.plperl tool and the UCSC genome browser annotation databases.Genes with a distance greater than 1,500 bp were not considered for the enrichment analysis (Kuhn, Haussler and James Kent, 2013).To conduct functional enrichment analysis we used Metascape (Zhou et al., 2019), Enrichr (Kuleshov et al., 2016), and WebGestalt (Liao et al., 2019), which integrates databases such as NCA TS BioPlanet (Huang et al., 2019), Panther (Mi et al., 2013), Gene Ontology Consortium (Gene Ontology Consortium, 2015), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000).As the enriched terms across all databases may share relationships or exhibit redundancy, we also performed enrichment analysis networks using Metascape to visualize the interaction among the enriched terms (biological pathways).For this, we calculate the Kappa-test score of the enriched terms and hierarchically cluster them into non-redundant groups.This process involves selecting the most significant terms of each cluster and determining connections among them, guided by a Kappa similarity threshold above 0.3.Each node is color-coded for distinction based on the biological pathway they belong to.Protein-protein interaction (PPI) enrichment analysis was also conducted to examine interaction across genes mapping significant 5mC/5hmC CpG sites.The molecular code detection (MCODE) algorithm was used to cluster enrichment ontology terms to identify neighborhoods where proteins are densely connected in the following databases: STRING, BioGrid, and OmniPath (Türei et al., 2016;Oughtred et al., 2019;Szklarczyk et al., 2023).
A gene-level overlap between Clark et al. ( 2022) (Clark et al., 2022) reported genes (p < 1 × 10 −5 ) and our significant findings was assessed.For this, we employed the liftOver tool to convert genome coordinates from hg37 to hg38.Fisher's exact test was applied using the GeneOverlap 1.38.0R package (Li Shen and Shenli, 2017) to evaluate whether the gene-level overlap is statistically significant.

AUD-associated 5mC and 5hmC differential CpG sites
After QC, 1908163 and 1707251 CpG sites were included for 5mC and 5hmC differential analyses, respectively.For 5mC, we identified 417 CpG sites after multiple testing corrections with a difference in the percentage of methylation between AUD and non-AUD groups higher than 2. Of these, 137 were hypomethylated and 280 were hypermethylated (Figure 1).For 5hmC, we identified 363 CpG sites after multiple testing corrections, with 213 hypo-and 150 hypermethylated CpG sites (Figure 2).Tables 2, 3 list the top GWS differential 5mC and 5hmC CpG sites.The 5mC GWS annotation showed that 59% of them were located at the promoter region, 21% at the intragenic region, 14% at introns, and 6% at exons.For the 5hmC GWS, 66% were located at the promoter region, 17% at the intragenic region, 11% at introns, and 6% at exons.

Protein-protein interaction analysis
The PPI network analysis of the differential 5mC genes showed as significant pathways cell-cell adhesion (GO:0098609), homophilic cell adhesion via plasma membrane adhesion molecules (GO: 0007156), cell-cell adhesion via plasma-membrane adhesion molecules (GO:0098742).The MCODE cluster algorithm identified GO pathways related to neurogenesis, chromatin organization, and cell adhesion (Supplementary Figures S6, S7; Supplementary Tables S10, S11).
For the 5hmC marks, the significant pathways identified in the PPI network analysis were cell-cell adhesion homophilic cell adhesion via plasma membrane adhesion molecules (GO: 0007156), embryonic organ development (GO:0048568), and cell-cell adhesion via plasma-membrane adhesion molecules (GO:0098742).The MCODE cluster algorithm identified GO pathways implicated in neurogenesis, cell adhesion, calcium ion transport, and Wnt signaling (Supplementary Figures S8, S9; Supplementary Tables S12, S13).

Gene overlap analysis
We compared our findings with those reported by (Clark et al., 2022), which evaluated AUD-associated 5mC and 5hmC in bulk tissue from the human postmortem PFC.For 5mC, we found an overlap of 14 genes with their reported 576 (Figure 4A; Supplementary Table S8).Fisher's exact test showed that the overlap is significant (p = 8.5e-25, odds ratio = 124.00).Similarly, for 5hmC, we found a significant overlap of 14 genes with their reported 1,023 genes (Figure 4B; p = 3.9e-22, odds ratio = 79.80;Supplementary Table S9).

Differentially methylated and hydroxymethylated regions (DMR/ DhMR) analysis
After QC, a total of 38757 and 37316 regions were examined in the 5mC and 5hmC analyses, respectively.With a q value < 0.05, a total of 5 DMRs (Supplementary Table S14) and 35 DhMRs (Supplementary Table S15) were significantly associated with AUD.One of the DhMR maps to the TDP1 gene (chr14: 89954001-89955000-hg38), a gene also observed with differential 5hmC and nominally associated with AUD in the (Clark et al., 2022).

Discussion
This study presents a neuronal-specific epigenomic investigation of AUD in the human brain.We profiled 5mC and 5hmC at the genome-wide scale and revealed differential CpG sites associated with AUD that map genes previously reported with this disorder.Our analysis identified 417 and 363 GWS AUD-associated CpG sites for 5mC (Supplementary Table S1) and 5hmC (Supplementary Table S2), respectively.Our 5mC findings revealed genes related to AUD or other alcohol traits reported in previous studies.For example, SYK, this gene has been previously associated with alcohol metabolism in the liver and its inhibition is linked to reduced liver inflammation (Qu et al., 2018;Kurniawan et al., 2020).Furthermore, previous studies have found that binge drinking can induce SYK activation, and pharmacological inhibition of SYK significantly decreases alcoholic liver disease in a mouse model of binge drinking (Bukong et al., 2016).
Another differential 5mC gene with prior alcohol-related evidence is the DNMT3A, which encodes for an enzyme involved in de novo DNA methylation (Pino et al., 2017;Fischer et al., 2021).Ethanol exposure can induce a prolonged upregulation of DNMT3A in neuronal precursor cell lines and primary mouse embryonic fibroblasts (Miozzo et al., 2018).Similarly, Dnmt3a has been found upregulated in the nucleus accumbens of alcohol-preferring rats exposed to intermittent ethanol exposure (Niinep et al., 2021).
For 5hmC, one of the identified differential CpG sites previously linked to AUD mapped to GAD1, a gene that encodes for one of the glutamate decarboxylases that catalyze the conversion of glutamate to Gamma-aminobutyric acid (GABA).Chronic alcohol exposure has been linked to a decrease in GABA levels and an increase in GABA receptors (Sytinsky et al., 1975;Tran et al., 1981;Dodd et al., FIGURE 3 Generalized Gene-Set Analysis of GWAS for 5mC and 5hmC.The bar plot shows the gene-set analysis of GWAS for 5mC and 5hmC, including alcohol use disorder (AUD), problematic alcohol use (PAU), cannabis use disorder (CUD), opioid use disorder (OUD), and post-traumatic stress disorder (PTSD).Overlapped differential 5mC and 5hmC genes in human brain tissue.(A) shows the overlap between our 5mC GWS hits and 5mC reported findings from Clark et al. (2022), where 14 genes overlap between both results.(B) shows the overlap between our 5hmC findings from OFC and 5hmC reported findings from Clark et al. (2022), where 14 genes overlap between both results.
1996; Behar et al., 1999).GAD1 has also been found to be upregulated in the dorsomedial thalamus of human subjects in individuals with AUD (Hade et al., 2021).However, a recent study did not find a difference in GAD1 mRNA levels in OFC between individuals with AUD and the non-alcoholic group AUD (Edenberg et al., 2010;Underwood et al., 2019).More research is needed to fully elucidate the role of these genes in AUD.
Additional AUD-associated CpG sites with differential 5hmC were those mapped to DLX1 and DLX2, which are also implicated GABA signaling, specifically in interneuron GABA synthesis (Le et al., 2017;Pla et al., 2018) and play a role in interneuron synaptogenesis and dendritogenesis (Pla et al., 2018).Moreover, these genes are suggested to directly promote the expression of Grin2b, a gene previously reported in human stem-cell-derived cortical neurons exposed to chronic alcohol consumption and also reported in the prefrontal cortex and hippocampus of mice treated with chronic alcohol consumption followed by withdrawal (Endele et al., 2010;Xiang et al., 2015;Pla et al., 2018;Myers et al., 2019).In the current study, both DLX1 and DLX2 showed hyper-5hmC and were located in the exon region.DLX2 CpG site was located in the promoter region, suggesting that it could be directly impacting gene expression regulation in individuals with AUD.The effect of 5hmC in the exon region is not well understood; however, in the case of 5mC, it is suggested that the density of 5mC in the exon can enhance gene expression (Li et al., 2018).
GATA4 was identified in both 5mC and 5hmC differential analyses and has been previously associated with alcohol dependence (Treutlein et al., 2009;Edenberg et al., 2010;Karpyak et al., 2014).GATA4 encodes the GATA-motif binding protein type 4, a transcription factor that controls the expression of proteins involved in drug metabolism (Karpyak et al., 2014).In addition, high doses of alcohol increase its expression (Zhong et al., 2010).In our study, this gene was hypomethylated and hypo hydroxymethylated (5mC location at 11697675bp, 5hmC location at 11703255).The differential 5mC CpG site was located in the intronic region, suggesting chromosomal instability.The 5hmC CpG sites were located in the promoter and intronic regions.The presence of 5hmC at the promoter region has been associated with protection of gene transcription in regions where 5mC is present (Ehrlich and Ehrlich, 2014).This suggests an interaction of these epigenetic mechanisms in the modulation of GATA4 expression in the OFC.
In addition to differential analysis at CpG sites, we also assessed differential 5mC and 5hmC regions.We found 5 GWS DMRs, but none of the genes mapping to these regions have been previously linked to AUD or related psychiatric traits.Among the GWS DhMR identified, chr1:65265001-65266000, mapped to intron 1 of DNAJC6, a gene highly expressed in the brain and involved in synaptic vesicle uncoating.When comparing the identified regions with the differential CpG analysis, we observed an overlap between DMR and differential 5mC CpGs with the OSBP2 gene.The EWAS catalog (www.ewascatalog.org)have reported CpG sites in this gene associated with chronic obstructive pulmonary disease and type 2 diabetes based on a blood-based EWAS in the Generation Scotland cohort (n = 18,413).For 5hmC, we found an overlap for RP11-44N22.3, but its functionality is unknown.Lastly, when comparing the identified regions with Clark et al., 2022, we observed an overlap with the LRRC1 and TDP1 gene.CpG sites at the TDP1 were also identified in the differential 5hmC analysis.This gene participates in DNA repair of neuronal cells and disruptions on this gene cause a neurodegenerative disorder known as spinocerebellar ataxia with axonal neuropathy.More research is needed to understand the role of these genes in AUD.
The findings of our enrichment analysis are consistent with previous AUD-related studies identifying development and neurogenesis.For instance, in a previous report examining fetal alcohol syndrome (FAS), these pathways were found to be highly significant (Fischer et al., 2021).Moreover, studies have also reported that AUD impacts cell adhesion and neurogenesis, which involves the development of new neurons and their integration into functional neural networks (Ramanathan, 1996;Arevalo et al., 2008;Pino et al., 2017;Poulose et al., 2017;Lees et al., 2020;Wooden et al., 2021).The effects on developmental processes and neurogenesis may contribute to the cognitive impairment reported in individuals with AUD and may also be linked to brain dysfunction in various regions, including the OFC (De Wilde et al., 2007;Arevalo et al., 2008).
GWAS enrichment analysis showed enrichment of differential 5mC genes with AUD and PAU GWASes (Zhou et al., 2020).No enrichment was identified with GWAS of other SUDs, suggesting a specificity of our differential 5mC marks with alcohol-related traits.When comparing our findings with those reported by Clark et al. (2022), we observed 14 overlapping genes for each epigenetic mark, 5mC and 5hmC.Several overlapping genes with differential 5hmC have been previously linked to SUDs.For instance, KCNQ1, a gene encoding a potassium ion channel, was identified in a GWAS of alcohol dependence (Edenberg et al., 2010;Feng et al., 2022).Similarly, APBB2 has been associated with opioid and amphetamine dependence (Gelernter et al., 2014;Liu et al., 2018).Among the 5mC overlapping genes, ASIC2 has been linked to addiction-related behavior in mice (Kreple et al., 2014).These findings underscore the importance of evaluating both 5mC and 5hmC to fully investigate the role of epigenetic mechanisms on AUD, at least in the brain where 5hmC is highly prevalent and enriched in neurons.Future studies should explore these genes in greater depth to better understand their involvement in AUD.
The limitations of this study are that the donors of the study cohort present heterogeneous comorbidities.All individuals in the OUD group were also diagnosed with PTSD.We controlled this by using PTSD as a covariate in the differential 5mC and 5hmC analyses.In addition, we conducted a GWAS enrichment analysis of our 5mC and 5hmC annotated genes to determine whether the reported genes were enriched for the comorbidity traits, including PTSD and other SUDs.In addition, the cohort size is limited; however, it is comparable to other recently published postmortem brain studies.Another limitation is that all samples are male, limiting to identify the effect of sex in our results.The analyses were only carried out on CpGs sites and it would be important to conduct a study on non-CpGs sites, because of the role of 5mC at non-CpGs on neuropsychiatric diseases from our group and others (Jang et al., 2017;Nagamatsu et al., 2022).Further, since the technologies to assess 5mC and 5hmC in our study and Clark et al., 2022 are different, we were only able to evaluate the convergence between the two studies at the gene-level, not at the CpG level, which limits the interpretation in gene regulatory processes of the differential epigenetic marks identified.To determine if the epigenetic marks in the reported genes of this study are a cause of AUD or an effect, research in model organisms is warranted.
Our study characterized the methylome and hydroxymethylome profiles of AUD in neurons from the human OFC.Our results replicate previous findings in certain genes and highlight new findings for both 5mC and 5hmC.This study reveals new insights into the epigenomic dysregulation of AUD in the human brain and pinpoints potential drug targets for the treatment of individuals suffering from AUD.

FIGURE 1
FIGURE 1 5mC differential CpG sites associated with AUD.(A) Volcano plot shows the 5mC differential CpG sites associated with AUD.(B) Pie chart depicts the gene location of the GWS CpG sites identified.
FIGURE 2 5hmC differential CpG sites associated with AUD.(A) Volcano plot shows the 5hmC differential CpG sites associated with AUD.(B) Pie chart depicts the gene location of the GWS CpG sites identified.

TABLE 1
Demographic and clinical information of the study cohort.