Cell type-specific Multi-Omics Analysis of Cocaine Use Disorder in the Human Caudate Nucleus

doi:10.21203/rs.3.rs-4834308/v1

Download PDF

Article

Cell type-specific Multi-Omics Analysis of Cocaine Use Disorder in the Human Caudate Nucleus

https://doi.org/10.21203/rs.3.rs-4834308/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

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 31,178 nuclei, identifying 13 cell types including D1- and D2-medium spiny neurons (MSNs) and glial cells. We observed 1,383 differentially regulated genes and 10,235 differentially accessible peaks, with alterations in MSNs and astrocytes related to neurotransmitter activity and synapse organization. Gene regulatory network analysis identified the transcription factor 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.

Biological sciences/Computational biology and bioinformatics/Gene regulatory networks

Biological sciences/Neuroscience/Molecular neuroscience

Health sciences/Diseases/Psychiatric disorders/Addiction

Cocaine is one of the most frequently consumed illicit psychostimulants worldwide and poses a significant public health challenge. In the US alone, about 4.8 million people, or 1.7% of the population, used cocaine in the past 12 months and 0.5% are diagnosed with a cocaine use disorder (CocUD)¹. Cocaine use is associated with a range of negative side effects, most prominently cardiovascular and respiratory complications². In recent years, there has been a rise in deaths caused by cocaine overdoses, especially in the African American population³.

A complex interplay between genetic and environmental factors is assumed to contribute to an individual’s susceptibility and resilience to CocUD⁴. In recent years, epigenetics has been increasingly studied as an important mechanism for how environmental factors influence the transcription of genetic information, thus altering brain function⁵. However, the role of epigenetic changes and their relevance for altering gene expression levels in brain regions implicated in CocUD remains unclear⁶, requiring joint profiling of epigenetics and gene expression and their integration using multi-omics analyses.

Our current understanding of the neurobiological basis of CUD implicates the involvement of a number of different brain regions involved in the development and maintenance of this disorder. Most notable are the ventral striatum, which is involved in reward processing as part of the mesocorticolimbic pathway⁷, and the dorsal striatum – consisting of the putamen and caudate nucleus – which plays a critical role in habit formation and decision-making^8,9. Most of our current knowledge comes from preclinical studies, with the majority investigating gene expression in the ventral striatum, particularly the nucleus accumbens. To date, only a few studies have focused on human postmortem brain tissue. An early study using microarray technology identified a reduction in myelin-related genes in the nucleus accumbens in CUD¹⁰. Studies in other brain regions, such as the hippocampus, highlight the relevance of metabolic alterations, particularly in mitochondria, in CocUD¹¹, as well as changes in ion- and potassium-channel expression and receptor density¹². While studies in the human striatum are sparse, more evidence is available from rodent models. On the bulk level, several studies highlight CREB-signaling^13,14, immediate early genes (IEGs)^13-15, and ion-channel expression^16-18. Interestingly, a recent study showed that transcriptomic signatures associated with CocUD widely overlap in the nucleus accumbens and caudate nucleus in postmortem human brain tissue¹⁸. Furthermore, differential gene expression in human postmortem brain tissue was concordant with a mouse model, highlighting the translational value of such model systems in substance use disorders¹⁸.

One major limitation of bulk-level studies results from the lack of cell type specificity of gene expression alterations. So far, cell type-specific gene expression studies have only been performed in the rat nucleus accumbens and drosophila models of CocUD. These studies suggested the importance of D1-medium spiny neurons (D1-MSNs)¹⁹ and confirmed cell type-specific alterations in IEG expression²⁰.

A second research gap depicts the lack of integrative multi-omics analyses in the same brain tissue specimens at single-cell resolution. Multi-omics analysis allows for a more mechanistic investigation of molecular alterations by combining evidence from multiple layers. With the 10x Genomics Multiome technology gene expression and open chromatin states can be profiled from the same cell, providing a powerful tool for linking epigenetic mechanisms to gene expression changes in CocUD. Here, we can directly infer transcription factor activity and construct gene regulatory networks to identify upstream regulators of transcriptional changes.

In our study, we aim to address these two research gaps by profiling the differential gene expression patterns and chromatin accessibility landscapes in a cell type-specific manner within the caudate nucleus of individuals with CocUD compared to non-affected subjects. By employing state-of-the-art sequencing techniques and data analysis, we sought to identify key upstream regulators of transcriptional changes that underlie the pathophysiology of CocUD and validate them in the well-established 3-crit model of CocUD in rats^21,22. To our knowledge, this is the first study that has profiled human CocUD at the single-cell multi-omics level, in a highly affected but underrepresented population.

Identification of cell types in the human caudate nucleus

We isolated nuclei from the human caudate nucleus of six individuals with CocUD (mean age 55) and eight controls (mean age 57), all of African American ancestry (Fig. 1A). CocUD cases and controls were balanced according to age and sex, see also Sup. Table S1. We performed single-nuclei ATAC-seq and RNA-seq using the 10x Genomics Multiome platform, which allows for simultaneous profiling of chromatin accessibility and gene expression from each cell. After quality control, we retained 31,178 nuclei with a median of 1,641 expressed genes and 14,418 ATAC fragments per nucleus. For normalization and dimension reduction, we used Seurat²³ and Signac²⁴. The R package Harmony²⁵ was used for the integration and reduction of batch effects (Fig. 1B-D, Sup. Fig. S1A-S1E). Weighted nearest neighbor analysis resulted in an integrated UMAP depicting 13 clearly delineated cell types, including astrocytes, microglia, oligodendrocytes, oligodendrocyte precursor cells (OPC), endothelial cells and lymphocytes, as well as different GABAergic neurons, such as D1 and D2 medium spiny neurons (MSNs) and a small population of cholinergic neurons (Fig. 1B). Distribution of cell types was comparable between samples, with one control sample having a higher ratio of astrocytes to oligodendrocytes (Fig. 1E). Annotation of clusters was supported by the expression of known marker genes (Fig. 1F), as well as their promoter accessibility (Fig. 1G). As expected in this brain region²⁶, the main neuronal cell types observed were GABAergic D1 and D2 medium spiny neurons (MSN).

Differential Expression & Accessibility Analysis

