Integrated Single-Cell Multiomic Profiling of Caudate Nucleus Suggests Key Mechanisms in Alcohol Use Disorder

Abstract

Alcohol use disorder (AUD) induces complex transcriptional and regulatory changes across multiple brain regions including the caudate nucleus, which remains understudied. Using paired single-nucleus RNA-seq and ATAC-seq on caudate samples from 143 human postmortem brains, including 74 with AUD, we identified 17 distinct cell types. We found that a significant portion of the alcohol-induced changes in gene expression occurred through altered chromatin accessibility. Notably, we identified novel transcriptional and chromatin accessibility differences in medium spiny neurons, impacting pathways such as RNA metabolism and immune response. A small cluster of D1/D2 hybrid neurons showed distinct differences, suggesting a unique role in AUD. Microglia exhibited distinct activation states deviating from classical M1/M2 designations, and astrocytes entered a reactive state partially regulated by JUND, affecting glutamatergic synapse pathways. Oligodendrocyte dysregulation, driven in part by OLIG2, was linked to demyelination and increased TGF-β1 signaling from microglia and astrocytes. We also observed increased microglia-astrocyte communication via the IL-1β pathway. Leveraging our multiomic data, we performed cell type-specific expression quantitative trait loci analysis, integrating that with public genome-wide association studies to identify AUD risk genes such as ADAL and PPP2R3C, providing a direct link between genetic variants, chromatin accessibility, and gene expression in AUD. These findings not only provide new insights into the genetic and cellular mechanisms in the caudate related to AUD but also demonstrate the broader utility of large-scale multiomic studies in uncovering complex gene regulation across diverse cell types, which has implications beyond the substance use field.

Extended Data Figure 1:. Proportion of each cell type for each of the 163 samples,

grouped by AUD classification, for all snRNA-seq barcodes used in cell clustering and cell type annotation (see Fig. 1).

Extended Data Figure 2:. D1 medium spiny neurons subtypes.

a, UMAP of D1 MSN cells, colored by compartment (either matrix or striosome). b, UMAP of D1 MSN cells, colored by expression of marker genes used to assign compartment. c, Dotpot of expression and prevalence of representative marker genes for matrix and striosome compartments.

Extended Data Figure 3.. D2 medium spiny neurons subtypes.

a, UMAP of D2 MSN cells, colored by compartment (either matrix or striosome). b, UMAP of D2 MSN cells, colored by expression of marker genes used to assign compartment. c, Dotpot of expression and prevalence of representative marker genes for matrix and striosome compartments.

Extended Data Figure 4:. Number of differentially expressed genes in multiple cell types.

Upset plot shows the number of genes with expression significantly associated with AUD (padj < 0.2) in different combinations of cell types.

Extended Data Figure 5:. snATAC-seq cell landscape.

UMAP plotting each cell for which snATAC-seq data was available, clustered based on snATAC-seq data. Cell type labels for each cell was provided based on the cell’s snRNA-seq data (see Fig. 2).

Extended Data Figure 6:. Pseudobulk accessibility profiles for each cell type at canonical marker genes.

5 kilobases on each side of the transcription start site of each gene are shown. Arrow indicates the direction of transcription.

Extended Data Figure 7:. Normalized pseudobulk expression of PPP2R3C (above) and ADAL (below) for each sample.

Plotted by genotype for the significant expression quantitative trait loci (eQTL) in the cell types in which the eQTL was significant.

Extended Data Figure 8:. Cell type-specific driver gene scores.

a, Dotplot of all driver genes found, using LINGER’s driver score, based on chromatin accessibility of target genes of that transcription factor. Size of dot corresponds to p-value and color indicates t-value, of change in driver score between AUD and control individuals. Left two columns correspond to membership in regulatory modules. Genes ordered by cell type with most significant difference in driver score. b, Dotplot of all driver genes found, using LINGER’s driver score, based on gene expression of target genes of that transcription factor. Size of dot corresponds to p-value and color indicates t-value, of change in driver score between AUD and control individuals. Left two columns correspond to membership in regulatory modules. Genes are ordered alphabetically.

Extended Data Figure 9:. Normalized pseudobulk expression of target genes with eQTLs.

Genes from modules 2 and 3 (as identified by LINGER) with a significant eQTL in astrocytes and oligodendrocytes. Plotted by genotype for the signficant eQTLs.

