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

Cocaine use disorder (CUD) is a prevalent substance use disorder characterized by compulsive cocaine use, craving, and high relapse rates. There is no FDA-approved pharmacotherapy, and treatment focuses on symptom reduction. Neuroimaging shows structural and functional alterations in frontal cortical areas, including the prefrontal cortex (PFC), which are implicated in inhibitory control and addiction neurocircuitry. Epigenetic and gene expression changes have been proposed as molecular bases for CUD-associated brain alterations. The study aims to characterize molecular underpinnings of CUD in human PFC (Brodmann Area 9) through integrated multi-omics analyses (DNA methylation, gene expression, and alternative splicing), replicate findings in independent datasets, contextualize results biologically, and explore drug repositioning options to reverse CUD-associated expression profiles.

Prior work in animal models identified differential DNA methylation in PFC and nucleus accumbens and gene expression changes involving immediate early genes (Egr1, Nr4a1, Fos), neurotransmission, ion channels, and metabolism. Human postmortem studies in striatal regions reported CUD-associated differentially methylated regions in nucleus accumbens and caudate nucleus (RRBS), and transcriptome-wide upregulation of synaptic transporter genes with downregulated immune processes. In human PFC, neuronal nuclei from BA46 showed 883 nominally significant DEGs with co-expression networks enriched for GTPase signaling and neurotransmitter secretion. The authors previously identified 20 CUD-associated DMRs and co-methylation modules enriched for synaptic signaling in BA9. Alternative splicing alterations have been reported in AUD and OUD human brain tissues and in mouse cocaine models, but not systematically in human CUD. The current study addresses the gap by integrating DNAm, gene expression, and alternative splicing in the same BA9 samples with replication.

- Cohort and tissue: Postmortem human BA9 tissue from Douglas Bell Canada Brain Bank (DBCBB). Inclusion: age >18, DSM-IV cocaine dependence (termed CUD), male, European American. Excluded: severe neurodevelopmental/psychiatric disorders (except depressive disorders), other SUDs (except AUD). N=42 total for DNAm; N=25 (RIN >5.5) for RNA-seq. Demographics assessed; groups balanced on pH, PMI, RIN, comorbid MDD/AUD; suicide more frequent in CUD. - DNA methylation: Illumina EPIC (850k) array on N=42. QC via CPACOR-based pipeline; neuronal fraction estimated (Houseman). Quantile normalized beta values, logit-transformed to M values. EWAS with linear regression adjusting for age, PMI, pH, neuronal fraction, MDD/AUD comorbidity, technical factors. Downstream: DMRs, GO enrichment (CpG Passoc<0.001), WGCNA co-methylation. - RNA-seq: Total RNA from ~5 mg tissue; rRNA depletion; NEBNext Ultra II Directional library prep; NovaSeq 6000, ~60M read pairs (2x100bp). QC FastQC; alignment STAR (GRCh38); quantification featureCounts (GENCODE v43). DESeq2 differential expression with model mRNA ~ CUD + age + PMI + pH + RIN. Nominal p<0.05; FDR q<0.05 transcriptome-wide significance. Absolute log2FC>0.07 threshold (~5% change). VariancePartition confirmed covariates; MDD/AUD not included due to minimal variance; sensitivity analysis including them showed stable results. - Cell type deconvolution: CIBERSORT using curated PFC marker genes (10x specificity threshold). Bayesian estimation (BEST) tested differences; none significant. Overlap analysis of DEGs with cell-type markers via GeneOverlap (Fisher test). - Functional enrichment: GSEA for GO terms using DESeq2 Wald statistic; FDR q<0.05; emapplot visualization. - WGCNA: Networks on vst-normalized counts; soft power 9; 27 modules detected. Module-trait correlations computed; GO enrichment for significant modules; hub genes via gene significance x module membership; PPI networks in Cytoscape/stringApp. - Replication: Two datasets: (1) BA9 bulk RNA-seq (NPBB; N=7 CUD, N=14 controls), DESeq2 with age, sex, PMI, RIN as covariates. (2) Neuronal nuclei BA46 RNA-seq (GSE99349; N=19 CUD, N=17 controls); processed identically; model with age, RIN, pH, PMI. Overlap of nominal DEGs (p<0.05), targeted look-up for key genes (e.g., ZBTB4, INPP5E), and RRHO2 for transcriptome-wide concordance. - Drug repositioning: Connectivity Map (CMap) query with top 150 up- and downregulated genes (by DESeq2 statistic). Extract normalized connectivity scores (NCS), FDR-adjusted p-values; assess perturbagen classes and targets. - Differential splicing: LeafCutter (annotation-free). Two-pass STAR mapping with junction augmentation; intron clustering; Dirichlet-Multinomial GLM including age, PMI, pH, RIN. Significance: |dPSI|>0.025 and FDR q<0.05. Visualization with leafviz; GO enrichment (clusterProfiler). - Integrative locus analysis: Spark visualization combining GWAS summary stats for cocaine dependence (EA N=6378), EWAS (from prior study), RNA-seq, splicing, and ENCODE dIPFC male ChIP-seq marks (H3K4me1/3, H3K27ac, H3K36me3, H3K9me3, H3K27me3). MAGMA gene-based analysis for GWAS signal. - MOFA: Multi-omics factor analysis on matched N=25 with top 20,000 variable promoter CpGs (TSS200/TSS1500) and top 20,000 variable genes. 10 factors, convergence_mode="slow", seed=42, maxiter=10000. Correlate factors with covariates; GSEA on positive/negative weights; missMethyl enrichment for CpG weights. - Integrative GO analysis: compareCluster across 10 CUD-associated gene/CpG lists from EWAS, WGCNA (DNAm and expression), DEGs, AS genes, MOFA weights; FDR q<0.05 pathways mapped into functional modules.

