Cell type-specific multi-omics analysis of cocaine use disorder in the human caudate nucleus

Abstract

Structural and functional alterations in the brain's reward circuitry are present in cocaine use disorder (CocUD), but their molecular underpinnings remain unclear. To investigate these mechanisms, we performed single-nuclei multiome profiling on postmortem caudate nucleus tissue from six individuals with CocUD and eight controls. We profiled 30,030 nuclei, identifying 13 cell types including D1- and D2-medium spiny neurons (MSNs) and glial cells. We observed 1485 differentially regulated genes and 10,342 differentially accessible peaks, with alterations in MSNs and astrocytes related to neurotransmitter activity and synapse organization. Gene regulatory network analysis identified transcription factors including ZEB1 as exhibiting distinct CocUD-specific subclusters, activating downstream expression of ion- and calcium-channels in MSNs. Further, PDE10A emerged as a potential drug target, showing conserved effects in a rat model. This study highlights cell type-specific molecular alterations in CocUD and provides targets for further investigation, demonstrating the value of multi-omics approaches in addiction research.

Fig. 1. Study workflow and characterization of single-nuclei multiome dataset.

A Sample workflow for N = 14 individuals, top left image was created with BrainNet Viewer (v.1.7; Xia et al., http://www.nitrc.org/projects/bnv/), and the bottom right image was created in BioRender. Espinola, M. (2025) https://BioRender.com/c45a312; and bioinformatics workflow relating each analysis to the corresponding figure. B Integrated UMAP for 30,030 nuclei, depicting 13 clearly delineated cell types. C UMAP of the RNA assay. D UMAP of the ATAC assay. E Cell type proportions per sample. F Dot plot depicting marker gene expression per cell type. G Coverage plot depicting chromatin peaks in promoter regions of cell type marker genes. CN = caudate nucleus, CocUD = cocaine use disorder, UMAP = uniform manifold approximation and projection.

Fig. 2. Differential expression and accessibility analysis of cocaine use disorder and GO overrepresenation.

A Heatmap depicting differentially expressed genes identified using a two-sided Wilcoxon rank sum test (Bonferroni-adjusted p < 0.05, log2FC > 0.5, expressed in at least 25% of cells), red = upregulated in CocUD, blue = downregulated in CocUD, top ten DE genes per cell type are labeled. Upset plot depicting the number of shared and distinct differentially expressed features per cell type, for B differential expression and C differential accessibility. D emapplot depicting the top 14 category results of the Gene Ontology overrepresentation analysis for DE genes in a cell type-specific manner. Each circle represents one GO term as a pie chart, indicating the strength of enrichment in the different cell types, represented by color. Edges indicate semantically similar GO terms, with overlapping genes. Circle size indicates the number of genes in the pathway. Overrepresentation was determined using a Fisher’s exact test and Benjamini-Hochberg adjustment for multiple testing. E emapplot depicting the top 14 category results of the Kyoto Encyclopedia of Genes and Genomes (KEGG) overrepresentation analysis for DE genes in a cell type-specific manner. Each circle represents one GO term as a pie chart, indicating the strength of enrichment in the different cell types, represented by color. Edges indicate semantically similar GO terms, with overlapping genes. Circle size indicates the number of genes in the pathway. F heatmap showing differentially accessible (DA) peaks identified using a two-sided Wilcoxon rank sum test (Bonferroni adjusted p < 0.05, log2FC > 0.25, expressed in at least 5% of cells), red = upregulated in CocUD, blue = downregulated in CocUD. G Annotation of DA peaks to their genomic location. H comparison between genomic location of DA peaks to the genomic background, using a two-sided chi-square test with Bonferroni adjustment for multiple testing, *** = p < 0.001, ** = p < 0.01, * = p < 0.05.

Fig. 3. Linked peaks and transcription factor motif and expression analysis.

A Heatmap depicting shared (brown), CocUD-specific (red), and control-specific (blue) links between chromatin accessibility and RNA expression by cell type. B Emapplot showing overrepresented gene ontology biological processes for genes characterized by a CocUD-specific linked peak in (A) color represents statistical significance and circle size the number of genes in the GO term; overrepresentation was determined using a Fisher’s exact test and Benjamini-Hochberg adjustment for multiple testing. C Heatmap depicting cell type-specific motif enrichment, yellow motif names indicate increased activity while motif names in blue represent reduced activity of the motif in cell types based on chromVAR activity scoring. D Feature plots confirming motif expression for MA1114.1/PBX3 in D1- and D2-MSNs, MA1577.1/TLX2 in astrocytes and MA0080.5/SPI1 in Microglia. Enrichment scores are color-coded from blue (no enrichment) to yellow (strong enrichment). E Heatmap depicting differential transcription factor motif activity between CocUD and controls, red = increased activity in CocUD, blue = decreased activity in CocUD. F Heatmap depicting corresponding transcription factor differential transcript expression in CocUD, red = upregulated in CocUD, blue = downregulated in CocUD. G Transcription factor motifs for ZEB1 and FOXP2, size of sequence represents probability. H Feature plots depicting ZEB1 and FOXP2 expression in cell types of the caudate nucleus (log-normalized counts).

Fig. 4. Gene regulatory network of D1- and D2-MSNs, candidate gene validation in the 3-CRIT model and sex-specific RRHO analysis.

A Gene Regulatory Network (GRN) from the D1- and D2-MSNs, color and size of nodes relate to centrality in the network, larger circles represent increased centrality, orange edges represent activation, and gray edges inhibition. B Feature plots depicting ZEB1 RNA expression (upper), module scores for ZEB1 suppressed downstream targets (middle), and module scores for ZEB1 activating downstream targets (lower), separately for CocUD (left) and controls (right). C Emapplot summarizing GO overrepresentation results from genes activated (red) and suppressed (turquoise) by ZEB1, PBX3, and ZNF148, the top TFs from regulon scoring. Each circle represents one GO term as a pie chart, indicating the strength of enrichment in the different regulons, represented by color. Edges indicate semantically similar GO terms, with overlapping genes. Circle size indicates the number of genes in the pathway. D Summary of validation analysis in 3-crit model of cocaine use disorder in rats, upper panel was Created in BioRender. Espinola, M. (2025) https://BioRender.com/c99g053, data are presented as mean values +/- SEM, light gray = controls (N = 6), gray = 0crit (N = 6), black = 3crit (N = 5), differential expression was determined using a Wilcoxon rank sum test, p_PDE10A = 0.015, p_CACNA1C = 0.026, p_KCNIP4 = 0.82, p_ZEB1 = 0.82. E RRHO plots showing the concordance between DE results in males and females, rank rank hypergeometric overlap tests were used to identify sex-specific concordance, upper left = upregulated in females and downregulated in males, upper right = upregulated in males and females, lower left = downregulated in males and females, lower right = upregulated in males and downregulated in females. MSNs = medium spiny neurons, GRN = gene regulatory network, CocUD = cocaine use disorder, Ctrl = control, (+) = activating, (-) = suppressing.