Abstract
Alcohol Use Disorder (AUD) is a major contributor to global mortality and morbidity. Postmortem human brain tissue enables the investigation of molecular mechanisms of AUD in the neurocircuitry of addiction. We aimed to identify differentially expressed (DE) genes in the ventral and dorsal striatum between individuals with AUD and controls, and to integrate the results with findings from genome- and epigenome-wide association studies (GWAS/EWAS) to identify functionally relevant molecular mechanisms of AUD. DNA-methylation and gene expression (RNA-seq) data was generated from postmortem brain samples of 48 individuals with AUD and 51 controls from the ventral striatum (VS) and the dorsal striatal regions caudate nucleus (CN) and putamen (PUT). We identified DE genes using DESeq2, performed gene-set enrichment analysis (GSEA), and tested enrichment of DE genes in results of GWASs using MAGMA. Weighted correlation network analysis (WGCNA) was performed for DNA-methylation and gene expression data and gene overlap was tested. In the dorsal striatum, we discovered differential expression (FDR<0.05) for a total of 50 genes. In the VS, DE genes at FDR<0.25 were overrepresented in a recent GWAS of problematic alcohol use. The ARHGEF15 gene was upregulated in all three brain regions. GSEA in CN and VS pointed towards cell-structure associated GO-terms and in PUT towards immune pathways. The WGCNA modules most strongly associated with AUD showed strong enrichment for immune response and inflammation pathways. Our integrated analysis of multi-omics data sets provides further evidence for the importance of immune- and inflammation-related processes in AUD.
Introduction
Alcohol Use Disorder (AUD) is a major contributor to the global disease burden, with a prevalence of ∼17% among 12-month alcohol users in the US1, 2 and an estimated heritability of 49%3. Knowledge about the molecular mechanisms can foster understanding of causes and promote prevention. Recent genome-wide association studies (GWASs) have identified 29 genetic loci associated with Problematic Alcohol Use (PAU), a proxy of AUD4. While GWASs identify increasing numbers of disease-associated loci, the functional interpretation of many of these findings remains inconclusive. Analyzing the transcriptome can extend the understanding of the molecular mechanisms underlying AUD, by identifying associated gene expression patterns. Findings can in turn be integrated with results from GWASs and epigenome-wide association studies (EWASs) to identify the pathomechanisms underlying disease.
Processes in the central nervous system are considered to play a major role in the etiology of addiction, and the transition from chronic alcohol consumption to AUD5. Therefore, it is of particular interest to examine molecular changes associated with addiction in brain tissue. So far, only few studies have been conducted in postmortem human brain tissue to identify transcriptional changes associated with AUD6-8. These studies mainly focused on the prefrontal cortex (PFC) one important part of the neurocircuitry of addiction9, 10. The first transcriptome-wide study in the PFC found DE genes implicated in neuronal processes, such as myelination, neurogenesis, and neural diseases, as well as cellular processes, such as cell adhesion and apoptosis11. In Brodmann Area 9 downregulation of calcium signaling pathways has been observed in individuals with AUD compared to controls12. In the same study, a weighted gene co-expression analysis (WGCNA) pointed towards modules associated with AUD case/control status, which were enriched for nicotine and opioid signaling, as well as immune processes. Another study in the PFC (Brodmann Area 8) showed that co-expression networks associated with lifetime alcohol consumption were enriched for GWAS signals of alcohol dependence6.
Despite the importance of striatal regions in addiction processes, genome-wide human omics studies of these brain regions are still missing. The striatum is divided into the ventral striatum (VS), consisting of the nucleus accumbens and olfactory tubercle; and the dorsal striatum, comprising the caudate nucleus (CN) and putamen (PUT)13. The nucleus accumbens is involved in mediating motivational processes such as aversion and reward, which play a significant role in the development and maintenance of substance use disorders (SUD)13. In addition to regulating motor function, the CN and PUT are involved in cognitive processes relevant for addiction, such as executive functioning and cognitive control, reinforcement learning and habit formation14. Analyses of omics data from striatal regions could complement the knowledge on global molecular changes in the neurocircuitry of addiction in AUD.
In a recent EWAS of AUD in postmortem brain tissue, we identified differentially methylated CpG-sites and regions in the ventral and dorsal striatum15. Previous studies have shown the utility of integrating epigenetic and transcriptomic data in postmortem brain tissue of SUDs using weighted correlation network analysis (WGCNA)16. WGCNA clusters genes or CpG-sites into co-expressed or co-methylated modules based on correlation matrices. By relating modules to each other, WGCNA can be used for data integration, providing more insights than descriptive overlap. For example, whereas a descriptive comparison of histone H3 lysine 4 trimethylation (H3K4me3) and mRNA expression in individuals with AUD and cocaine use disorder revealed no consistent overlap between H3K4me3 trimethylation and gene expression17, a network analysis identified overlapping modules pointing towards co-expressed genes associated with H3K4me3 trimethylation6. Modules associated with AUD were enriched for CNS functions, such as synaptic transmission and regulation of neurogenesis6. WGCNA has also been used for integrating epigenetic and transcriptomic data and investigating their association with opioid use disorder (OUD) in postmortem human brain, identifying immune-related transcriptional regulation to be enriched in co-expressed and co-methylated modules18.
The aim of the present study was to investigate differential gene expression associated with AUD status in the ventral and dorsal striatum, relate these to GWAS findings, and to integrate the findings with DNA-methylation data using a network approach (WGCNA) in order to identify functionally relevant molecular mechanisms in AUD.
Materials and Methods
Samples
Postmortem human brain tissue from CN, PUT and VS of a total of 48 individuals with AUD and 51 control individuals was obtained from the New South Wales Tissue Resource Centre at the University of Sydney. After quality control (QC), the sample sizes for each brain region ranged between 63 and 77. Phenotypic information was assessed by next-of-kin interviews. Inclusion criteria for this study were: age>18 years, Western European Ancestry, no history of severe psychiatric or neurodevelopmental disorders, or SUDs other than AUD and nicotine use disorder. AUD was defined as meeting DSM-IV criteria for alcohol dependence and consuming 80g of alcohol a day or more (control group: <20g/day). Descriptive information can be found in Table 1 and Supplementary Table S1.
DNA extraction and Methylation Profiling
DNA extraction, methylation profiling, and QC was performed as described in Zillich, Frank 15. In brief, DNA was extracted using the DNeasy extraction kit (Qiagen, Hilden, Germany); the Illumina HumanMethylation EPIC BeadChip and the Illumina HiScan array scanning system (Illumina, San Diego, CA) were used to determine DNA-methylation levels. We used an updated and customized version of the CPACOR pipeline to extract beta values from raw intensities19. Criteria for the removal of samples and probes can be found in Zillich, Frank 15.
RNA extraction and -sequencing
RNA was extracted from frozen tissue according to the manufacturer’s protocol using the Qiagen RNeasy microKit (Qiagen, Hilden, Germany). The RNA Integrity Number (RIN) of all samples was determined using a TapeStation 4200 (Agilent, Santa Clara, CA). RIN values of 273 samples were larger than 5.5, for which libraries were prepared using the TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA) and which were sequenced on the NovaSeq 6000 (Illumina) at the Life & Brain Center in Bonn, Germany with read lengths of 2×100bp and a sequencing depth of 62.5 M read pairs per sample on average. Technical replicates were sequenced for all but four samples.
Statistical Analyses
All analyses apart from QC and read mapping were performed using R version 3.6.120. An overview of the analysis workflow can be found in Figure 1.
Mapping and Quantification
Sequencing quality was determined using FastQC21 and 24 samples (11 cases and 13 controls) were excluded due to insufficient sequencing quality (e.g. strong overrepresentation of sequences, GC distribution). Raw reads were mapped to the human genome (hg38) using HISAT2 (v.2.1.0)22. Quantification was performed with the featureCounts function of the Rsubread package (v.2.0.1)23, with hg38 annotation.
Differential Gene Expression Analysis
Differential gene expression was determined using DESeq2 (v.1.26.0)24. Minimal pre-filtering was applied, removing genes with normalized counts <10 for more than two samples. Technical replicates were merged prior to differential expression analysis using the collapseReplicates function as implemented in DESeq2. For the differential gene expression analysis, we included age, sex, RIN, pH-value of the brain, and postmortem interval (PMI) as covariates, because of their known influence on gene expression25-27. Results were filtered for differentially expressed (DE) genes with an absolute log2 fold change larger than 0.02. P-values were adjusted for multiple testing using Benjamini-Hochberg correction (FDR<0.05) 28. All downstream analyses were performed for DE genes with FDR<0.25.
Pathway Analysis
Gene-set enrichment analysis was performed using the R package fgsea (v.1.12.0)29, for which DE genes were ranked according to p-value. Enrichment analysis was performed for Gene-Ontology (GO) terms30 and Hallmark genesets31 and the results were adjusted using FDR correction.
WGCNA
Weighted correlation network analyses (WGCNA, v.1.70-3)16 were performed to identify modules of co-expressed genes and co-methylated CpG-sites. We assessed the relationship of these modules with AUD case/control status and tested the overlap between associated modules. WGCNA clusters the input matrix according to a dynamic tree cutting algorithm, using a soft power threshold that approximates the criterion of scale-free topology (Rsigned 2>0.80). Resulting soft power thresholds for expression networks were 6 for CN, 5 for PUT, and 14 for VS; for methylation networks, all power thresholds were 2.
To identify methylation networks associated with gene expression, beta values from normalized intensities of all samples from which gene expression data were available were filtered for promoter-associated CpG-sites based on the manufacturer’s manifest (Illumina, San Diego, CA). The resulting 105 796 CpG-sites were used as input.
For the RNA-seq data, count matrices were normalized using the DESeq2 function normalizeCounts and variance stability transformation was applied.
Networks were constructed using following settings: minimum module size=30, mergeCutHeight=0.25, maxBlockSize=36 000. In WGCNA, modules are labeled using colors. In the results section modules are labeled according to type of data, brain region, and color assigned in the analysis, e.g. “ e-VS-pink” for module “ pink” from the WGCNA analysis of gene expression data in the ventral striatum. For each module, its eigengene was calculated and correlated with AUD status. For modules associated with AUD status, we performed enrichment analysis using the GOenrichmentAnalysis function implemented in the WGCNA package for expression data and the R package missMethyl (v.1.20.4)32 for methylation modules.
Expression and Methylation Data Integration
To identify genes both DE and differentially methylated, we analyzed the overlap of DE genes (FDR<0.25) with the results of an EWAS (p<0.001) in the same sample15. We prioritized CpG-sites based on their functional relevance in gene expression regulation. Thus, promoter-associated CpG-sites were used in the analysis.
At the module level, gene-set overlap tests were performed using the R package GeneOverlap (v.1.22.0)33. Here, Fisher’s exact test is used to identify significant overlap. For each brain region, the overlap of the AUD-associated co-expression and co-methylation modules was tested.
GWAS Enrichment Analysis
We analyzed enrichment of DE genes with an FDR<0.25, and genes in AUD-associated WGCNA modules in GWAS summary statistics using Multi-marker Analysis of GenoMic Annotation (MAGMA, v.1.08b)34. We performed GWAS enrichment analysis for several SUDs, such as alcohol use disorder and problematic alcohol use4, cannabis use disorder35, and a recent GWAS comparing individuals with opioid use disorder with unexposed controls36. Results with p<0.05 were considered statistically significant.
Results
Differential Gene Expression
Gene expression analysis of postmortem brain tissue from AUD cases and controls revealed DE genes in all three investigated brain regions. DE genes were assessed at two different FDR thresholds (0.05 and 0.25). In the caudate nucleus, 1 326 DE genes were identified at FDR<0.25 (49 at FDR<0.05). Tubulin Tyrosine Ligase Like 4 (TTLL4, log2FC=0.11, p=2.3*10−8) and GATA Binding Protein 2 (GATA2, log2FC=-0.27, p=8.6*10−7) were the most significantly upregulated and downregulated genes, respectively. In the putamen, 84 DE genes were detected at FDR<0.25 (1 at FDR<0.05). Top DE genes were found to be Transcription Elongation Factor A Like 2 (TCEAL2, log2FC=0.09, p=5.8*10−5) and Desmin (DES, log2FC=-0.86, p=2.6*10−6). One gene, cardiomyopathy associated 5 (CMYA5) showed an upregulation in both caudate nucleus and putamen. Nine genes were downregulated in both dorsal striatal regions, with HLA-DOB having the highest log2FC in both regions. In the ventral striatum, 5 genes were DE at FDR<0.25, but none at FDR<0.05. Strongest differential gene expression in the ventral striatum was observed for Ankyrin Repeat And Ubiquitin Domain Containing 1 (ANKUB1) which was upregulated in AUD cases (log2FC=1.35, p=5.8*10−5). The most significantly downregulated gene in the VS of AUD cases was Caseinolytic Mitochondrial Matrix Peptidase Chaperone Subunit B (CLPB, log2FC=-0.11, p=5.2*10−6). ARHGEF15 (Rho Guanine Nucleotide Exchange Factor 15) was upregulated in all three brain regions at FDR<0.10. The Top 5 DE genes are listed in Table 2; DE genes significant at FDR<0.10 are listed in Supplementary Table S2 (CN), S3 (PUT), and S4 (VS). Overlap between DE genes in the different brain regions is shown in Figure 2A.
Pathway Analysis
Pathway analysis using a pre-ranked enrichment analysis revealed significant enrichment of dorsal striatum DE genes for several GO terms and Hallmark gene-sets. Genes in the CN were found to be related to cilia- and microtubule-associated GO-terms, while none of the Hallmark gene-sets was significantly enriched. GO-term and Hallmark gene-set analysis in PUT samples showed enrichment for immune processes, such as “acute inflammatory response to antigenic stimuli” (padj=0.006) and “adaptive immune response” (padj=0.006). In the VS the most significantly enriched GO-terms were also related to cilia and microtubules, as well as antigen processing. All GO-terms and Hallmark gene-set with FDR <0.10 are listed in Supplementary Tables S5 (CN), S6 (PUT), and S7 (VS).
WGCNA
Expression
In the CN, 21 modules with a median size of 352 genes (range: 64-7 259) were identified. Module “e-CN-magenta”, consisting of 328 genes, showed the strongest association with AUD status (r=0.42, p=2.89*10−4). In PUT module “e-PUT-black”, of the 25 modules (median size 249 genes, range: 33-5 381) identified, was most strongly correlated with AUD (r=0.41, p=2.31*10−4). For expression data from the ventral striatum, 16 modules with a median size of 429 genes (range: 35-9 708) were identified; module “e-VS-pink” had the strongest association with AUD (r=0.41, p=0.009). Interestingly, in a GO-term analysis the three AUD associated modules were all enriched for immune processes, such as “defense response” and “inflammation response”. There was also a wide overlap of the genes in the three modules: 174 (22.54%) were partially shared between all three modules corresponding to the three brain regions, while another 21.76% were shared between at least two modules (Figure 2B).
In the GWAS-enrichment analyses, the WGCNA modules “e-CN-magenta” and “e-PUT-black” were not enriched for signals from SUD-GWASs. “e-VS-pink” was enriched for genes associated with Cannabis Use Disorder (p=0.043).
Methylation
In the CN, WGCNA resulted in 36 modules with a median size of 346 CpG-sites (range: 66-41 423). Module “m-CN-red”, consisting of 2 117 CpG-sites, showed the strongest association with AUD case control status (r=-0.27, p=0.021). This module was most highly enriched for the biological processes “cell activation” (p=1.52*10−5) and “leukocyte activation” (p=2.09*10−5). In PUT 177 modules were identified (median size=57 CpG-sites, range: 30-42 248). Module “m-PUT-plum” consisted of 70 CpG-sites and was significantly associated with AUD case/control status (r=-0.29, p=0.023) and enriched for the biological processes “positive regulation of I-κB kinase/NF-κB signaling” (p=0.002) and “regulation of I-κB kinase/NF-κB signaling” (p=0.005). WGCNA in the VS methylation data resulted in 85 modules (median size=178 CpG-sites, range: 35-30 370). The module with the strongest association with AUD was “m-VS-lavender” (r=-0.29, p=0.023), which consisted of 117 CpG-sites and was enriched for the molecular function “natural killer cell lectin-like receptor binding” (p=3.43*10−4) and the biological process “susceptibility to natural killer cell mediated cytotoxicity” (p=3.65*10−4). The top 10 enriched GO-terms for all AUD-associated modules can be found in Supplementary Tables S8-S10.
Expression and Methylation Data Integration
In the CN, 12 genes showed both differential methylation and differential gene expression (see Supplementary Table S11). No overlap was observed in the VS and PUT.
At the module-level, co-expression module “e-CN-magenta” showed significant overlap with the methylation modules “m-CN-red” (p=0.003) and “m-CN-midnightblue” (p=0.014) (Figure 2C), while expression module “e-CN-purple” did not show significant overlap with the methylation modules in CN. Of the 3 AUD-associated expression modules in the VS, only “e-VS-salmon” showed significant overlap with the methylation module “m-VS-turquoise” (p=0.003), but not “m-VS-lavender”. No overlap was observed for gene expression and DNA-methylation in PUT.
GWAS Enrichment Analysis
In the VS, but not in the dorsal striatum, we observed enrichment of DE genes in GWAS signal of PAU. In the dorsal striatum DE genes were enriched for GWAS signal from a study comparing individuals with OUD to unexposed controls. None of the DE genes in any of the brain regions showed enrichment for signals from a GWAS of Cannabis Use Disorder or AUD. Results from the respective analyses are depicted in Figure 3.
From the WGCNA modules showing the strongest association with AUD, only module e-VS-pink showed significant enrichment for GWAS signals of CUD. Results of this analysis can be found in Supplementary Table S12.
Discussion
In the present study, we identified DE genes, co-expression networks, and pathways associated with AUD in the dorsal and ventral striatum. The results were integrated with DNA-methylation data and results from GWASs of SUDs.
We discovered that one gene (ARHGEF15) was consistently upregulated in all investigated brain regions of AUD cases compared to controls. ARHGEF15 encodes a specific guanine nucleotide exchange factor for the activation of Ras homolog family member A (RhoA), a GTPase, which has been linked to higher blood pressure and hypertension over the Rho/ROCK signaling cascade37. It is postulated that the Rho Guanine Nucleotide Exchange Factor 15 negatively regulates excitatory synapse development by suppressing the synapse-promoting activity of EPHB2. EPHB2 deficiency has been linked to depression-like behaviors and memory impairments in animal studies38. In line with this, genetic variation within ARHGEF15 has been associated with hematocrit, red blood cell count, and hemoglobin concentration39, but also with psychiatric traits, such as neuroticism and worries40 as well as bipolar disorder41.
Among the genes that were downregulated in both dorsal striatal regions, HLA-DOB displayed the highest fold change. HLAs of the Major Histocompatibility Class II are an essential part of the acquired immune system presenting antigens to T-lymphocytes (for review: Howell, Carter 42). The most significantly downregulated gene in the VS is CLPB, a mitochondrial chaperone, which has been associated with progressive brain atrophy43 and with the cellular response to alcohol-induced stress44. In a recent GWAS, CLPB was associated with the amount of alcohol consumed on a typical day (p=9.67*10−5, N=116 163)45.
DE genes in the ventral striatum were enriched for GWAS signals of PAU, but not AUD. This could be a result of the larger sample size of the PAU GWAS, but also point towards differences in genetic variation as responsible for differential expression.
Our results from the pathway and network analyses further underline immune-related effects of chronic alcohol exposure; the pathway and network modules most strongly associated with AUD case-control status were also enriched for immune system and inflammation processes. This was observed for all three brain regions, and both in expression and methylation data, providing further evidence for the important role of immune processes in AUD.
These results strongly reflect the well-described effect of chronic alcohol exposure on different aspects of the innate and acquired immune systems46. Chronic alcohol exposure accelerates the inflammatory response and reduces anti-inflammatory cytokines46. An activated immune response in response to chronic alcohol exposure has been shown on the cell level47, as well as on the transcription47, and protein levels48, 49. In a previous EWAS, we found strong enrichment of immune processes in differentially methylated CpG-sites associated with alcohol withdrawal50. Neuroinflammation has been repeatedly associated with AUD and both the glutamate excitotoxicity and the production of acetaldehyde, key processes in AUD metabolism, have been suggested to produce an inflammatory response in the brain51. On a phenotypic level, there is also widespread overlap between symptoms of inflammation and of SUDs, such as anhedonia, depression, and decreased cognitive functioning52. In addition, in postmortem human brain studies in the PFC, hippocampus, and orbitofrontal cortex, increased mRNA levels of HMGB1 encoding a proinflammatory cytokine and toll-like receptor genes have been associated with alcohol consumption in AUD cases, providing evidence for chronic neuroinflammation in response to alcohol53-55. Notably, there is an overlap of findings not only on the single-gene level but also on the level of pathways and networks/modules. This overlap underlines that alcohol consumption has common biological effects in different brain regions, i.e., most prominently, effects on immune and inflammation processes.
Several limitations apply to our study. First, we cannot distinguish between effects being a consequence of chronic alcohol consumption or addiction. Second, although we corrected for PMI, which can influence tissue quality as a confounding factor, it cannot be ruled out that other characteristics not easily accounted for, such as cause of death, or blood alcohol level for which the majority of individuals have missing data, influenced gene expression.
In summary, the present study provides further evidence from multi-omics data sets for the importance of immune-and inflammation-related processes in AUD. Notably, drugs that reduce neuroinflammation to reduce drinking, such as phosphodiesterases, may be promising approaches for novel treatment options for AUD. Recently published randomized controlled trials suggest that a phosphodiesterase inhibitor reduces heavy drinking whereas an antibiotic compound was not effective56, 57. A deeper understanding of the underlying mechanisms will enhance the discovery of drug targets and drive forward the development of precision medicine within in this field.
Data Availability
Raw data and summary statistics are available on request.
Conflict of Interest
The authors have nothing to disclose.
Supplementary Information
Supplementary information is available at MP’s website.
Data Availability
Raw data and summary statistics are available on request.
Acknowledgments
The study was supported by the German Federal Ministry of Education and Research (BMBF), “A systems-medicine approach towards distinct and shared resilience and pathological mechanisms of substance use disorders” (01ZX01909), “Towards Targeted Oxytocin Treatment in Alcohol Addiction” (031L0190). ERA-NET: Psi-Alc (01EW1908), and the Deutsche Forschungsgemeinschaft (DFG) – Project-ID 402170461 – TRR 26558. We thank Elisabeth Röbel and Claudia Schäfer-Arnold for technical assistance.
Tissues were received from the New South Wales Brain Tissue Resource Centre at the University of Sydney which is supported by the University of Sydney. Research reported in this publication was supported by the National Institute of Alcohol Abuse and Alcoholism of the National Institutes of Health under Award Number R28AA012725. The content is solely the responsibility of the authors and does not represent the official views of the National Institutes of Health.