We performed within cell type differential expression analysis comparing CocUD cases with controls. Here, we identified 1,383 differentially expressed (DE) genes (Sup. Table S2), clustering distinctly according to glial or neuronal cell identity (Fig. 2A). Strongest deregulation of gene expression levels was observed in astrocytes (305 DE genes) and medium spiny neurons (D1: 27, D2: 70, D1/D2: 224 DE genes). The top upregulated DE gene in astrocytes was MBP (log2FC=1.76, p<2.2e^-16), the myelin basic protein involved in myelination of the nervous system, while the top downregulated gene was LSAMP (log2FC=-0.78, p<2.2e^-16), which encodes the limbic system associated preprotein involved in axon guidance and neuronal growth. In both D1- and D2-MSNs, the top upregulated gene was the ribosomal gene RPS27A (log2FC_D1=2.09, p_D1<2.2e^-16;log2FC_D2=2.32, p_D2<2.2e^-16) and the top downregulated gene SLC35F3, a thiamine transporter (log2FC_D1=-1.76, p_D1<2.2e^-16; log2FC_D1=-1.04, p_D1<2.2e^-16). Interestingly, we observed an overall large overlap between DE genes in D1- and D2-MSNs (Fig. 2B). Gene Ontology (GO) enrichment analysis revealed an overrepresentation of DE genes from all clusters in ribosomal pathways, further a neuronal-specific cluster related to the respiratory chain and a synaptic cluster evident in astrocytes, microglia, D1- and D2-MSNs (Fig. 2C). In a subsequent KEGG analysis, we observed enrichment for addiction-related pathways, such as morphine addiction, and pathways related to neurodegenerative diseases (Fig. 2D).

Differential accessibility analysis revealed 10,235 differentially accessible (DA) peaks (Sup. Table S3), again with a clear distinction between glial and neuronal cell types (Fig. 2E). DA peaks were overrepresented in promoter and exonic regions, compared to the genomic background (Fig. 2F & 2G). In a GO term enrichment analysis performed separately for promoter and distal peaks, we observed an overrepresentation of voltage-gated ion channel and serine-threonine kinase activity in promoter peaks (Sup. Fig. S2A) and ligand- and voltage-gated channel and methyltransferase activity in distal peaks (Sup. Fig. S2B).

CocUD-specific peak-linked genes are associated with ribosomal activity, synaptic pathways, and nervous system development

To identify CocUD-specific linked peaks, we linked peaks to gene expression separately for samples from individuals with and without CocUD, revealing N=821 CocUD-specific, N=3,624 common, and N=1,159 control-specific gene-peak pairs (Fig. 3A). As multiple peaks can be linked to a single gene, we then further investigated the functional relevance of linked genes. GO term analysis of CocUD-specific peak-linked genes showed overrepresentation of ribosomal and synaptic pathways, as well as a pathway cluster related to nervous system development (Fig. 3B).

Family-wise motif enrichment of transcription factors in cell types

Altered open chromatin states in CocUD can increase or decrease the availability of transcription factor motifs thereby resulting in transcription factor-mediated downstream effects on gene expression. To evaluate differential chromatin accessibility related to transcription factor motifs, we performed a cell type-specific motif enrichment analysis (Fig. 3C) suggesting cell type-specific motif enrichment in CocUD (Fig. 3D). Here, we observed major downregulation of transcription factor (TF) motifs in CocUD of immediate early genes, such as JUND and FOS, in D1- and D2-MSNs, while the same motifs were upregulated in astrocytes, oligodendrocytes and OPCs (Fig. 3E and Sup. Table S4). Strong effects were also present for the SOX and RFX families. To determine whether there was differential TF expression, alongside differential motif availability, we specifically investigated differential TF expression (Fig. 3F). Here, the strongest differential expression patterns were observed for ZEB1, RFX3, FOXP2, and NFIB. For ZEB1, we observed reduced availability in CocUD of its motif (Fig. 3G) in neurons while increased availability in astrocytes was found. At the same time, statistically significant downregulation of the ZEB1 transcript was found across neuronal and astrocyte clusters in CocUD. We confirmed robust expression of ZEB1 in neuronal and astrocyte clusters across conditions (Fig. 3H) implicating it is a promising candidate transcription factor for inducing downstream differential expression in these cell types. Another deregulated transcription factor in CocUD was FOXP2, for which we found altered motif availability and transcript expression especially in D1-MSNs and microglia (Fig. 3G+3H).

Gene Regulatory Networks reveal TFs accounting for transcriptional changes in CocUD

Following up on the observation of common transcriptional changes in D1- and D2-MSNs and the identification of transcription factor candidates such as ZEB1 and FOXP2, we were interested in gene regulatory networks (GRN) related to CocUD. GRNs can be used to infer the relationship between TFs and their downstream targets (Fig. 4A). An MSN-specific GRN was constructed using pando²⁷ and metrics are depicted in Sup. Table S5 and Sup. Fig. S3A. As immediate early genes are known to influence gene expression of a large set of downstream targets after exposure to multiple drugs of abuse, we were interested if there are conserved sets of TFs and their targets in our CocUD dataset that might be of particular importance for the induction of altered transcriptional programs. Appropriately, two subclusters emerged in the GRN, one consisting mainly of immediate early genes, such as FOSB or JUND, targeting heat shock proteins (e.g. HSPA1B, HSPA1B) and downstream map kinases (Fig. 4A). A larger cluster encompasses the TFs highlighted in the motif analysis, targeting calcium-, potassium- and ion-channels. To identify the most relevant TFs, we performed a scoring approach wherein the size of the TF regulon is weighed by the log2FC from a differential regulon expression analysis (Sup. Fig. S3B & Sup. Table S6). The regulon refers to all downstream targets that are either activated or suppressed by the associated TF (Sup. Fig. S3C). Hence, one TF can have two regulons for activated downstream genes, e.g. ZEB1(+), and for suppressed ones, e.g. ZEB1(-). By investigating the regulons, we observed that a small number of TFs account for a large number of genes, many of which have been previously identified in the DE analysis (Sup. Fig. S3D). A TF with particularly large effects in both activating and suppressing downstream gene expression was again ZEB1 that is widely expressed in D1- and D2-MSNs (Sup. Fig. S3B & S3C; Sup. Table S6). Interestingly, we observed that ZEB1 has specific suppressing effects in controls and activating effects in individuals with CocUD: CocUD-specific subclusters of MSNs express downstream targets activated by ZEB1, while suppressed ZEB1 targets are preferably expressed in control individuals (Fig. 4B). GO ORA analysis of the downstream targets of the three most important TFs derived from the GRN, ZEB1, PBX3, and ZNF148, revealed common suppressing effects on genes related to ribosomal pathways while activated target genes were consistently involved in synaptic signaling pathways (Fig. 4C). To further explore the GRN in MSNs, we characterized the chromatin structure of highly connected TF targets, such as PDE10A and CACNA1C by visualizing differential accessibility in TF binding sites between individuals with CocUD and controls (Sup. Fig. S3E & S3F).

Druggability analysis suggests phosphodiesterase-inhibitors as potential pharmacological targets

When identifying potential pharmacological targets for disorders with complex molecular alterations, it is crucial to preserve this complexity in druggability analysis. Therefore, we performed a druggability analysis of highly connected target genes from the D1-/D2-MSN GRN, using the drug-gene interaction database (DGldb). Here, we observed several pharmacological agents with the potential to treat CocUD (Sup. Table S7). Interestingly, among the top drug targets with large effect size in the DE analysis, we observed PDEs (PDE10A and PDE7B) several times, suggesting brain-expressed phosphodiesterases as a potential treatment target in CocUD.

