Multi-omics profiling of DNA methylation and gene expression alterations in human cocaine use disorder

Abstract

Structural and functional changes of the brain are assumed to contribute to excessive cocaine intake, craving, and relapse in cocaine use disorder (CUD). Epigenetic and transcriptional changes were hypothesized as a molecular basis for CUD-associated brain alterations. Here we performed a multi-omics study of CUD by integrating epigenome-wide methylomic (N = 42) and transcriptomic (N = 25) data from the same individuals using postmortem brain tissue of Brodmann Area 9 (BA9). Of the N = 1 057 differentially expressed genes (p < 0.05), one gene, ZFAND2A, was significantly upregulated in CUD at transcriptome-wide significance (q < 0.05). Differential alternative splicing (AS) analysis revealed N = 98 alternatively spliced transcripts enriched in axon and dendrite extension pathways. Strong convergent overlap in CUD-associated expression deregulation was found between our BA9 cohort and independent replication datasets. Epigenomic, transcriptomic, and AS changes in BA9 converged at two genes, ZBTB4 and INPP5E. In pathway analyses, synaptic signaling, neuron morphogenesis, and fatty acid metabolism emerged as the most prominently deregulated biological processes. Drug repositioning analysis revealed glucocorticoid receptor targeting drugs as most potent in reversing the CUD expression profile. Our study highlights the value of multi-omics approaches for an in-depth molecular characterization and provides insights into the relationship between CUD-associated epigenomic and transcriptomic signatures in the human prefrontal cortex.

Fig. 1. Differential expression analysis in CUD suggests synaptic signaling and immunological alterations in Brodmann Area 9.

a Volcano plot of the differential expression (DE) analysis showing the N = 378 upregulated (red) and N = 679 downregulated genes (blue) at nominal significance (p < 0.05). Solid black line indicates nominal significance (p < 0.05), dashed gray line indicates transcriptome-wide significance (FDR q < 0.05). b Results of the overlap analysis for upregulated (up) and downregulated (down) DEGs among cell type-specific marker genes. Green color depicts the odds ratio (OR) of overlap, p-values inside the panels indicates significance of overlap based on Fisher-Test. Gene-set enrichment analysis (GSEA) was performed for the DEGs in BA9 ranked by the Wald test statistic from DESeq2. Statistically significant results (q < 0.05) from GSEA separated by c negative and d positive normalized enrichment scores (NES) are shown in an enrichment map visualization. N.S. = not significant, OPC = oligodendrocyte progenitor cell.

Fig. 2. Differential alternatively spliced genes in CUD are related to neuron morphogenetic processes.

a Volcano plot of the differential alternative splicing (AS) results in Brodmann Area 9 (BA9). Statistically significant intron clusters (N = 108) identified by LeafCutter (q < 0.05) were annotated by gene name while dots represent individual introns of an intron cluster. Introns highlighted in red (dPSI>0) are more abundant in CUD while introns highlighted in blue (dPSI<0) are less abundant in CUD. b Results of the differential AS analysis at the cluster and gene level for one of the top findings, an intron cluster (clu_14172_-) in the Bridging Integrator 1 (BIN1) gene. Upper panel: visualization of BIN1 exons and introns with percent spliced in (PSI) measures related to the significant cluster clu_14172_-. The table indicates delta percent spliced in (dPSI) values from the CUD vs. Ctrl comparison. Lower panel: gene-level summary of all intron clusters detected in BIN1. FDR q-values are shown below cluster names. c GO enrichment analysis for the N = 98 AS genes harboring differentially statistically significant (q < 0.05) intron clusters in CUD. d Overlap of findings from differential expression (DE) analysis (N = 1057 DEGs at p < 0.05) and differential AS analysis.

Fig. 3. Replication analysis of CUD associated transcriptomic alterations in independent datasets.

a Overlap of nominally significant (p < 0.05) differentially expressed genes across datasets reveals two shared DEGs, HSPA6 and FKBP4. The replication datasets were based on N = 21 BA9 samples from the National PTSD Brain Bank (NPBB) and neuronal-specific transcriptomic data of N = 36 BA46 samples available under GEO accession number GSE99349 (BA46 replication). b HSPA6 is the strongest spliceosome-associated DEG in BA9. c Results of the look-up analysis for shared DEGs HSPA6 and FKBP4 - log2-fold change and p-value: association p-value from the DESeq2-based differential expression results. Significant associations are highlighted in bold. Rank-rank hypergeometric overlap (RRHO) visualization for d the BA9 replication dataset and e the neuronal-specific BA46 dataset indicating convergent and divergent expression patterns across studies using full differential expression statistics as the input datasets. Color scale represents -log₁₀(p) of the hypergeometric testing procedure in RRHO. Convergent expression across datasets: up-up and down-down, divergent: up-down and down-up.

Fig. 4. Convergence of DNA methylation, alternative splicing, and gene expression alterations in CUD at the ZBTB4 and INPP5E gene loci.

a Overlap of differential expression (DE), differential DNA methylation (DiffMeth), and differential alternative splicing analyses suggest two genes, ZBTB4 and INPP5E, where alterations are consistently associated with CUD. b Relationship between DE and DiffMeth genes in Brodmann Area 9 based on log2FC (y-axis) from DE analysis and effect size β from linear regression in the EWAS of CUD (x-axis). For both genes, ZBTB4 and INPP5E (highlighted in red), hypomethylation of the strongest significant CUD-associated CpG site and increased transcript levels are observed. c Integrated visualization of functional genomic datasets for ZBTB4 and INPP5E gene loci. CUD-associated genomic variants (SNPs p < 0.05 from [49]), CUD-associated CpG sites (p < 0.05 from [22]), RNA-seq data and intron clusters (q < 0.05) from the present study were visualized together with ENCODE ChIP-seq data for different chromatin marks in human dorsolateral prefrontal cortex.