Extended Data Figure 10:. Downstream genes of AUD-associated ligand-receptor pairs.

a. Microglia to astrocytes ligand-receptor pairs from MultiNicheNet cell-cell communication. All ligand-receptor pairs and target genes with high expression correlation (Spearman or Pearson > 0.50), having some prior knowledge to support their link (in the top 250 predicted target genes for the ligand, ‘prior score’ as predicted by MultiNicheNet), and being within the top 50 ligand-receptor pairs associated with AUD (as calculated by MultiNicheNet) are shown. Size of dots indicate Pearson correlation between expression of ligand-receptor pair and target gene. Color of dot indicates prior score for link between ligand-receptor and downstream gene. b, As a, with astrocytes to oligodendrocytes pairs (above) and microglia to oligodendrocyte pairs (below)

Figure 1:. Experimental and bioinformatics workflow.

a, Sequencing experiments included single-nucleus RNA-seq (snRNA-seq) and single-nucleus multiome (sn-multiome) of 163 postmortem brain samples from individuals with and without alcohol use disorder (AUD). b, Bioinformatics approaches and tools used for analysis. See Online Methods for details.

Figure 2:. Cell type landscape of the caudate nucleus in alcohol use disorder.

a, UMAP plot of the 1,307,323 nuclei profiled in the snRNA-seq and sn-multiome assays; visualization shown is based on the snRNA-seq profile. Nuclei are labeled by cell type and cell type proportion among all snRNA-seq cells. Cell types: cholinergic neurons (Ach), astrocytes (Astro), cholecystokinin-expressing interneurons (CCK), calretinin-expressing interneurons (CR), D1-type medium spiny neurons (D1), D2-type medium spiny neuron (D2), medium spiny neurons expressing both D1 and D2 receptors (D1/D2), endothelial cells (Endo), ependymal cells (Epend), fast-spiking interneurons (FS), glutamatergic neurons (Glut), low-threshold spiking interneurons (LTS), non-microglial macrophages (Macro), microglia (Micro), oligodendrocytes (Oligo), oligodendrocyte progenitor cells (OPCs), and vascular smooth muscle cells (vSMCs). b, Normalized expression in each cell type of the marker genes used to identify cell types. Dot size corresponds to the percentage of cells expressing the gene; dot intensity indicates average gene expression level. c, Left, UMAP of microglial cells, colored by subcluster. Right, dotpot of expression and prevalence of representative marker genes for each microglial subcluster. d, Scatter plot showing the proportion of inflammatory microglia in each sample and the subject's age (red, AUD; blue, control). Bars on the left quantify the ratio of individuals with AUD to without AUD among those with >= 50% (above) or < 50% (below the solid black line) of microglia in the inflammatory state. e, Left, UMAP of astrocyte cells, colored by subcluster. Right, Dotpot of expression and prevalence of representative marker genes for each subcluster.

Figure 3:. Characterization of AUD-associated changes in gene expression in the caudate nucleus.

a, Barplot showing number of genes differentially expressed in individuals with AUD for the eight cell types which have over 100 differentially expressed genes. Red and blue indicate positively and negatively differentially expressed genes, respectively. See Fig. 2A caption for cell type abbreviations. b, Heatmap of Pearson correlation of gene expression changes (log2 fold changes) between cell types, hierarchically clustered by Pearson correlation. Black-outlined squares indicate groups of cell types with moderate correlation, namely, D1, D2, and D1/D2 neurons, and OPCs, astrocytes, and oligodendrocytes. c, Heatmap of biological pathways from the Reactome database enriched in brain samples from individuals with AUD in each cell type. The top 100 enriched pathways (based on smallest Benjamini-Hochberg adjusted p values) across all cell types are shown and were hierarchically clustered based on the number of genes shared between the pathways. Heatmap cell color indicates adjusted p value. Asterisk indicates negative enrichment score; all other pathways have positive enrichment scores. The pathways are divided into 25 clusters, which are manually labeled with a summary of pathways making up that cluster.

Figure 4:. Characterization of AUD-associated changes in chromatin accessibility in the caudate nucleus.