Validation of differential expression in PDE10A and CACNA1C in the 3-crit model of CocUD

We prioritized the DE genes from the CocUD DE analyses based on their relationship with the DA peaks (CACNA1C), their importance in GO terms (KCNIP4), their role in the GRN (ZEB1) and in the druggability analysis (PDE10A). To validate these findings, we performed RT-qPCR analysis on mRNA extracted from the caudate nucleus of rats that underwent the 3-crit model of CocUD, a cocaine self-administration procedure after which three addiction-like criteria are assessed: persistence of drug seeking, motivation for drug taking, and resistance to punishment. Animals fulfilling all three criteria are classified as “addicted-like”. We compared five 3-crit animals, to six 0-crit animals that were exposed to cocaine, but did not develop addiction-like symptoms and six cocaine-naïve control animals (N_total= 17, Fig. 4D). We observed differential expression in PDE10A (t=2.89, p=0.022), consistent with an upregulation in the non-neuronal cell types (Sup. Table S2), and in CACNA1C (t=-2.75, p=0.024), in line with a downregulation in MSNs (not significant in non-neuronal cell types, Table S2). These effects were significant for 0-crit, but not 3-crit animals, compared to controls, indicating either a cocaine-consumption effect or a habituation in the 3-crit animals. No effect was observed for ZEB1 (t=0.15, p=0.88) and KCNIP4 (t=0.05, p=0.96).

Multi-Omics Factor Analysis reveals a sex-specific factor characterized by metabolic pathway enrichment

In addition to the data-informed approaches, we performed a data-driven multi-omics factor analysis to identify a dimensionality-reduced set of factors associated with CocUD. Of the 15 resulting factors, four were associated with CocUD (Sup. Fig. S4A and S4B). Factor6 is characterized by high factor loadings for microglia, with positive enrichment for immune pathways and negative enrichment for pathways related to neuronal development (Sup. Fig. S4B & S4C).

Factor7 was mainly associated with GABAergic neurons and enriched for GO terms related to metabolic processes and synaptic signaling (Sup. Fig. S4D). Overall, the results for the CocUD-associated factors confirmed the findings of the DE analysis.

Of particular interest is the CocUD-associated Factor8 that indicated sex-specific subclusters in the majority of cell types, marked by high factor-loadings in subclusters of OPCs, oligodendrocytes, MSNs, and astrocytes. Characterized by an enrichment for metabolic GO terms, this factor suggests a male-specific effect. Interestingly, in the promoter-associated motifs of Factor8, we observe a negative enrichment of major regulators from the GRN, such as ZEB1 and TCF4 (Sup. Fig. S4E).

Sex-stratified DE analysis shows concordant effects in astrocytes, microglia, and D1-/D2- MSNs

Because we observed high Factor8 loadings in the majority of cell types, we performed exploratory analyses comparing the effects in the DE analysis between male and female participants using rank-rank hypergeometric overlap (RRHO) analysis. Interestingly, we observed high concordance of DE genes between sexes in astrocytes, GABAergic D1MSN, and microglia, and moderate concordance between GABAergic D2 MSNs, and OPCs (Fig. 4E). At the same time, strong sex-specific effects were evident in GABAergic CALB2 MSNs and oligodendrocytes, underlining the results from MOFA.

No cell type-specific enrichment of GWAS signals of CocUD

To test whether GWAS signals of substance use disorders and CocUD in particular, were overrepresented in cell type gene expression patterns, we used single-cell disease-relevance scoring (scDRS). While we observed low enrichment scores of CocUD in microglia and cholinergic neurons, feature plots did not support a clear cell-type specific effect (Sup. Fig. S4F, Sup. Table S8).

Our single-nuclei multi-omics study of CocUD in the human caudate nucleus reveals first insights into cell type-specific alterations of gene expression and the underlying chromatin states, directly linking transcription factor activity to downstream gene expression. To our knowledge, this is the first study that investigated human CocUD at single-nuclei resolution combining ATAC- and RNA-seq in a multi-omics approach in the same cell. Our study highlights synaptic pathways, ribosomal activity, and metabolic pathways as altered in CocUD, often in a cell type-specific manner.

We present different lines of analysis from data-informed differential gene expression and DNA accessibility analyses to data-driven GRN and factor analyses, all with convergent results.

In both the DE and DA analyses, we observed initial evidence for altered ion- and calcium-channel activity, implicating altered synaptic signaling in CocUD. In addition, we found an overrepresentation of genes in pathways related to synapse organization in CocUD. Changes in synaptic plasticity following repeated cocaine exposure have been well described ^28,29. This is also in line with the large effects observed for ribosomal genes, both in the DE, DA and GO ORA analyses, as the transcription of ribosomal DNA is an important prerequisite to neural plasticity^30,31. Interestingly, we also observed KEGG-enrichment of neurodegenerative diseases and addiction-related processes in the neuronal cell types, providing evidence that our results reflect known pathologies. In addition, we observed vast expression differences in astrocytes, related to synaptic processes. Consistent with our results, it has been shown that cocaine can trigger astrocyte-mediated synaptogenesis in the nucleus accumbens shell ³².

To exploit the strength of the multiome framework, the simultaneous profiling of RNA and open chromatin in the same cell, we calculated a GRN in D1- and D2-MSNs. Here, we observed two clusters, the smaller one consisting of immediate early genes, whose effect has been described extensively^20,33. The other cluster consisted of a set of major regulators, with downstream effects mainly on calcium-, ion-, and potassium-channels that have also been found in the differential expression analysis. Interestingly, we were able to show CocUD-specific effects, not on TF activity per se, but on the activation and suppression of downstream gene expression. We observed that ZEB1 activates genes related to synaptic signaling in a CocUD-specific subcluster of MSNs, ZEB1 could therefore potentially mediate this well-described effect ²⁸. Interestingly, ZEB1 was also recently identified as a major transcriptional regulator in Alzheimer’s disease³⁴, providing another link between neurodegenerative diseases and SUDs. We further explored the druggability of the highly connected downstream targets, revealing phosphodiesterase inhibitors as a candidate drug, consistent with previous studies in rodent models showing reduced acquisition of cocaine place preference and cue-induced reinstatement of cocaine seeking ^35,36. Because of the large overlap of DE genes in the D1- and D2-MSNs, we constructed the GRN using information from both cell types. It has to be noted that while the transcriptional signatures of the two cell types are largely overlapping, a recent study has shown that transcriptional changes in D1-MSNs are more prevalent in the acute response to cocaine and in D2-MSNs more reflective of a cumulative effect of chronic exposure ³⁷.