- Differential expression in BA9: 1,057 DEGs at p<0.05 (378 up, 679 down). ZFAND2A significantly upregulated at FDR q=0.04 (log2FC=0.43, p=1.98e-06). Upregulated DEGs enriched for neuronal markers; downregulated DEGs enriched for astrocyte, endothelial, and oligodendrocyte markers. - GSEA: 276 GO terms with positive NES and 782 with negative NES at FDR q<0.05. Top positive: vesicle-mediated transport in synapse (NES=2.63, q=5.31e-15). Top negative: superoxide metabolic process (NES=-2.43, q=2.08e-06). Major clusters: neurotransmission/synaptic signaling; transmembrane transporter activity; inflammatory/immune processes; angiogenesis; ECM organization; gliogenesis. - WGCNA: Module yellow (2,517 genes) negatively correlated with CUD (r=-0.47, p=0.02); enriched for small molecule and carboxylic acid catabolism with a prominent fatty acid metabolism cluster. PPI hub nodes: APOE (9 edges), ERBB2 (8), ALDH7A1 (7), PPARA (7), TLR4 (7). - Alternative splicing: 108 significant intron clusters across 98 genes (FDR q<0.05). Top gene: BIN1 (q=7.8e-04). Enriched GO terms include cell junction assembly and neuron projection extension (both q=3.62e-03). Eight genes overlapped between AS and DEGs: ITPKB, CPLX1, HLA-F, INPP5E, GALNT8, IGFBP6, ZBTB4, BCAT2. - Replication: Two shared upregulated DEGs across BA9 discovery, BA9 replication, and BA46 replication: HSPA6 and FKBP4, with consistent effect sizes (e.g., HSPA6 log2FC: 2.59/2.45/3.90; p=0.002/0.005/0.02). RRHO showed strong convergence between BA9 discovery and BA9 replication; selective convergence for shared upregulated genes with BA46. - Drug repositioning: Glucocorticoid receptor agonist medrysone showed strongest negative connectivity (NCS=-1.78, q=2.2e-16). Glucocorticoid receptor-targeting drugs (agonists and antagonists) overall had significant negative connectivity scores, suggesting potential to reverse CUD expression signatures. The GO term response to glucocorticoid had negative NES and was FDR-significant (NES=-1.54, q=0.019). - Gene-level multi-omics convergence: ZBTB4 and INPP5E showed hypomethylated CpGs, increased expression (ZBTB4 log2FC=0.08, p=0.015; INPP5E log2FC=0.17, p=0.025), and significant AS clusters (ZBTB4 chr17:clu_10246_ q=0.028; INPP5E chr9:clu_25078- q=0.015). Strongest CpG associations: ZBTB4 cg03443505 (β=0.84, p=1.01e-05); INPP5E cg18558462 (β=0.93, p=8.55e-03). MAGMA indicated a stronger gene-based GWAS signal for ZBTB4 (Z=1.74, p=0.04) than INPP5E (Z=-1.14, p=0.87). ENCODE marks indicated active chromatin at both loci. - MOFA: Factor 9 correlated with CUD (r=-0.48, p=0.02) and age (r=0.47, p=0.02); lower factor values in CUD. Top weights included CpGs in MTA3 (positive) and TES (negative); genes RAB6A (positive) and HIVEP2 (negative). GSEA on negative expression weights enriched for synaptic signaling, cell junction organization, neurogenesis; positive weights enriched for cellular respiration and small molecule metabolism. DNAm weight enrichments suggested intracellular calcium regulation, synaptic vesicle processes, lactate transport, and ER stress. - Integrative GO modules: Convergent modules included neurotransmission/synaptic signaling (FM1), neuron/glial differentiation and morphogenesis (FM2), synapse and cell junction organization (FM3), and fatty acid metabolism (FM4).