a, Venn diagram of overlap between the union of open chromatin regions from all neuronal cell types and the union of regions from all non-neuronal cell types. b, Heatmap of Jaccard similarity between open chromatin regions in each cell type, hierarchically clustered by Jaccard similarity. c, Barplot showing number of differentially accessible regions identified in oligodendrocytes, astrocytes, D1, and D2-type MSNs; red and blue indicate positively and negatively differentially accessible regions, respectively, and lighter and darker coloring indicate regions in promoter and enhancer regions of genes, respectively. Promoter regions were defined as 1 kilobase surrounding the transcription start side of each gene. d-g, Top, scatterplot of ATAC peak log(2) fold changes and RNA-seq log(2) fold changes for genes with at least one differentially accessible region (padj < 0.2). Genes are colored based on whether the gene is also differentially expressed (padj < 0.2). Bottom, GSEA enrichment plot of enrichment of the same ATAC-significant genes, split into two sets based on positive or negative effect size, across genes ranked by differential expression fold change. Normalized enrichment scores (NES) and Benjamini-Hochberg adjusted p values (padj) for each GSEA test are shown; d, oligodendrocytes; e, astrocytes; f, D1-type MSNs; g, D2-type MSNs.

Figure 5:. Integration of eQTL analysis with GWAS data and differential gene expression.

a, Overview of DEG-GWAS-eQTL integration. Potential driver genes in AUD are defined as genes that are differentially expressed, contain GWAS loci associated with an alcohol-related trait, and contain an eQTL. b, Upset plot showing number of differentially expressed genes (FDR < 0.2) containing a cis-eQTL (FDR < 0.2) and a GWAS locus for combinations of cell types. The phenotype trait of the GWAS loci overlapping the gene is indicated in color. PAU = Problematic alcohol use. c, For oligodendrocytes, differential expression log fold change (individuals with AUD vs those without) plotted against GWAS effect size (using the variant within the locus with the smallest p-value), multiplied by the eQTL effect size, for each gene with a significant eQTL (Benjamini-Hochberg adjusted p-value < 0.2). Shape indicates with which phenotypic trait the GWAS locus overlapping the gene is associated. d, e, f, as (c), for oligodendrocytes (d), D1 neurons (e), and D2 neurons (f).

Figure 6:. Cell type-specific gene regulatory networks associated with AUD.

a, Gene regulatory network analysis of samples from individuals with and without AUD comparing the average expression of genes within 10 regulatory modules (M1-M10; see Online Methods). Dot size represents Benjamini-Hochberg adjusted p-value and color represents t-value of the difference in average gene expression. Asterisk indicates significance (<0.05). Left two columns display enrichment of module genes in the set of genes significantly associated with problematic alcohol use (GWAS(PAU)) and drinks per week (GWAS(Drkwk)). b, Astrocyte co-expression modules. Left two columns are the same as (a). Top (red heatmap), Jaccard index between genes belonging to the 10 regulatory modules, and 10 co-expression modules (Co.E0-Co.E9), as calculated by WGCNA. Bottom (blue heatmap), Benjamini-Hochberg adjusted p-value of a t test of the difference in average expression of genes in each co-expression module, between individuals with and without AUD. DEG_1>0 indicates higher expression in samples from individuals with AUD, and DEG_0>1 indicates higher expression in those without AUD. Bold number “75” indicates number of genes overlapping between the regulatory module and co-expression module. c, Gene Ontology (GO) functional enrichment for the 75 genes overlapping regulatory module 2 and co-expression module 3 shown in (b). Numbers marked by asterisk indicates Benjamini-Hochberg adjusted adjusted p-value of enrichment. d, UMAP of astrocytes from individuals without AUD (left plot) and with AUD (right plot). Dot color indicates the enrichment of the JUND motif (red) and the log-normalized C3 expression (green). Yellow indicates high expression of C3 and high JUND motif enrichment. e, Boxplot of the log of chromvar motif activity score for the JUND motif in astrocytes for samples from individuals with and without AUD. f, Boxplot of the log of C3 expression in astrocytes for samples from individuals with and without AUD. g, Left, UMAP of oligodendrocytes, clustered and annotated into three subclusters using graph-based clustering. Right, dotplot of MBP and OLIG2 expression for each of the three oligodendrocyte subclusters. Dot size indicates percentage of cells expressing the gene, and dot color indicates average expression of the gene. h, Boxplot of the log of chromvar motif activity score for the OLIG2 motif in oligodendrocytes for samples from individuals with and without AUD. i, Circos plot showing the top five ligand-receptor interactions (determined by scaled ligand activity score from MultiNicheNet) between astrocytes, oligodendrocytes, and microglia.