We explored whether broad gene expression changes of GRN-identified targets could be observed in the 3-crit model of cocaine use disorder. Here, we validated the effects for PDE10A and CACNA1C in comparing 0-crit, but not 3-crit, to control animals. This highlights an important challenge in addiction research: based on postmortem human brain tissue alone, we cannot conclude whether the observed effects result from substance intake or reflect neuroadaptations as observed in substance use disorder.

In addition to the GRN, we performed another data-driven analysis, Multi-Omics Factor Analysis (MOFA). MOFA confirmed the results of the DE analysis, with the neuronal-specific factor pointing towards increased metabolic processes and synaptic signaling in CocUD. We also observed a microglia-specific factor associated with control status, providing further evidence that neuroinflammation is not elevated in CocUD¹⁸, as opposed to other SUDs, such as alcohol use disorder^38,39 and opioid use disorder¹⁸. Interestingly, we also observed a sex-specific factor associated with CocUD and enriched for ribosomal and metabolic processes. While sex-specific effects have been observed at large in rodent models, suggesting that female rats are more susceptible to cocaine consumption, evidence from human studies suggests an effect specific to the luteal phase of the menstrual cycle⁴⁰. In our exploratory analyses of potential sex-effects in the gene expression data, we observed strong sex-specific effects in oligodendrocytes, endothelial cells, and CALB2-GABAergic neurons, providing first evidence into potential drivers of sex-specific effects in CocUD.

At the same time, the sex-specific effects need to be interpreted with caution. While the single-nuclei framework allowed for a high statistical power, because each cell was treated as one observation, the overall sample size was still small and results could lack generalizability.

We did not observe cell type-specific enrichment of GWAS signals of CocUD. It has to be noted that the currently available GWAS samples of CocUD are small, comprising 2,482 cases and 836 controls of African American ancestry, and have not yet identified genome-wide significant SNPs⁴¹. Future studies increasing the sample size of CocUD GWAS are needed⁴, providing the basis to interpret GWAS signals in multi-omics studies.

Our results replicate known neurobiological mechanisms in CocUD and identify new targets, with complex CocUD interaction effects, particularly ZEB1, for future investigation. The present analysis highlights the potential of using a multi-omics framework allowed us to analyze the data in a more mechanistic approach compared to single-omics datasets.

Human Tissue

Postmortem human brain tissue

Postmortem human brain tissue samples of the caudate nucleus from six individuals with and eight individuals without CocUD of African American descent were obtained from the University of Texas Health Science Center at Houston (UTHealth Houston) Brain Collection (Sup Table S1). CocUD donors had age > 18 and a lifetime diagnosis of cocaine use disorder based on next-of-kin interviews and psychological autopsy⁴². Controls had no history of substance use or neuropsychiatric disorders. For all cases, exclusion criteria included a known diagnosis of severe neurodevelopmental disorders.

Nuclei isolation

For each sample, nuclei were isolated from 10 mg of frozen postmortem tissue, in order to minimize the content of free-floating debris. Brain samples were processed following the 10x Genomics Chromium Nuclei Isolation kit with RNase Inhibitor user guide. Briefly, the tissues were dissociated in lysis buffer, filtered and the cellular debris eliminated with the designated buffer. Then, alternate steps of centrifugation and wash resulted in clean nuclei suspensions. Nuclei were automatically counted on LUNA-FL™ Dual Fluorescence Cell Counter and particles smaller than 5 µm of diameter and bigger than 15 µm of diameter were excluded from the counting.

Library preparation and sequencing

Single-nuclei RNA-seq and ATAC-seq libraries were generated using the 10x Genomics Chromium Next GEM Single Cell Multiome ATAC + Gene Expression kit and user guide. In brief, single nuclei were loaded into chromium chips J with a capture target of 10,000 nuclei per sample. Final libraries were prepared as indicated in the 10x Genomics protocol and sequenced at NGS Competence Center Tübingen (Tübingen, Germany) using an Illumina NovaSeq6000 sequencer with a depth of 50,000 read pairs per nuclei for the RNA libraries and 25,000 read pairs per nuclei for the ATAC libraries.

Single nuclei RNA- and ATAC sequencing data analysis

Preprocessing and quality control

Raw data obtained from single-nuclei RNA- and ATAC-sequencing was preprocessed using the Cell Ranger ARC software (10x Genomics, v.2.0.2) using the human genome built GRCh38 as reference. Resulting filtered feature barcode matrices were imported in R (v.4.2.1) and converted to a Seurat object using Seurat v.4.3.0²³. In Seurat, quality control was performed and nuclei were excluded if they did not meet the following parameters: number of RNA features between 200 and 6,500, number of ATAC features between 1,000 and 25,000, percentage of mitochondrial genes < 10%, TSS enrichment score between 2 and 8, and nucleosome signal between 0.5 and 2.

Data integration and clustering

After quality control, individual datasets from the 14 individuals were merged to a single Seurat object using the merge function. ATAC peaks were called in the merged object using MACS2⁴³ to obtain a uniform set of peaks. After normalization, variable feature identification, and data scaling of the RNA assay, we applied dimensionality reduction using the RunPCA function and integrated the expression values across the 14 individuals with harmony²⁵ (v.0.1.1) using the first 20 principal components from PCA, followed by UMAP reduction. In a similar approach using Signac (v.1.9.0)²⁴, the ATAC assay was normalized using TFIDF, and linear dimensionality reduction was performed using SVD. Based on the principal components 2-30, UMAP reduction was performed followed by data integration using harmony. Finally, the RNA and ATAC assays were combined using a weighted nearest neighbor methodology as implemented in the FindMultiModalNeighbors function in Seurat. Using a resolution threshold of 0.25, a total of 15 clusters were identified in the single-nuclei multiomic dataset.

Identification of cell types

The expression of known marker genes for cell types in the human striatum⁴⁴ and different types of expected interneurons was evaluated in the RNA assay. In addition, we utilized the FindAllMarkers function to generate differential expression profiles distinguishing expression in one cluster from all other clusters (Sup Fig. S5). We cross-checked highly expressed genes in clusters with the Transcriptomics Explorer of the Allen Human Brain Atlas⁴⁵ and the Human Protein Atlas⁴⁶ leading to the identification of 15 different cell type clusters. Two minor clusters were characterized by strong expression of mitochondrial genes as evaluated by the FindAllMarkers function in R. We excluded these clusters expressing oligodendrocytes and excitatory neuron markers from the object and performed further analysis in the 13-cluster object.

Differential expression analysis

Differentially expressed (DE) genes were determined based on a Wilcoxon rank sum test as implemented in the FindMakers function in Seurat. As differential expression cut-offs, absolute log2-fold change > 0.5, minimal expression in at least 25% of cells in either of the compared clusters and an adjusted p-value < 0.05 were used. DE analyses were restricted to clusters consisting of at least 100 nuclei in each condition (CocUD yes/no) to ensure sufficient statistical power in the comparison. Visualization of DE genes was performed using the R packages ComplexHeatmap⁴⁷ (v.2.14.0) and UpSetR⁴⁸ (v.1.4.0).