The integrated multi-omics approach in BA9 reveals that CUD is associated with coordinated alterations in DNA methylation, gene expression, and alternative splicing, particularly impacting synaptic signaling, neuronal morphology, and metabolic pathways. The transcriptome identifies ZFAND2A as a robust upregulated gene, and replication across independent PFC datasets highlights HSPA6 and FKBP4 as conserved DEGs, implicating heat shock/spliceosomal machinery and glucocorticoid receptor signaling. The convergence of DNAm, expression, and splicing at ZBTB4 and INPP5E suggests potential regulatory nodes: ZBTB4 as a methyl-DNA-binding transcriptional repressor with possible downstream effects on neuronal endocytosis pathways, and INPP5E as a cilia-specific phosphatase potentially influencing GPCR trafficking and ciliary signaling relevant to psychostimulant responses. WGCNA and MOFA point to fatty acid metabolism and electron transport alterations alongside neuronal morphogenesis and cell junction organization, consistent with observed neuroplastic and metabolic changes in addiction models. Drug repositioning findings suggest that modulating glucocorticoid receptor signaling could reverse CUD-associated expression profiles, aligning with preclinical data showing complex roles for glucocorticoid agonists and antagonists in cocaine-related behaviors. Overall, the findings support the hypothesis that epigenetic regulation and transcriptional reprogramming contribute to PFC dysfunction in CUD and identify candidate therapeutic targets and pathways.

This study provides an in-depth molecular characterization of CUD in human BA9 by integrating DNA methylation, gene expression, and alternative splicing data from the same individuals. It identifies ZFAND2A as an FDR-significant upregulated gene, replicates HSPA6 and FKBP4 across independent PFC datasets, and highlights ZBTB4 and INPP5E as convergent multi-omics loci in CUD. Pathway-level integration implicates synaptic signaling, neuronal morphogenesis, cell junction organization, and fatty acid metabolism as core biological processes. Connectivity mapping suggests glucocorticoid receptor-targeting drugs as candidates to reverse CUD-associated transcriptional changes. Future research should include larger, more diverse cohorts; longitudinal and cell-type-specific analyses; functional validation including proteomics; and preclinical/clinical studies testing modulation of glucocorticoid signaling and metabolic pathways.

- Cross-sectional postmortem design captures end-stage changes, limiting inference on temporal dynamics of DNAm, expression, and splicing. - Modest sample size with few FDR-significant DEGs; limited statistical power reduces generalizability. - Cohort homogeneity (male, European American) improves internal validity but limits assessment of sex- or ancestry-specific effects. - Important covariates (e.g., cocaine at death, cause of death) could not be modeled due to power constraints; suicide differed between groups. - Bulk tissue analyses may miss cell-type-specific effects despite deconvolution showing no major composition differences. - mRNA changes may not reflect protein-level alterations; proteomic validation is needed.