Differential chromatin accessibility analysis

Identification of differentially accessible (DA) peaks in the ATAC assay for each cluster was performed using FindMarkers based on peaks that were at least expressed in 5% of cells in each of the compared clusters. Differential accessibility cut-offs were absolute log2 fold-change > 0.25 and adjusted p-value < 0.05. DA analyses also were restricted to clusters with at least 100 nuclei per condition.

Genomic annotation and overrepresentation analysis of ATAC peaks

Annotation of DA peaks to genomic features based on chromosomal coordinates of the GRCh38 reference genome was performed using the annotatePeak function of the R package ChIPseeker⁴⁹ (v.1.34.1). We further tested the enrichment of DA peaks in different genomic feature categories using the annotation of all peaks in the ATAC assay as the reference. A chi-squared test was used to test for significant enrichment within a category and Bonferroni-correction was used to adjust resulting p-values for multiple testing.

Gene ontology enrichment analysis

To test the overrepresentation of upregulated and downregulated DE genes in GO terms for the different cell type clusters, we used the compareCluster function of the R package clusterProfiler⁵⁰ (v.4.6.2) with the list of significant DE genes as input. Results for GO terms with a significant enrichment (adjusted p-value < 0.05) were visualized as pie plots in GO term networks using enrichplot (v.1.18.3). The same approach was applied on gene lists resulting from the annotation of CocUD-associated differentially accessible distal and promoter peaks to genes with the intention of identifying GO terms related to differential accessibility.

Motif enrichment analysis

To assess the enrichment of transcriptional motifs within DA peaks associated with CocUD status, we performed motif enrichment analysis (MEA) using chromVAR⁵¹ (v.1.20.2). First, the RunChromVar function was applied to the full ATAC assay to identify motif activities (z-scores) in the dataset. Up- and downregulated differentially active motifs that are associated with CocUD status were assessed for each cell type cluster using the FindMarkers function resulting in an estimate for the average difference in ChromVAR z-scores. Differentially active motifs with an adjusted p-value < 0.05 were considered statistically significant.

Identification of gene regulatory networks

We used pando²⁷ to infer gene regulatory networks (GRN) based on snRNA and snATAC data for each cluster using the respective DE genes as the candidate gene list for inferring GRNs. Pando uses multi-omic data together with inferred transcription factor binding sites to generate a GRN consisting of transcription factors and their inferred targets. After initiating the GRN, transcription factor modules were identified using the following parameters: p-value threshold = 0.05, number of variables = 5, minimum number of genes per module = 5, and R²=0.1. Network graphs for TFs and targets in each of the evaluated clusters were generated using the force-directed Fruchterman-Reingold layout.

Druggability analysis of GRN targets

Reference data from the drug-gene interaction database⁵² (DGIdb) was used to evaluate drug-gene interactions and potential druggability of GRN targets. Gene targets from the GRN were prioritized based on the number of statistically significant TF-target associations inferred by Pando. We determined quantiles of the TF-target connectivity distribution and selected targets with a minimum number of N=17 significant TF-target connections corresponding to the value related to the 95% quantile. All of the identified target genes were connected to multiple TFs (mean=16 upstream TFs) in the GRN network highlighting their potential use as convergent downstream drug targets. The prioritized set of target genes was queried using the DGIdb online tool (https://www.dgidb.org, v.5.0.5) with “Interactions by Gene” mode. Results for unique matches of gene-drug interactions were downloaded and sorted by gene-drug interaction score and FDA approval state of connected drugs.

Multi omics factor analysis

We applied multi omics factor analysis (MOFA)⁵³ to identify factors that are correlated with CocUD consisting of features from both the ATAC and the RNA assay. We separated the ATAC assay by peak annotation to “distal” and “promoter” category as performed in Cell Ranger Arc to investigate the relationship between ATAC peaks in different genomic regions with RNA expression levels in the factor analysis. After the identification of highest variable features in the RNA (N=5,000), ATAC promoter (N=16,700) and ATAC distal (N=42,651) assays, a MOFA object based on 15 factors was trained in R using MOFA2 (v.1.8.0). In the resulting object, we tested for association of the resulting factors with CocUD status and performed MOFA-implemented downstream analyses such as gene set enrichment analysis (GSEA) and motif enrichment analyses in the CocUD-associated factors.

Rank-Rank Hypergeometric Overlap

We utilized the R package RRHO2⁵⁴ (v.1.074) to conduct rank-rank hypergeometric overlap (RRHO) analysis, aiming to assess convergent and divergent expression patterns at the transcriptome-wide level between male and female donors. RRHO scores were derived from comprehensive differential expression statistics from both male and female datasets. Subsequently, we evaluated overlapping signatures using the hypergeometric testing procedure implemented in RRHO2.

GWAS enrichment analysis

To evaluate the expression signatures of genes containing risk variants associated with substance use disorders including CocUD we performed a GWAS enrichment analysis using single-cell disease-relevance scoring (scDRS). Here, we aimed to identify cell types of the caudate nucleus that are enriched for the expression of risk genes and might therefore be especially susceptible to the effects of genetic risk variants. The following GWAS summary statistics for African American (AA) populations were included in the analysis: cocaine dependence (N=3,318, Gelernter et al., 2014)⁴¹, alcohol dependence⁵⁵ (N=6,280, Walters et al., 2018), cannabis use disorder⁵⁶ (N=9,745, Johnson et al., 2020), and opioid use disorder⁵⁷ (N=2555, Polimanti et al., 2020). First, a gene-based analysis was performed for each summary statistics file using MAGMA⁵⁸ v.1.10. Next, scDRS was performed according with default settings as outlined in https://martinjzhang.github.io/scDRS/notebooks/quickstart.html. First, a gene-p-value matrix was created including information on the strength of the association between the gene and the phenotype. The scdrs munge-gs function was applied to convert the GWAS-based gene statistics to scDRS format. Next, single-cell expression data was extracted from the Seurat object and converted to h5ad file-format. Finally, using the h5ad-file and the munged scDRS input file, disease-relevance scoring was performed using the compute-score function in scDRS. We used the scdrs perform-downstream function to evaluate group-level (i.e. cell cluster) associations with the inferred disease score.

Validation

3-crit model of cocaine use disorder and rodent tissue

RT-qPCR was performed in tissue from the rat dorsal striatum, corresponding to the human caudate putamen, in male Sprague-Dawley rats that underwent the multisymptomatic 3-CRIT animal model of CocUD. Briefly, following extended intravenous self-administration (45 days), rats were tested on 3 measures modeled from criteria used to diagnose human use disorder in DSM-IV/V: the persistence of drug-seeking (responding when cocaine is signaled as unavailable), motivation for drug-seeking (progressive ratio), and drug taking despite adverse consequences (responding for cocaine in the presence of foot shock). The 3-CRIT model results in a distribution of animals ranging from those showing no addicted-like behavior (0-crit) to those showing addicted-like behavior on all three criteria (3-crit) (further described in Pohořalá, et al. ²²). For RT-qPCR in the current experiment, we analyzed tissue from six 0-crit and five 3-crit rats, as well as six control rats (N_total = 17) that were confined to the home cage for the entirety of experiment and thus had no exposure to cocaine or behavioral procedures. All rats were sacrificed 2d following the completion of the 3-CRIT protocol, at which point brains were flash-frozen in isopentane and stored at -80° C until further processing. Experimental procedures were conducted according to the NIH ethical guidelines for the care and use of laboratory animals, in compliance with the German Animal Welfare Act, and approved by the local animal care committee (Regierungspräsidium Karlsruhe, Germany).

Quantitative Real-Time PCR in the 3-crit model

Total RNA was isolated from rat brain tissues using the RNeasy Mini Kit (Qiagen) and reverse transcribed into cDNA using the High-Capacity cDNA Reverse Transcription Kit (4368813; ThermoFisher). Quantitative real-time PCR (qPCR) was performed on the QuantStudio™ 7 Flex Real-Time PCR System (Applied Biosystems) with the following reaction setup: 10 µL TaqMan™ Gene Expression Master Mix (Applied Biosystems, catalog number: 4369016), 1 µL TaqMan™ Gene Expression Assay (Thermo Fisher Scientific), and 9 µL cDNA. Specific TaqMan™ Gene Expression Assays were used to detect the genes of interest: PDE10A (Rn00673152_m1), KCNIP4 (Rn00710466_m1), CACNA1C (Rn00709287_m1), and ZEB1 (Rn01538408_m1). The housekeeping gene GUSB (Rn00566655_m1) served as internal control. Gene expression analysis was performed in triplicates for the genes of interest and in duplicates for the internal control, with all samples analyzed together on the same plate per gene.

Statistical Analysis of 3-crit model qPCR data

We first generated the delta CT (ΔCt) values by subtracting the Ct mean values of the triplicates of the genes of interest (ZEB1, CACNA1C, PDE10A, KCNIP4) from the Ct mean values of the duplicates of the internal control GUSB. Then, we used Student's t-tests to compare the ΔCt values between different groups.

Acknowledgment

We are grateful for the invaluable donations and participation from families, as well as for the generous collaboration of the medical examiners at the Harris County Institute of Forensic Sciences. We also thank all colleagues who have contributed to the UTHealth Brain Collection and whose work continues to grow and broaden the collection.

We thank Joel Gelernter and his group for sharing the summary statistics of the cocaine dependence GWAS⁴¹.

Figures 1A, 4E and the graphical abstract were created with BioRender.com

P.K., L.Z., A.A. and A.J. acknowledge the generous support of the Hector Stiftung II.

Ethics statement

Postmortem tissue was sampled at the University of Texas Health Science Center at Houston based on written informed consent from the next-of-kin of each donor. The study was conducted in line with the guidelines of the Declaration of Helsinki and was approved by the Ethics Committee II of the Medical Faculty Mannheim under the register number 2021-681 and by the University of Texas Health Science Center at Houston (HSC-MS-15-0247).

Data and code availability statement

Raw data is currently being uploaded to the European Genome Archive (EGA). The code that has been used to analyze and visualize the dataset is deposited in the Github repository https://github.com/lzillich/CN_multiome_cocaine.

Funding

Financial support for this work was provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 402170461 (B10) – TRR 265⁵⁹ to S.H.W., A.M.M.O., and R.E.B., and Project-ID 43081111 to P.K.. The study was funded by the German Federal Ministry of Education and Research (BMBF) through the e:Med consortium “A systems-medicine approach towards distinct and shared resilience and pathological mechanisms of substance use disorders - SysMedSUDs” (01ZX01909 to R.S., P.K., M.R., A.C.H., S.H.W.). FS is supported by a Families for Borderline Personality Disorder Research (Beth and Rob Elliott) 2023 NARSAD Young Investigator Grant (# 31537). R01DA044859 and the John S. Dunn Foundation to C.W.B..

Conflicts of Interest

All authors report no conflicts of interest.

1 Substance Abuse and Mental Health Services Administration. Key substance use and mental health indicators in the United States: Results from the 2021 National Survey on Drug Use and Health (Center for Behavioral Health Statistics and Quality, Substance Abuse and Mental Health Services Administration, 2022).

2 Riezzo, I. et al. Side effects of cocaine abuse: multiorgan toxicity and pathological consequences. Curr Med Chem 19, 5624-5646 (2012). https://doi.org:10.2174/092986712803988893

3 Shiels, M. S., Freedman, N. D., Thomas, D. & Berrington de Gonzalez, A. Trends in U.S. Drug Overdose Deaths in Non-Hispanic Black, Hispanic, and Non-Hispanic White Persons, 2000-2015. Ann Intern Med 168, 453-455 (2018). https://doi.org:10.7326/m17-1812

4 Fernàndez-Castillo, N., Cabana-Domínguez, J., Corominas, R. & Cormand, B. Molecular genetics of cocaine use disorders in humans. Mol Psychiatry 27, 624-639 (2022). https://doi.org:10.1038/s41380-021-01256-1

5 Gibney, E. R. & Nolan, C. M. Epigenetics and gene expression. Heredity (Edinb) 105, 4-13 (2010). https://doi.org:10.1038/hdy.2010.54

6 Nestler, E. J. Epigenetic mechanisms of drug addiction. Neuropharmacology 76 Pt B, 259-268 (2014). https://doi.org:10.1016/j.neuropharm.2013.04.004

7 Volkow, Nora D. & Morales, M. The Brain on Drugs: From Reward to Addiction. Cell 162, 712-725 (2015). https://doi.org:https://doi.org/10.1016/j.cell.2015.07.046

8 Malvaez, M. & Wassum, K. M. Regulation of habit formation in the dorsal striatum. Curr Opin Behav Sci 20, 67-74 (2018). https://doi.org:10.1016/j.cobeha.2017.11.005

9 Belin, D. & Everitt, B. J. Cocaine seeking habits depend upon dopamine-dependent serial connectivity linking the ventral with the dorsal striatum. Neuron 57, 432-441 (2008). https://doi.org:10.1016/j.neuron.2007.12.019

10 Albertson, D. N. et al. Gene expression profile of the nucleus accumbens of human cocaine abusers: evidence for dysregulation of myelin. J Neurochem 88, 1211-1219 (2004). https://doi.org:10.1046/j.1471-4159.2003.02247.x

11 Zhou, Z., Yuan, Q., Mash, D. C. & Goldman, D. Substance-specific and shared transcription and epigenetic changes in the human hippocampus chronically exposed to cocaine and alcohol. Proc Natl Acad Sci U S A 108, 6626-6631 (2011). https://doi.org:10.1073/pnas.1018514108

12 Mash, D. C. et al. Gene expression in human hippocampus from cocaine abusers identifies genes which regulate extracellular matrix remodeling. PLoS One 2, e1187 (2007). https://doi.org:10.1371/journal.pone.0001187

13 Walker, D. M. et al. Cocaine Self-administration Alters Transcriptome-wide Responses in the Brain's Reward Circuitry. Biol Psychiatry 84, 867-880 (2018). https://doi.org:10.1016/j.biopsych.2018.04.009

14 McClung, C. A. & Nestler, E. J. Regulation of gene expression and cocaine reward by CREB and DeltaFosB. Nat Neurosci 6, 1208-1215 (2003). https://doi.org:10.1038/nn1143

15 Gao, P., Limpens, J. H., Spijker, S., Vanderschuren, L. J. & Voorn, P. Stable immediate early gene expression patterns in medial prefrontal cortex and striatum after long-term cocaine self-administration. Addict Biol 22, 354-368 (2017). https://doi.org:10.1111/adb.12330

16 Ribeiro, E. A. et al. Transcriptional and physiological adaptations in nucleus accumbens somatostatin interneurons that regulate behavioral responses to cocaine. Nat Commun 9, 3149 (2018). https://doi.org:10.1038/s41467-018-05657-9

17 Freeman, W. M. et al. Gene expression changes in the medial prefrontal cortex and nucleus accumbens following abstinence from cocaine self-administration. BMC Neurosci 11, 29 (2010). https://doi.org:10.1186/1471-2202-11-29

18 Mews, P. et al. Convergent abnormalities in striatal gene networks in human cocaine use disorder and mouse cocaine administration models. Sci Adv 9, eadd8946 (2023). https://doi.org:10.1126/sciadv.add8946

19 Phillips, R. A., 3rd et al. Distinct subpopulations of D1 medium spiny neurons exhibit unique transcriptional responsiveness to cocaine. Mol Cell Neurosci 125, 103849 (2023). https://doi.org:10.1016/j.mcn.2023.103849

20 Savell, K. E. et al. A dopamine-induced gene expression signature regulates neuronal function and cocaine response. Sci Adv 6, eaba4221 (2020). https://doi.org:10.1126/sciadv.aba4221

21 Deroche-Gamonet, V., Belin, D. & Piazza, P. V. Evidence for addiction-like behavior in the rat. Science 305, 1014-1017 (2004). https://doi.org:10.1126/science.1099020

22 Pohořalá, V., Enkel, T., Bartsch, D., Spanagel, R. & Bernardi, R. E. Sign- and goal-tracking score does not correlate with addiction-like behavior following prolonged cocaine self-administration. Psychopharmacology (Berl) 238, 2335-2346 (2021). https://doi.org:10.1007/s00213-021-05858-z

23 Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573-3587.e3529 (2021). https://doi.org:10.1016/j.cell.2021.04.048

24 Stuart, T., Srivastava, A., Madad, S., Lareau, C. A. & Satija, R. Single-cell chromatin state analysis with Signac. Nat Methods 18, 1333-1341 (2021). https://doi.org:10.1038/s41592-021-01282-5

25 Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 16, 1289-1296 (2019). https://doi.org:10.1038/s41592-019-0619-0

26 Driscoll, M. E., Bollu, P. C. & Tadi, P. Neuroanatomy, nucleus caudate. (2020).

27 Fleck, J. S. et al. Inferring and perturbing cell fate regulomes in human brain organoids. Nature 621, 365-372 (2023). https://doi.org:10.1038/s41586-022-05279-8

28 Wolf, M. E. Synaptic mechanisms underlying persistent cocaine craving. Nat Rev Neurosci 17, 351-365 (2016). https://doi.org:10.1038/nrn.2016.39

29 Mameli, M. et al. Cocaine-evoked synaptic plasticity: persistence in the VTA triggers adaptations in the NAc. Nat Neurosci 12, 1036-1041 (2009). https://doi.org:10.1038/nn.2367

30 Gomes, C. et al. RNA polymerase 1-driven transcription as a mediator of BDNF-induced neurite outgrowth. J Biol Chem 286, 4357-4363 (2011). https://doi.org:10.1074/jbc.M110.170134

31 Leighton, L. J. et al. Experience-dependent neural plasticity, learning, and memory in the era of epitranscriptomics. Genes Brain Behav 17, e12426 (2018). https://doi.org:10.1111/gbb.12426

32 Wang, J. et al. Cocaine Triggers Astrocyte-Mediated Synaptogenesis. Biol Psychiatry 89, 386-397 (2021). https://doi.org:10.1016/j.biopsych.2020.08.012

33 Teague, C. D. & Nestler, E. J. Key transcription factors mediating cocaine-induced plasticity in the nucleus accumbens. Mol Psychiatry 27, 687-709 (2022). https://doi.org:10.1038/s41380-021-01163-5

34 Anderson, A. G. et al. Single nucleus multiomics identifies ZEB1 and MAFB as candidate regulators of Alzheimer's disease-specific cis-regulatory elements. Cell Genom 3 (2023). https://doi.org:10.1016/j.xgen.2023.100263

35 Liddie, S., Anderson, K. L., Paz, A. & Itzhak, Y. The effect of phosphodiesterase inhibitors on the extinction of cocaine-induced conditioned place preference in mice. J Psychopharmacol 26, 1375-1382 (2012). https://doi.org:10.1177/0269881112447991

36 Zhong, P. et al. Phosphodiesterase 4 inhibition impairs cocaine-induced inhibitory synaptic plasticity and conditioned place preference. Neuropsychopharmacology 37, 2377-2387 (2012). https://doi.org:10.1038/npp.2012.93

37 Browne, C. J. et al. Shared and divergent transcriptomic regulation in nucleus accumbens D1 and D2 medium spiny neurons by cocaine and morphine. bioRxiv (2023). https://doi.org:10.1101/2023.09.19.558477

38 Zillich, L. et al. Epigenome-wide association study of alcohol use disorder in five brain regions. Neuropsychopharmacology 47, 832-839 (2022). https://doi.org:10.1038/s41386-021-01228-7

39 Zillich, L. et al. Multi-omics signatures of alcohol use disorder in the dorsal and ventral striatum. Transl Psychiatry 12, 190 (2022). https://doi.org:10.1038/s41398-022-01959-1

40 Evans, S. M. & Foltin, R. W. Does the response to cocaine differ as a function of sex or hormonal status in human and non-human primates? Horm Behav 58, 13-21 (2010). https://doi.org:10.1016/j.yhbeh.2009.08.010

41 Gelernter, J. et al. Genome-wide association study of cocaine dependence and related traits: FAM53B identified as a risk gene. Mol Psychiatry 19, 717-723 (2014). https://doi.org:10.1038/mp.2013.99

42 Meyer, T. D., Godfrey, C. J. & Walss-Bass, C. The UT health Psychological Autopsy Interview Schedule (UTH- PAIS) - Description and reliability of diagnoses and transdiagnostic personality measures. J Psychiatr Res 156, 221-227 (2022). https://doi.org:10.1016/j.jpsychires.2022.10.014

43 Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol 9, R137 (2008). https://doi.org:10.1186/gb-2008-9-9-r137

44 Brenner, E. et al. Single cell transcriptome profiling of the human alcohol-dependent brain. Hum Mol Genet 29, 1144-1153 (2020). https://doi.org:10.1093/hmg/ddaa038

45 Shen, E. H., Overly, C. C. & Jones, A. R. The Allen Human Brain Atlas: comprehensive gene expression mapping of the human brain. Trends Neurosci 35, 711-714 (2012). https://doi.org:10.1016/j.tins.2012.09.005

46 Uhlén, M. et al. Proteomics. Tissue-based map of the human proteome. Science 347, 1260419 (2015). https://doi.org:10.1126/science.1260419

47 Gu, Z., Eils, R. & Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847-2849 (2016). https://doi.org:10.1093/bioinformatics/btw313

48 Conway, J. R., Lex, A. & Gehlenborg, N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics 33, 2938-2940 (2017). https://doi.org:10.1093/bioinformatics/btx364

49 Wang, Q. et al. Exploring Epigenomic Datasets by ChIPseeker. Curr Protoc 2, e585 (2022). https://doi.org:10.1002/cpz1.585

50 Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141 (2021). https://doi.org:10.1016/j.xinn.2021.100141

51 Schep, A. N., Wu, B., Buenrostro, J. D. & Greenleaf, W. J. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nat Methods 14, 975-978 (2017). https://doi.org:10.1038/nmeth.4401

52 Cannon, M. et al. DGIdb 5.0: rebuilding the drug-gene interaction database for precision medicine and drug discovery platforms. Nucleic Acids Res 52, D1227-d1235 (2024). https://doi.org:10.1093/nar/gkad1040

53 Argelaguet, R. et al. Multi-Omics Factor Analysis-a framework for unsupervised integration of multi-omics data sets. Mol Syst Biol 14 (2018). https://doi.org:10.15252/msb.20178124

54 Cahill, K. M., Huo, Z., Tseng, G. C., Logan, R. W. & Seney, M. L. Improved identification of concordant and discordant gene expression signatures using an updated rank-rank hypergeometric overlap approach. Sci Rep 8, 9588 (2018). https://doi.org:10.1038/s41598-018-27903-2

55 Walters, R. K. et al. Transancestral GWAS of alcohol dependence reveals common genetic underpinnings with psychiatric disorders. Nat Neurosci 21, 1656-1669 (2018). https://doi.org:10.1038/s41593-018-0275-1

56 Johnson, E. C. et al. A large-scale genome-wide association study meta-analysis of cannabis use disorder. Lancet Psychiatry 7, 1032-1045 (2020). https://doi.org:10.1016/s2215-0366(20)30339-4

57 Polimanti, R. et al. Leveraging genome-wide data to investigate differences between opioid use vs. opioid dependence in 41,176 individuals from the Psychiatric Genomics Consortium. Mol Psychiatry 25, 1673-1687 (2020). https://doi.org:10.1038/s41380-020-0677-9

58 de Leeuw, C. A., Mooij, J. M., Heskes, T. & Posthuma, D. MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput Biol 11, e1004219 (2015). https://doi.org:10.1371/journal.pcbi.1004219

59 Spanagel, R. et al. The ReCoDe addiction research consortium: Losing and regaining control over drug intake-Findings and future perspectives. Addict Biol 29, e13419 (2024). https://doi.org:10.1111/adb.13419

There is NO Competing Interest.

FigureS1final.pdf
Figure S1. UMAPs colored by A) CocUD case-control status, B) donor, C) sex, D) PMI (hours), E) overdose. CocUD = cocaine use disorder, Ctrl = control, PMI = postmortem interval, 0 = no, 1 = yes.
FigureS2final.pdf
Figure S2. emapplot depicting the results of the Gene Ontology overrepresentation analysis in a cell type specific manner, for A) promoter peaks and B) distal peaks, circle size indicates the number of genes in the pathway.
FigureS3final.pdf
Figure S3. A) violin plots depicting the number of target genes per transcription factor for D1-/D2-MSN GRN. B) regulon of transcription factors from D1-/D2-MSN GRN weighted by log2foldchange from the differential expression analysis. C) regulon size of highly connected transcription factors from D1-/D2-MSN GRN. D) goodness of fit indices of D1-/D2-MSN GRN. E) coverage plot of PDE10A indicating TF binding sites. F) coverage plot of CACNA1C indicating TF binding sites, red = ctrl, blue = CocUD.
TF = transcription factor, GRN = gene regulatory network, MSN = medium spiny neurons.
FigureS4final.pdf
Figure S4. A) association of MOFA factors with CocUD, number of RNA features, number of ATAC features, postmortem interval, sex and celltypes (upper) and variance explained by RNA, distal ATAC peaks and promoter ATAC peaks for each factor. B) factor loadings for CocUD by cell type for CocUD-associated factors 6, 7, and 8. C) MOFA UMAP projection colored by cell types (left), factor6 loadings in MOFA UMAP (middle left), factor7 loadings in MOFA UMAP (middle right), factor8 loadings in MOFA UMAP (left). D) results from gene-set enrichment analysis of factor6 (left), factor7 (middle), and factor8 (right). E) enriched promoter motifs of factor8, positive loadings (left), negative (right). F) projection of scDRS results of cocaine dependence (Gelernter et al., 2014), alcohol dependence (Zhou et al., 2020), cannabis use disorder (Johnson et al., 2021), opioid use disorder (Polimanti et al., 2020).
MOFA = Multi-Omics Factor Analysis, CocUD = cocaine use disorder, nFeatureRNA = number of RNA features, nFeature ATAC = number of ATAC features, PMI = postmortem interval.
FigureS5final.pdf
Figure S5. Heatmap depicting the highly expressed genes in each cell type cluster compared to gene expression in all other cell type clusters. Yellow = high expression, pink = low = expression.
SuppTablesCNmultiome.xlsx

Download PDF

Version 1

posted

You are reading this latest preprint version

Cell type-specific Multi-Omics Analysis of Cocaine Use Disorder in the Human Caudate Nucleus

Status:

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1