Central amygdala single-nucleus atlas reveals chromatin and gene transcription dynamics in human alcohol use disorder

Domain: neuroscience
Published: 2026-01-19
Editorial outcome: minor_revision

Abstract
Regulation of gene expression is a highly coordinated process in both the healthy and pathological brain with unique patterns across a multitude of cell types. Here we present a multi-omic single nucleus study of ~175,000 nuclei from 50 donors with alcohol use disorder (AUD) and control donors without AUD, profiling cell type specific gene expression and chromatin accessibility in the human central amygdala. We identify all major CNS cell types and neuronal subtypes and find inhibitory neurons are particularly affected by AUD. We find high numbers of differentially expressed genes (DEGs) including GABRA2, GRM8, and NCAM1 and show significant enrichment for AUD risk genes within these DEGs. We identified 51,431 cell type-specific, disease associated candidate cis-regulatory elements including an interneuron-associated set of chromatin loops at the AUD risk gene CALN1. Transcription factor footprinting identified Kruppel-like factors upstream of AUD GWAS genes and DEGs. Finally, we also perform cell type-specific fine mapping for AUD GWAS to prioritize variants within functional genomic elements.

Introduction
Alcohol use disorder is a common neuropsychiatric disorder that is accompanied by significant social and vocational dysfunction, numerous medical co-morbidities and increased mortality. AUD, including alcohol abuse and dependence, is moderately heritable and is caused by genetic and environmental factors1,2. AUD has a prevalence of about 6.2% in the adult American population and genetic studies suggest heritability of 50–60%3–6.

Genome-wide association studies (GWAS) have demonstrated that genetic variation also contributes to individual vulnerability to AUD7. There are many published AUD GWAS that have identified genome wide-significant risk loci1,5,8–17. However, the link between genetic risk and molecular brain pathology is still far from clear. The pursuit of the molecular underpinnings of AUD risk is greatly complicated by the direct and indirect consequences of acute and chronic alcohol exposure and other AUD-related effects (e.g. malnutrition)18.

Current clinical and basic research suggests important, complex, and currently inadequately understood interactions between alcohol use and AUD neuropathology19,20. Alcohol crosses the blood brain barrier and triggers changes in the CNS, including reduced glucose uptake and disturbances in neurotransmitter release including glutamatergic, GABAergic and dopaminergic alterations. Although the mechanistic basis for these interactions is currently not well understood, candidates from previous studies have included alcohol effects on neurons and astrocytes21,22, inflammatory cytokines and immune mediators23, brain inhibitory function24,25 and synaptic signaling26 specifically at GABA receptors, where ethanol can directly bind and modulate behaviors including anxiolysis, disinhibition, and hyperlocomotion. This suggests alcohol may alter the epigenome resulting in long lasting aberrant effects on the brain related to these processes13.

Single-cell atlases have emerged as a promising tool for determining both the gene regulatory mechanisms and functions of disease related genes in specific brain regions. Consistent with preclinical studies, neuroimaging studies in AUD patients implicate multiple brain regions in impulsivity, alcohol reward, and craving including the medial prefrontal cortex (PFC)27, nucleus accumbens28, amygdala29 and hippocampus30. Recent studies have demonstrated alcohol-use disorders are associated with negative reinforcing states and that individuals may drink alcohol to reduce anxiety or withdrawal symptoms31. A subset of individuals self-medicate with alcohol to reduce anxiety and stress and these stress-dampening properties are mediated through the amygdala32. The central nucleus of the amygdala (CeA) is the primary output nucleus of the amygdala and is largely made up of inhibitory neurons. It has also been implicated in the reinforcing effects of alcohol and other drugs of abuse33. To date there have been no single-cell multi-omic studies of the human CeA in AUD. Therefore, we performed single nucleus (sn)ATAC-seq and snRNA-seq from the same nuclei (snMultiome) to define AUD-associated gene regulatory programs at the epigenomic and transcriptomic levels. This work provides a powerful look into the cellular heterogeneity of the CeA and allows us to unravel critical biological systems underlying alcohol dependence in the amygdala.

Results
Single cell genomic atlas of AUD central amygdala
We used single-nucleus multiome-sequencing (snMultiome-seq) to profile nuclei isolated from the central amygdala of post mortem brain of 50 donors from University of Pittsburgh Medical Center brain bank. The sample included 22 donors with a history of alcohol abuse and dependence (alcohol use disorder) and 28 CONs with no history of alcohol abuse. Complete demographics for the entire cohort are included in Supplementary Data 1. To dissect the single cell-type gene expression and gene regulatory changes in the AUD brain, we used a custom bioinformatic pipeline for analysis of each genomic layer and integration34. After rigorous quality-control and filtering, we report results from a discovery cohort of 174,188 droplet based nuclei (2081 genes and 6028 transcripts per nuclei) from snMultiome (10X Chromium Single Cell Multiome ATAC + Gene Expression) (Supplementary Figs. 1, 2)35. The snMultiome technology uniquely tags DNA and RNA from the same cell to measure chromatin accessibility and transcript levels from the same nuclei.

After UMAP clustering, marker genes analysis and comparison with previous studies of the whole amygdala, we annotated all major cell types in the CNS including inhibitory neurons (INH; 23,876 nuclei, 13.71%), excitatory neurons (EXC; 8643 nuclei, 4.96%), oligodendrocytes (OLI; 80,549 nuclei, 46.24%), oligodendrocyte progenitor cells (OPC; 15,341 nuclei, 8.81%), endothelial cells (END; 3316 nuclei, 1.90%), astrocytes (AST; 27,313 nuclei, 15.68%) and microglia (MIC; 15,150 nuclei, 8.70%). Notably, INH neurons represent ~75% of the central amygdala neurons (with a median of ~80% INH neuron composition across all samples) consistent with previous neuroanatomical counts36. To visualize our data across genomic modalities, we generated UMAP clusters for each using the same nuclei. RNA and ATAC UMAPs were created based on gene expression levels (Fig. 1a) and chromatin accessibility (Fig. 1b), respectively. Gene expression and chromatin accessibility from the same nuclei were integrated to generate a snMutliome UMAP (Fig. 1c). We calculated cell type proportions based on the number of nuclei isolated (Fig. 1d) and found no significant proportion shifts due to disease (AUD versus CON): INH P-value=0.718; EXC P-value=0.324; OLI P-value=0.977; OPC P-value=0.406; END P-value=0.551; AST P-value=0.945; and MIC P-value=0.945. We plotted average gene expression levels for each cell type marker transcripts (Fig. 1e), and normalized chromatin signals for open ATAC peaks for each marker gene (Fig. 1f). Transcriptomic subtypes were largely conserved across individuals (Supplementary Fig. 2a, b).

With the aid of previously observed cell type markers of the amygdala36–38 and using de novo marker detection, we were able to identify subtypes for both our neuronal and non-neuronal cells. We identified six INH neuron subtypes including DRD1, LHX9, SST, VIP, PENK, and FREM1 positive, two oligodendroctye subtypes, differentiated endothelial and mural cells, and separated microglia and macrophages, with numbers of cells significantly higher than what has been previously reported38 (Fig. 2a). Average gene expression levels for individual marker transcripts were used to classify each subtype (Fig. 2b). We found the highest percentage of INH neurons were the DRD1 (31.6%) and LHX9 (21.7%) positive neurons. The percentage of the remaining INH subtypes were SST, 12.8%; VIP, 9.7%; PENK, 14.2%; FREM1, 10.0% (Fig. 2c). We identified two OLI subtypes of unequal proportion proportions including oligo-1, 78.2%; oligo-2, 21.8%) (Fig. 2d). We identified two vascular subtypes (endothelial cells, 39.8% and mural cells 60.2%) (Fig. 2e); and two immune subtypes (microglia, 97.7% and macrophages, 2.3%) (Fig. 2f). To confirm our immune cell subtypes, we examined expression of additional microglia-specific markers P2RY12, TMEM119, and found high co-expression with APBB1IP. Additionally, we examined the expression of infiltrating macrophage-specific marker SIGLEC1 (aka CD169) and found high co-expression with our F13A1 macrophage cluster (Supplementary Fig. 3a). Cell subtypes were largely conserved across individuals, as most of the identified cellular subsets contained nuclei from the majority of donors (Supplementary Fig. 2a). To confirm selected interneuron markers, we used immunohistochemistry to colabel each marker with the GABAergic marker GAD1. We used acetylcholinesterase histochemistry to confirm the location of the CeA (lightest stain) compared with surrounding amygdala structures (Supplementary Fig. 3b). We confirmed co-expression of PENK (Supplementary Fig. 3c) and FREM (Supplementary Fig. 3d) with GAD1 in amygdala tissue sections adjacent to those used for acetylcholinesterase staining.

Systematic differential gene expression implicates inhibitory neuron vulnerability in AUD
To better determine how individual cell types are affected by alcohol abuse, we analyzed differential gene expression across all cell type clusters. To ensure the highest rigor in our differential gene expression analysis, we combined two extensively benchmarked methods for DEG detection: parametric-based MAST39 and permutation-based Wilcoxon Test40 (Methods). For each cell cluster, we report all DEGs that were identified as overlapping between the two pipelines (Supplementary Data 2). For our DEG analysis, we analyzed canonical cell types (Fig. 3a, left panel) but also examined gene expression in specific INH neuronal subtypes (Fig. 3a, right panel). DEGs significantly associated with AUD were identified in all cell types. We also found that DEGs (769) were enriched in INH neurons compared to the EXC neurons, with the majority (466, 60.6%) occurring in INH PENK neurons. We identified 4364 unique protein-coding AUD DEGs using MAST (Supplementary Data 3) and 3290 unique protein-coding DEGs from the Wilcoxon test (Supplementary Data 4). Of those, 2270 DEGs overlapped in the same direction, shared directional fold change (FC)>1.2 with FDR<0.01, and included 1805 unique genes (Supplementary Data 5). In general, we found that neuronal cell types shared the highest gene expression changes, with 5.93% DEGs shared between both EXC and INH classes (Fig. 3b). No gene was regulated across all 7 major cell types. Using cosine similarity, we compared AUD DEGs across cell types and found the highest concordance of gene expression among neuronal subtypes (Supplementary Figs. 4a, 3a). We also performed DEG analysis on all cell subtypes (Supplementary Fig. 5b, c). In general, we found that the subtypes shared the majority of DEGs with their major cell type clusters (~89% of upregulated and ~98% of downregulated DEGs) (Supplementary Fig. 5a) with 95% of DEGs identified in subtypes also found in the major cell types (Supplementary Fig. 5d).

To assess the impact of drug use and comorbidities, we performed additional analyses that included covariates for manner of death, lifetime alcohol consumption, illicit drug use, comorbid psychiatric disorders, and medications at time of death (Supplementary Data 6). Across cell types, ~80% of DEGs overlapped with our model (Supplementary Fig. 6a). Additionally, to ensure rigor, we tested a sample-based DEG method (i.e. pseudobulk) using DESeq241 (Supplementary Data 7). We identified additional DEGs: 14 in inhibitory neurons, 26 in excitatory neurons, 1 in oligodendrocyte, 10 in astrocytes, and 1 in microglia (Supplementary Fig. 6b). Finally, we also observed sex-specific responses (Supplementary Data 8–10), consistent with previous observations in AUD42,43. Specifically, males have the most DEGs in the INH cell type, whereas female DEGs are more uniformly distributed across cell types (Supplementary Fig. 6c, e). We also found pathway enrichment differences, specifically in WNT5A signaling in females and ERBB4 signaling in males (Supplementary Fig. 6d, f). In addition, male pathways tended to be more cell type-specific, while females’ tended to be enriched across cell types. 710 DEGs were common across the male and female specific comparisons or approximately 55% or 25% overlap, respectively.

Among top DEGs per major cell type, we found little overlap (Fig. 3c). Due to the high sensitivity of our analyses, we were able to identify several genes previously implicated in alcohol addiction including KCNIP1 in both INH neurons and MIC, as well as VAMP2 which has previously been shown to stimulate GABA neurotransmitter release after alcohol intake44. To better understand the global molecular processes affected by AUD, we performed gene set enrichment analysis on all 1805 AUD DEGs (Supplementary Fig. 4b and Supplementary Data 11). We observed significant enrichment for neuronal pathways related to axon guidance, glutamate receptors and cell to cell signaling. Individual cell type-specific gene set enrichment analysis (Fig. 3d, Supplementary Fig. 4c, Supplementary Data 12, 13) identified pathways related to GABA synthesis, glutamate receptor binding, and P2Y receptor signaling. Taken together, differential gene expression changes in AUD amygdala are mostly neuronal and affect systems related to neurotransmitter signaling and immune processes.

Interestingly, we observe a statistically significant enrichment of 19 PAU/AUD GWAS risk genes (P-val=5.32e−3, O.R.=2.23) (Fig. 3e and Supplementary Fig. 7a, b). We highlight the gene expression changes in the INH neurons of the amygdala (Fig. 3e) as these are the most abundant neuronal cell type and have the most number of DEGs. We found significant decreases in GRM8 (logFC=−0.485 and FDR=5.50e−5) and increases in GABRA2 (logFC=0.588 and FDR=6.73e−19 in DRD interneurons. We also observed an increase in GRM8 (logFC=1.26 and FDR=1.73e−4) and a decrease in NCAM1 (logFC=−0.353 and FDR=1.06e−12) in PENK neurons and we confirmed these expression changes using fluorescent in situ hybridization (Supplementary Fig. 7c-j). There is extensive evidence for a role of glutamate neurotransmission dysfunction in alcohol dependence45,46 and related addiction behaviors and several metabotropic glutamate receptors have been implicated in these processes47.

Cis-gene regulation of calcium activity and glutamatergic signaling genes disrupted in AUD
In addition to measuring gene expression differences in CeA of the AUD brain, we measured chromatin accessibility within the same nuclei. This allowed us to identify cell type-specific open chromatin peaks and cis-gene regulatory mechanisms. Since the majority of common risk variants for neuropsychiatric disorders are located within non-coding regions of the genome, they likely functionally disrupt gene expression through cis-regulating elements (CREs). To identify open chromatin peaks we utilized MACS248, and found 789,247 peaks across all cell types, the majority of which were located in introns (50.0%) or distal intergenic regions (24.9%) (Fig. 4a, Supplementary Data 14). We then used the merged peak set of fixed-width 501-bp union peaks (378,931 total) based on the iterative overlap merging procedure which preserves cell type-specific peak annotation (Methods). This union peakset consists of peaks across all seven canonical cell types and is used in all downstream analyses (Supplementary Data 15). We found that nearly half of all peaks are cell type-specific, highlighting the importance of studying individual cell type epigenomic signatures in the brain (Supplementary Fig. 8a). We found the highest overlap in peaks between EXC and INH neurons (36% Jaccard index overlap). We identified robust peaks for cell type markers for the 7 major cell types (Supplementary Fig. 8b). Additionally, we identified peaks in cell subtypes and found that a significant portion (55.1%) differed from their major cell types (Supplementary Fig. 8c, d) with most peaks in INH DRD+ and PENK+ neurons, suggesting a high degree of specificity in gene expression control among INH subtypes. Finally, we performed gene ontology enrichment on cell type marker peaks and found that they were enriched for expected cell type-specific biological processes (Supplementary Fig. 8e).

We then found 1,526,912 peak-to-gene links (correlation>0.45, FDR<1×10−4) using the union peaks and plotted side by side heatmaps of linked ATAC regions and gene expression (Fig. 4b). Peak-to-gene links are correlations between chromatin peak accessibility and nearby gene expression that can be calculated by integrating the modalities of snMultiome data. Since snMultiome simultaneously profiles the gene expression and chromatin accessibility of the same nuclei, our peak-to-gene links profiles are matched from the same nuclei. Previous studies have shown that not all peaks are CREs controlling gene expression in human brain tissue49. To determine which peaks are most likely to regulate gene expression, we calculated the correlation between peak accessibility and target gene expression (threshold of correlation > 0.4 and FDR < 0.05). We found that across all cell types approximately 92.64% of peaks showed strong gene expression regulatory activity (Fig. 4c).

Because we identified an extensive number of DEGs across all cell types and significant enrichment of DEGs near risk variants for AUD, we reasoned that cis-regulatory architecture may be altered in the AUD brain. To identify disease gene regulatory mechanisms, we intersected our CRE linked genes (CLGs) with DEGs for each cell type. 36.05% of CREs were linked to an upregulated DEG. We found the most overlap in CLGs in INH neurons corresponding to the high number of DEGs we observe in those cells (Fig. 4d). We also found that CLG target genes that overlap between cell types are typically linked to different CREs. For example, GABRA2 is upregulated in INH neurons and downregulated in AST (Fig. 3e and Supplementary Data 16) with each having unique, non-overlapping CREs. We also highlight the peak-to-gene links for AUD risk gene CALN1. CALN1 is a calcium binding gene that is highly and exclusively expressed in the neuronal cell types (Fig. 4e, f) and was significantly upregulated in INH (Fig. 3e). We find an extensive number of ATAC peak-to-gene links (strong CREs, correlation>0.85) to CALN1’s TSS which are all cell type-specific (Fig. 4g dashed box).

We also identified peak-to-gene links (strong CREs, correlation>0.8) for genes previously implicated in addiction in other cell types. For example in astrocytes, we find an extensive number of ATAC peak-to-gene links to NKAIN3 and FGFR3. NKAIN3 is a sodium/potassium pump gene that is highly expressed in astrocytes (Supplementary Fig. 9a, d) and changes in FGFR3 transcript have been linked to alcohol exposure of offspring of binge drinking behavior50 (Supplementary Fig. 9b, e); both TSS’s are cell type-specific to astrocytes (Supplementary Fig. 9c, f dashed box). Additionally, both genes were significantly down regulated in AST (Supplementary Fig. 9g). Taken together, these results demonstrate unique epigenomic processes regulating gene expression in CeA and point to mechanisms affecting excitatory and calcium signaling driven by AUD.

Cell type-specific TF gene regulation alters GWAS and DE genes in AUD
We mapped transcription factor binding motifs within our CREs that correspond to differentially expressed genes in cells from our AUD cohort to understand the upstream regulatory mechanisms governing gene expression. We used chromVAR51 to compute TF motif variability in individual cell clusters by estimating the enrichment of TF binding motifs in accessible chromatin regions by cell type (Supplementary Fig. 10a). To obtain motifs with high enrichment in the CeA, we prioritize peak-to-gene links connecting to AUD DEGs (Supplementary Fig. 10b). We highlight TF RFX2 which has high enrichment in the INH and EXC neurons and AST (Fig. 5a). We then performed TF footprinting analysis to investigate binding dynamics of RFX2 within the genome of these cell types (Fig. 5b) and found significant enrichment for RFX2 motifs in both EXC and INH neurons.

To calculate disease-specific differences in chromatin accessibility, we focused on peak-to-gene linkages specific for AUD. After calculating cell type-specific peak-to-gene links across all samples, we isolated peaks with strong differential peak-to-gene correlation (Pearson Correlation > 0.5) in the AUD cohort and considered these ‘AUD-specific peaks’ for downstream analysis (Fig. 5c and Methods). We found the majority of AUD-specific peaks in OPCs and END (24.13 and 24.09%, respectively) and INH (21.77%) (Supplementary Data 17–23). We built cell-type-specific gene networks (Supplementary Data 24–30) and INH motifs were chosen based on the number of DEGs and GWAS genes (Supplementary Data 31). We then performed enrichment analysis using these AUD-specific peaks and found an overrepresentation of TFs belonging to the KLF (Kruppel-like) family (Fig. 5d and Supplementary Data 32). We found significant enrichment for KLF16 motifs in INH, END, and MIC cells (Fig. 5e). Transcription factor footprinting of KLF16 revealed stronger enrichment in AUD INHs compared to CON INHs (Fig. 5f), reflecting the abundance of KLF motifs in the AUD peak-to-gene links motif enrichment described above. We constructed a gene regulatory network for AUD INHs (Fig. 5g) and found KLF16 was linked to 629 DEGs and 64 AUD risk genes including GRIN2A (upregulated), GRIK1 (upregulated), SMAD3 (downregulated), and OPRM1 (risk gene) (Fig. 5h). We found significant enrichment for DEGs in predicted downstream targets of the KLFs (P-val 1.26e−53; O.R.=4.052) and nominal significance for GWAS genes downstream (P-val 4.052; O.R.=0.643).

In order to functionally validate our genomic findings, we employed a mouse model of alcohol abuse. We isolated frozen amygdala from mice treated with an 2g/kg injection of 100% EtOH. We performed single nucleus isolation and snMultiome in the same manner as our human donor cohorts. We identified all major cell types across all modalities: RNA (Supplementary Fig. 11a), ATAC (Supplementary Fig. 11b) and combined Multiome (Supplementary Fig. 11c) for 19,072 individual nuclei. We identified 1814 total DEGs across all major cell types with most falling in the INH (Supplementary Fig. 11d). When compared with the human AUD DEG dataset we found overlap for 678 DEGs (37.6%) (Supplementary Fig. 11e). Most DEGs overlapped in neuronal cell types and this was reelected in gene set enrichment analysis which identified terms related to neuron projection, synaptic transmission, and glutamatergic signaling (Supplementary Fig. 11f). In addition, we intersected our INH DEGs with a previous central amygdala mouse model of alcohol withdrawal52 and found 38 overlapping DEGs in INH and EXC neurons and MIC (Supplementary Fig. 11g). We performed peak-to-gene linkage analysis on our mouse snMultiome data and identified 7485 total in the mouse AUD model. Of these, 4710 overlapped with the human AUD CRE genes (62.9%) (Supplementary Fig. 11h). To confirm the function of our KLF gene regulatory network in human AUD (Fig. 5g), we constructed a GRN for these same TFs in our mouse model (Supplementary Fig. 11i). We found that combined KLF16, KLF7 and KLF6 form a GRN that regulates 549 DEGs. We found 50 DEGs overlapped between the human and mouse GRNs including GWAS gene CALN1 confirming its likely role in alcohol abuse. Taken together, these data support the functional significance of KLF TFs in regulating calcium and glutamatergic signaling pathways identified in our dataset.

Cell-type specific variant fine mapping within the AUD gene regulatory landscape
In order to assess the role genetic risk variants play in AUD brain pathology, we performed linkage disequilibrium score regression (LDSC) for AUD and PAU (problematic alcohol use disorder includes individuals with AUD and alcohol related traits) across GWAS’ including European ancestry, trans-ancestry and sex comparisons. We also included other commonly comorbid psychiatric disorders with AUD, and additional control traits to better understand the landscape of genetic risk within the CeA (Fig. 6a, Supplementary Data 33). We found significant enrichment for AUD and PAU GWAS signals in both INH and EXC neurons across all datasets. Only one GWAS (Audit-C, alcohol use disorders identification test-consumption)14 showed significant enrichment in INH alone. Finally, we performed an LDSC enrichment by genic region and observed that AUD risk SNPs were enriched within introns of INH (Fig. 6b, Supplementary Data 34).

Because CREs regulate expression of nearby genes (functioning as promoters and enhancers), we reasoned that GWAS risk variants within them were more likely to be biologically impactful on the AUD brain transcriptome. We performed fine-mapping using SuSiE53 to identify risk variants with epigenetic information from our single cell ATAC peaks for CREs across all cell types. We integrated cell type specific CRE information with 131 PAU genome-wide risk genes from the largest alcohol-related GWAS54 across all cell types. For each gene we calculated posterior inclusion probability to both identify and/or confirm lead SNPs from the original study. We identified 223 credible lead SNPs (168 located within our open chromatin peak set) for 98 genes across cell types (Supplementary Data 35). We confirmed 52 SNPs previously reported (7 existing within our chromatin peak set) for 63 genes across cell types (Supplementary Data 36). Of the 131 GWAS DEGs, we found lead SNPs for 84 in INH. We identified credible lead SNPs for several genes implicated in neurotransmitter release and excitatory/inhibitory (E/I) balance including CACNA1C (rs12424245 PIP=0.39, P-value=1.29e−34) (Fig. 6c), SEMA6D (rs1045688 PIP=0.18, P-value=2.41e−13) (Fig. 6d), PPP1R13B (rs11160760 PIP=0.33, P-value=3.79e−40) (Supplementary Fig. 12a), and BCL11B (rs1405238 PIP=0.67, P-value=3.47e−57) (Supplementary Fig. 12b). We also confirm the lead SNP for DRD2 (rs61687445 PIP=0.50, P-value=6.58e−31) in INH neurons (Fig. 6e). In AST, of the 131 GWAS DEGs we found 74 lead SNPs. We highlight the gene AGBL2 in AST (Fig. 6f) where we identified a credible lead SNP (rs2293578, PIP=0.37, P-value=2.23e−21). Finally, of the 131 GWAS DEGs we found lead SNPs for 84 in MIC. We highlight the risk gene NF1 in MIC (Supplementary Fig. 12c) where we identified a credible lead SNP (rs2214538, PIP=0.22, P-value=1.12e−17). The rs2214538 SNP was the only one within a MIC ATAC peak and maintain significant predicted looping. By fine-mapping PAU risk variants within cell type-specific ATAC data with posterior probability, we have nominated or confirmed putative causal risk loci and highlight these as candidates for additional functional genomic follow up.

Methods
Human tissue donors
Human postmortem brain tissue samples were obtained from the National Center for PTSD Brain Bank (with consent of next of kin) and the University of Pittsburgh Tissue Donation Program. Brain specimens were obtained during autopsies conducted at the Allegheny County Medical Examiner’s Office (Pittsburgh, PA) after consent for donation was obtained from the next-of-kin. Individuals were a mix of European– and African–American descent. Males and females were group-matched for age and PMI. Sociodemographic and clinical details are listed in Supplementary Data 1 and include the manner of death, alcohol-related phenotypes, tobacco use, additional psychiatric disorders, and illicit drug use. Inclusion criteria for AUD and CON cases were as follows: PMI < 48 h, age range > 18 to < 85 years. A total of 50 individuals (22 AUD: 9 females, 13 males; and 28 control donors without AUD: 17 females, 11 males) were used in this study. To control for the high comorbidity of mood disorders in our AUD sample (100% are comorbid for either major depression or PTSD) we included additional donors (N=15, 53.5%) with no history of alcohol abuse or dependence but did have MDD or PTSD in our control group (control donors without AUD; Supplementary Data 1) to aid in identifying AUD-specific changes. The brain tissue was fresh-frozen and samples from the central nucleus of the amygdala79 were collected (~5 mg) from each postmortem sample.

Psychiatric history and demographic information were obtained by psychological autopsies performed postmortem as well as a review of medical records and toxicology reports. Trained clinicians conducted structured interviews with informants (usually the next of kin) with knowledge of the deceased individuals. To avoid systematic biases, AUD and control cases from each source were characterized by the same psychological methods. Consensus DSM-IV diagnoses for each subject were made by trained clinicians who did not conduct the psychological autopsies. Severity ratings are also based on DSM-IV.

All brain tissue underwent physical examination by a board-certified neuropathologist at both the macroscopic and microscopic level. All tissue was screened for confounding neuropathologies, including amyloid plaques, Lewy bodies (Parkinson’s disease and dementia with Lewy bodies), and trans active response DNA binding protein-43 (TDP-43) (ALS, frontotemporal dementia, and limbic predominant age-related TDP-43 encephalopathy). All cases were also evaluated at a gross (macroscopic) level to detect evidence of atrophy, trauma, or infarction. Importantly, cases were excluded if massive trauma occurred (from head injury) that severely damaged the tissue, if the tissue was damaged from a stroke involving a large portion of the brain, if the decedent was on a respirator for an extended period of time, if the donor had brain cancer, or if the donor had a history of HIV, AIDS, COVID-19, or other communicable disease.

Brain nuclei isolation for multiome assay
Regions of interest were dissected on cryotome (leaflets of ~100–300 micrometers) from frozen Central Amygdala and stored at −80°C. Cell nuclei isolation from brain sections were treated similarly to already established protocols80–82 with some modifications needed to have nuclei suitable for two separate assays: gene expression and ATAC. To avoid experimental bias, nuclei isolation was done by the same person blinded for the metadata for the particular sample. All reagents were molecular biology grade and sourced from Sigma unless stated otherwise. Small amounts of tissue (10–20 mg) were added into 1 ml of ice-cold lysis buffer (“Buffer A” is 250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 10 mM Tricine-KOH (pH 7.8), protease inhibitors w/o EDTA (Roche), RNAse inhibitor (80 U/ml) (Roche), 1 mM DTT, 1% BSA (m/v) (Gemini Bio-Products, Woodland, CA), 0.3% NP-40 (v/v), 0.15 mM spermine, 0.5 mM spermidine, water–DTT, RNAse Protector, protease inhibitors, spermine, spermidine and NP-40 were added immediately before use. The suspension was transferred to 2 ml Dounce tissue homogenizer (Kimble) and lysed with constant pressure and without introduction of air with pestle A (2×5) and pestle B (4×5) in cycles of 5. The homogenate was strained through pre-wetted 40 μm tubetop cell strainer (Corning). All subsequent centrifugation was performed in a refrigerated, bench-top centrifuge with swing-out rotor (Eppendorf, Hamburg, Germany). Lysates were centrifuged at 1000g, 10 min, 4°C, pellets were saved, and resuspended in 0.4 ml lysis buffer “Buffer A”. Final 0.4 ml of solution was mixed with 0.4 ml (1:1) of Optiprep solution (Buffer “B” is iodixanol 50% (v/v), 25 mM KCl, 5 mM MgCl2, 20 mM Tricine-KOH (pH 7.8), protease inhibitors w/o EDTA, RNAse inhibitor (80 U/ml), 1 mM DTT, water). The suspension (25% iodixanol final) was mixed 10x head over tail. The 2 ml tube was filled with 0.6 ml of 40% iodixanol cushion (appropriate mix of Buffer “A” and “Buffer “B”), then overlayed with 0.6 ml of 30% iodixanol (appropriate mix of Buffer “A” and “Buffer “B”), and sample suspension was overlayed over at the top. The tubes were then centrifuged at 3000g, for 30 min at 4°C without brakes. After centrifugation ended the interphase 30%/40% ring was collected (300 microliters). 500 microliters of “ATAC Lysis Buffer” was added to the sample, mixed well and left on ice for 5 min. “ATAC Lysis Buffer” is: 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.01% (v/v) Tween-20, 0.01% (v/v) NP-40, 0.001% (v/v) Digitonin, 1 mM DTT, RNAse inhibitor (80 U/ml) (Roche), 1% (m/v) BSA, Nuclease-free Water. After 5 min, 500 microliter “ATAC Wash Buffer” was added to ATAC sample, mixed well and left on ice for 5 min. “ATAC Wash Buffer” is: 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% (v/v) Tween-20, 1 mM DTT, RNAse inhibitor (80 U/ml) (Roche), 1% (m/v) BSA, Nuclease-free Water. The samples were centrifuged at 1000g, 10 min, 4°C, supernatants were carefully and completely removed. After last centrifugation the samples were resuspended in “Diluted Nuclei Buffer” (1X Nuclei Buffer (PN-2000207), 1 mM DTT, 1 U/microliter RNase Protector, water) counted on hemocytometer at the concentration of 3.3 million/mL and submitted to YCGA where the protocol was followed for targeting 10000 nuclei. Libraries were sequenced with paired-end 150bp reads on an Illumina NovaSeq 6000 to a target depth of up to 500 million read pairs per sample for gene expression (RNA), and up to 500 million reads per sample for ATAC.

Animals
All animal care and procedures were performed according to the ethical guidelines of the Institutional Animal Care and Use Committee at Yale University School of Medicine and Yale Animal Resources Center. Adult 2-month-old male and female C57BL/6 mice (n = 8 mice in total) were used for this study. Until treatment, all mice were group housed and maintained in standard environmental conditions (23°C; 12h–12h light–dark cycle) with ad libitum food and water.

Ethanol treatment
Mice were delivered at seven weeks of age and allowed two weeks to accommodate to the vivarium, where they were housed in a reverse light-dark cycle housing room. Two hours into the dark cycle, body weights were measured. Three hours into the dark cycle, mice received two sequential intraperitoneal injections each with a 10 ml/kg injection volume: the first was a solution of 5% DMSO in saline directly followed by 2 g/kg ethanol dissolved in water (n=2 M, 2 F) or saline (n=2 M, 2 F). Animals were randomly assigned to control or ethanol treated groups. Three hours following these injections, mice were euthanized by rapid decapitation and brains removed and immediately snap frozen on dry ice and then stored at −80°C. For tissue analyses, frozen brains were sliced into 2-mm coronal sections on a metal brain matrix (Plastics One, Roanoke VA) and amygdala samples were collected bilaterally with a 1-mm tissue punch (Fine Science Tools, Foster City, CA). Bilateral punches were collected from a coronal section targeted at −1.43 mm from bregma, according to 83. Nuclei were isolated from rodent amygdala tissue as described in Brain Nuclei Isolation for Multiome methods section.

Fluorescence Immunohistochemistry
Immediately following autopsy, tissue blocks for immunohistochemistry were cut and submerged in 4% paraformaldehyde for 48 hrs at 4°C, followed by cryoprotection in sucrose. Sections were cut on a microtome at 40 um and stored in cryoprotectant at −20°C until use.

Following washes in phosphate buffer, tissue sections were pretreated for antigen retrieval (0.01M sodium citrate for 45 min at 90°C), followed sequentially by incubation in 1% sodium borohydride (30 min), PB wash, incubation in 0.3% Triton X-100 (30 min), blocked in 20% normal donkey serum (1.5 h), and incubated at 4°C for 72 hr with the following primary antibodies: goat anti-Gad1 1:100 dilution (R&D Systems AF2086), rabbit anti-PENK 1:500 dilution (Atlas Antibodies HPA013138), anti-Rb 647; rabbit anti- FREM1 1:50 dilution (Proteintech 13086-1-AP). Secondary antibody incubations were overnight at 4°C with anti-goat Alexa fluor-488 or anti-Rb AlexaFluor-647 (Invitrogen, Life Technologies, Carlsbad, CA). Following washing, sections were incubated in DAPI and mounted with ProLong Glass Antifade reagent (Life Technologies, Carlsbad, CA). Z-stack images were collected through the depth of each section on an Olympus FV-3000 confocal scanning laser microscope using a 20×, 0.75 na objective.

Acetylcholinesterase staining
A modification of the method of Lim et al., 2004 was used for acetylcholinesterase staining to aid in delineating amygdala subregions84. Briefly, slide mounted sections adjacent to sections used for marker identification were incubated for 2 h in 0.12% acetylthiocholine iodide, 0.075% glycine, 0.05% cupric sulfate and 0.68% sodium acetate. Slides were washed 3 times and incubated in 0.77% sodium sulfide nonahydrate for 15 min. Following 3 washes, slides were exposed to 1% silver nitrate for 10 min and washed again. Slides were dehydrated in ethanol, cleared and mounted.

RNAScope in situ hybridization
Blocks of fresh frozen human amygdala from three AUD subjects and three controls were sectioned at 10 μm and stored at −80°C. RNAScope in situ hybridization was performed with the RNAscope Multiplex Fluorescent Detection Kit v2 (Cat# 323100, ACD, Newark, CA, USA) according to manufacturer’s instructions. For each subject, one brain section was stained. Briefly, tissue sections were fixed with 10% neutral buffered formalin (Cat# 15740, Electron Microscopy Sciences, Hatfield, PA, USA) for 30 min at RT, series dehydrated in ethanol, pretreated with hydrogen peroxide for 10 min, and treated with protease IV for 30 min. Sections were incubated with probes for PENK, GRM8, NCAM, CALN1, and TCF4 (Cat# 548301-C2, 563351, 421461, 1317671-C1, and 557411 ACD, Newark, CA, USA). Probes were fluorescently labeled with TSA dyes (TSA-Vivid 650, Cat# 7536, and TSA-Vivid 570, Cat# 7535, Tocris Biosciences, Bristol, UK). Nucleus were labeled with DAPI (4′,6-diamidino-2-phenylindole).

Fluorescent images were acquired in z-series using an Olympus FV3000 confocal microscope with 60× 1.42NA objectives. For each brain section, images were taken from 3 regions of interest within the central amygdala. All images for the same combination of probes were acquired with the same imaging settings and laser power intensities. Subsequent automated quantification was performed utilizing the dotdotdot algorithm for MATLAB85. DAPI stained nuclei, transcripts and lipofuscin were segmented from the background. The segmented RNA transcripts were masked with segmented lipofuscin before the colocalization of transcripts and nuclei were counted. Statistics were performed with Prism 10 (Graphpad).

Single-nucleus multiome processing
Count matrix generation. We ran Cell Ranger ARC (v2.0.2) with the hg38 genome as the reference to process the initial reads from the raw FASTQ files. Quality control and filtering. We used Signac86 (v1.11.0), a comprehensive R package for the analysis of single-cell chromatin and gene expression data, to preprocess each snMultiome sample. Signac creates a Seurat object containing both RNA and ATAC assays. For quality control, we filtered cells using the following parameters: nCount_ATAC>1000, nCount_RNA>200, and TSS.enrichment>2.

RNA – Doublet detection. Next, doublets were identified using a combination of two computational methods: Scrublet87 and DoubletDetection88 and removed from each sample. Sample aggregation. We merged high-quality cells and high-quality samples in Pegasus (v1.5.0)40, a python tool for analyzing transcriptomes of single cells, resulting in a total of 55 samples (Supplementary Fig. 13). After aggregating these samples into a single data object, we filtered cells based on the following criteria: at most 10% mitochondrial genes, at least 200 genes, and at least 500 UMIs. Mitochondrial, sex, and ribosomal genes were excluded, and only the robust genes were included in the final gene set. Five samples were further removed due to a surprisingly high EXC neuron composition, not expected in the central amygdala. The final snRNA-seq data object had 183,391 nuclei and 26,919 genes. Clustering and annotation. After cell and gene filtration, the data was transformed into logarithmic space by normalizing each cell by total counts. Next, the top 2000 highly variable features were found and principal component analysis (PCA) was applied on the highly variable features using the top 50 components. Batch correction was performed with Harmony. The top 50 components were then used to build a k-nearest-neighbors graph with k=100 neighbors, and Leiden clustering was used to identify cell clusters. After the first round of clustering, marker genes were used to inspect and remove clusters that had mixed or unmatched marker gene expression, and another round of clustering was performed. Cell type annotation was adapted using markers from Ma et al. with 7 cell types80. Further subtypes were discovered through either NSForest de Novo or confirmed via existing literature on the amygdala36–38. Most notably, we discovered 6 INH subtypes: DRD1, LHX9, SST, VIP, PENK, and FREM1.

ATAC—Quality control and filtering. Continuing from the 183,391 nuclei in the RNA QC step, we created Arrow files for each sample with cells filtered based on the quality control parameters filterTSS=4 and filterFrags=1000 in ArchR (v1.0.3)89, an R package for analyzing single-cell ATAC-seq data with multi-omic support. ArchR considers three characteristics: 1) the number of unique nuclear fragments, 2) the signal-to-background ratio, and 3) the fragment size distribution. Due to nucleosomal periodicity, we expect to see depletion of fragments that are the length of DNA wrapped around a nucleosome, around 147bp. We used the addDoubletScores function to infer potential doublets that can confound downstream results. Clustering and annotation. The addIterativeLSI function was used to implement iterative LSI dimensionality reduction. Cell type annotations were transferred from the labels in the RNA modality, leveraging the multi-omic nature of our dataset. We visualized the embedding using addUMAP and plotEmbedding. We produced a final embedding of 174,188 high-quality nuclei in both RNA and ATAC modalities (Supplementary Fig. 13, Supplementary Data 37).

Differential gene expression analysis
snRNA-based differential expression analysis was performed for each cell type and the 6 INH subtypes using four tests: 1) parametric-based MAST39 with covariate correction including covariates such as age, sex, ancestry, PMI, and RIN, 2) permutation-based Wilcoxon Rank Sum Test using the de_analysis function from Pegasus, 3) parametric-based MAST39 with covariate correction including covariates such as age, sex, ancestry, PMI, RIN, and additional alcohol clinical covariates 4) pseudobulk using DESeq241 (v1.46.0) with covariate correction. Only genes expressed in at least 5% of the cells, the robustly expressed genes, were included for these tests. The FDR significantly differentially expressed genes (DEGs) were those that passed the thresholds of fold change (FC)>1.2 and FDR<0.01 and overlapped in the same direction for MAST and Wilcox. For INH, we used the FC value of the subtype that had the greatest magnitude for analyses. These DEG lists were utilized in downstream analyses such as TF-CRE-DEG links.

Peak calling on snATAC clusters
Pseudo-bulk replicates were created for each cell type with the addGroupCoverages function to resolve inherent sparsity in the snATAC-seq data. MACS248 was then used to call peaks on each pseudo-bulk replicate. ArchR uses an iterative overlap peak merging procedure with the addReproduciblePeakSet function to create a final peak set of fixed-width peaks using the following procedure: 1) peaks are first ranked by their significance, 2) the most significant peak is retained and any peak that directly overlaps with the most significant peak is removed from further analysis, and 3) of the remaining peaks, this process is repeated until no more peaks exist. ArchR analyzes all the pseudo-bulk replicates from a single cell type together, performing the first iteration of iterative overlap removal. It then checks to see the reproducibility of each peak across pseudo-bulk replicates and only keeps peaks that pass a threshold indicated by the reproducibility parameter. This resulted in a single merged peak set for each of the seven cell types. This procedure is then repeated to merge the cell type-specific peak sets by renormalizing the peak significance across the different cell types and performing the iterative overlap removal, resulting in a single merged peak set of fixed-width (501bp) peaks, denoted the union peaks. We utilized both cell type-specific peaks and the union peaks depending on the analysis we performed. The addPeakMatrix function adds the Peak Matrix using the union peaks to the data object, which is used for downstream analyses such as finding peak-to-gene links and performing TF motif enrichment analysis.

Characterization of CREs
Since snMultiome already simultaneously profiles the gene expression and chromatin accessibility of the same nuclei, no integration between the two modalities was needed. Directly using the GeneExpression Matrix, we identified peak-to-gene links using the addPeak2GeneLinks function, which leverages integrated snRNA-seq data to look for correlations between peak accessibility and gene expression. We set maxDist=1,000,000 to consider long-range interactions. Then we identified putative enhancers or cis-regulatory elements (CREs) based on the thresholds: peak-to-gene correlation > 0.4 and FDR<0.05.

TF motif enrichment analysis and TF regulatory network construction
TF motif enrichments were calculated to predict which regulatory factors are most active in a given cell type. The addMotifAnnotations function indicates motif presence in the union peaks using a binary matrix, and the peakAnnoEnrichment function tests either cell type-specific marker peaks or differential peaks for enrichment of various motifs. Besides the previous TF motif enrichments, the R package chromVAR predicts enrichment of TF activity on a per-cell basis while controlling for known technical biases. The addDeviationsMatrix function, which adds the Motif Matrix to the data object, was used to compute deviation z-scores for each motif by comparing the number of fragments that map to peaks containing the motif to the expected number of fragments in a background peak set that accounts for confounding factors such as GC content bias, PCR amplification, and Tn5 tagmentation. Then we used the getPositions, addGroupCoverages, and getFootprints functions to perform TF footprinting analysis with pseudo-bulk aggregates of cells in the same cell type. In the plotFootprints function, ArchR performs normalization by subtracting the Tn5 bias from the footprinting signal. We plot the footprinting results for TFs RFX2 and KLF16 comparing AUD and CON in INH neurons. We then created TF regulatory networks based on the TF-CRE-DEG links. We identified several TFs with high enrichment and plotted their target genes using the R package igraph (v1.2.6), where each vertex represents a TF or target gene, and each edge represents a linked enhancer binding event, overlaying additional information onto the network such as DEG (up/down) or AUD GWAS gene (True/False) status. For the TF networks, we used more conservative correlation cutoffs of 0.9 for visualization purposes.

Disease-specific peak-to-gene calculation
To combat the sparsity of peak data, which makes it difficult to calculate disease-specific differences, we focused on peak-to-gene linkages for our multi-omic data. Cell-type-specific peak-to-gene links were first calculated for each of the 7 canonical cell types via the cellsToUse parameter in the addPeak2GeneLinks function. Then, for each cell type, peaks with strong peak-to-gene correlation in AUD (corr_AUD>0.5) and also a large difference between AUD and Control (corr_diff>0.5) were kept. Finally, motif enrichment was performed on these subsetted peaks for each cell type to highlight disease-specific motifs.

Estimating GWAS enrichment using cell type-specific peaks
To estimate heritability of a variety of complex traits, we used LDSC (v1.0.1)90. GWAS summary statistics for Problematic Alcohol Use/Alcohol Use Disorder, psychiatric, and non psychiatric disorders were correctly formatted with munge_sumstats.py and lifted to hg38 coordinates using UCSC liftover. Cell type-specific peaks, peaks from the disease-specific peak-to-gene calculation, and INH peaks separated by genomic region were formatted for LDSC using make_annotation.py, and LD scores were computed for each set using ldsc.py. Benjamini-Hochberg multiple-testing correction was applied to the enrichment P-values.

GWAS fine-mapping
Peak enrichment scores as prior weights. We performed genetic fine-mapping with the sum of single-effects (SuSiE) regression model53 to infer risk variants with epigenetic information from our single-nuclei ATAC peaks. With some prior information on SNPs, we can improve the accuracy of detecting the causal variants91. For SNPs inside ATAC peaks, we used the peak enrichment scores (−log10 pval) which indicate the strength of each peak, as the prior weights for SuSiE to prioritize SNPs with strong ATAC signals. For SNPs outside ATAC peaks, we set the prior weights to 0.1. We then subset the union peaks to the CREs (peak-to-gene correlation > 0.4 and FDR<0.05) to perform fine-mapping in SuSiE. We performed SuSiE regression with the susie_rss function using AUD GWAS summary statistics (PAU EUR)16. SuSiE calculates the posterior inclusion probability (PIP), the probability for a given SNP being causally associated with the trait of interest, for each variant. The LD information was derived from the 503 European samples of the 1000 Genomes Project. Input files. We used the 1000 Genomes Project (1KG) Phase 3 dataset92 as the reference panel for our analysis. This dataset can be accessed via the Resources section on the PLINK 2.0 website at Resources-PLINK2.0 (cog-genomics.org). We selected the European samples based on the super-population information provided by 1KG and excluded all duplicate and ambiguous SNPs. We applied quality control to the 1KG data with European ancestry using PLINK93. To be more specific, we only included SNPs with a 99% genotyping rate and the minor allele frequency (MAF) no less than 0.005. We also exclude samples with more than 5% missing genotypes and markers that failed to pass the Hardy-Weinberg test (PLINK command: --geno 0.01 –hwe 1e-10 –mind 0.05 –maf 0.005). The final reference panel includes 503 non-overlapping European samples, genotyped at 8,190,311 SNPs. To obtain the coordinates for risk SNPs, we turned to the UCSC Genome Browser94,95 and focused on the human genome version GRCh38/hg38.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Lee et al., reports an interesting and important multi-omic single nucleus profiling of postmortem human\namygdala from donors with alcohol use disorder (AUD). The study aims to identify cell-type specific dysregulated gene\nregulatory programs. The approach is novel, and the study could improve our understanding of the molecular underpinnings\nof AUD; however, the design, the analysis and the reporting contain significant shortcomings that should be addressed.\nInformation on the metadata of the donors are significantly lacking, the authors do not report on the cause of death, lifetime\nconsumption of alcohol, co-use of other substances of abuse (other than tobacco), presence or absence of psychiatric\ndisorders, and related medications. All of which should have been included as covariates in the differential expression\nanalyses.\nA major limitation of the study is the lack of validation either in bulk tissues or alternative methods confirming cell specificity.\nSingle cell transcriptome approaches are strengthened by including some experimental validations since many of the\nfindings are computationally based.\nThere is a discrepancy between the reported method and the nuclei isolation protocol, 10x Genomics Chromium Next GEM\nSingle Cell Multiome ATAC + Gene Expression profiling which simultaneously profiling ATAC and RNA sequencing in the\nsame nucleus and the reported nuclei isolation method implies producing two aliquots from each sample one of which is\nprocessed for single nucleus ATAC and one for gene expression profiling.\nThe profiling yielded 183,391 cells from 50 samples (average 3600 nuclei per sample), the process as described targets\n10,000 cells. Details should be provided whether the nuclei were eliminated during quality control, and report on quality\ncontrol metrics for individual samples, or cells were lost during the library preparation process. How many cells were\nrecovered from each sample?\nHow many cells were taken out with Scrublet and DoubletDetection? ~183,000 cells is quite low for 55 samples. What does\nthe barcode rank plot look like across all samples, authors denote high and low quality but don’t show QC metrics for these\nsamples.\nThere is a discrepancy in reporting of the Wilcoxon test (de_analysis function from Pegasus\nin the methods and Seurat in the discussion section). The scripts used for performing MAST was not available at the time of\nreview since link for the Github is not working. Please confirm active links to provide access to the scripts used in the study.\nMAST and wilcoxon rank sum testing are quite lenient differential methods in single cell studies. This is evident in the higher\nfold change and FDR cutoffs used in the study. Were these directly compared to pseudo bulk methods? Sample size can\nhave unexpected impacts on most current DEG methods. Using two methods and then overlapping does help ensure that\nwhat DEGs are being identified are consistent but it would be good see pseudo bulking as well to determine the overlap\nbetween single nuclei based methods and a separate, more conservative approach for DE.\nP values should be reported for proportional change showing no difference in cell types between AUD and Control.\nIt is stated that microglia and macrophages were found but do not show markers to differentiate and collapse them further in\nthe paper. F13a1 is not exclusive to macrophages and can be representative of monocyte-derived infiltrating macrophages\nand in some cases microglia. Cd44, CD169 and Siglech have been reported as good marker for infiltrating macrophages\nrather than resident microglia, in addition to higher expression of CD9. What is the expression of a more robust microglia\nmarker like TMEM119 or P2RY12? The identification of a true macrophage population is not possible without extensive\nmarkers and reporting the differences in key microglia homeostatic genes.\nThe authors report sex specific AUD changes and provide supplementary tables for male AUDvsCON and female\nAUDvsCON; however, there was no discussion of how the analysis were performed and specifics were not provided. Given\nthe sample size this analysis is likely underpowered but potential sex differences should be examined and noted in the text.\nAuthors discuss the changes in chromatin accessibility that were linked into differentially expressed genes, these are\npotentially important findings, and additional details should be provided discussing instances where the changes of\nchromatin accessibility did not match observed gene expression changes and vice versa. Additional discussion of these\nfindings is warranted.\nAdditional detail and discussion about the changes in cell types other than inhibitory neurons should be provided and any\noverlaps between the identified differentially regulated elements in these different cell types. This could be presented in the\nform of a potential supplemental table.\nThe authors note a nominally significant up regulation of OPRM1 in INH neurons. This finding is somewhat surprising since\nOPRM1 is often not dysregulated in transcriptome studies specifically examining OUD. The donor metadata should be\nexamined for comorbid opioid use or reported OUD/MMD medication usage that could potentially explain this interesting\nfinding.\nMinor comments\nLine 96: (mean reads/nuclei), how it is written it suggests that this was the target sequencing depth.\nFigure 4b is not referenced in the text", "reviewer_2": "The manuscript by Lee et al provides insight into alcohol use disorder in humans by investigating single cell chromatin and\ngene expression data from post-mortem brain donors. The paper is well written and provides new data to the field. I’m not\nfully comfortable to provide an in-depth review of the bioinformatic analyses and therefore focus on clinical aspects primarily.\nThe title is slightly incorrect because the technique used is single nucleaus multiome.\nI’d suggest including the sample size in the abstract.\nLine 89: A quick google search identifies snRNAseq studies in human CeA, I may be wrong.\nThe clinical spectrum includes individuals with alcohol abuse and alcohol dependence, two very different phenotypes. What\nis the rationale and justification to combine both beyond the notion that DSM-V collapsed those two (can they be combined\nin this dataset given that the clinical heterogeneity within AUD is large)? The donor’s age at death is young, what led to their\ndemise? Do you have information on the severity of AUD (mild, moderate, severe) and the timeline across life? Did donors\nconsume alcohol at death (neurotox data?)? Do you have information on the neuropathological assessments of the donors\n(cases and controls), especially for those >60yrs [neurodegen panel]?\nLooking at the Extended Data Fig 1a, the details given regarding the post-mortem samples is very limited. The U Pitt Brain\nBank has a large amount of data on their donors that is critical to include and to evaluate the data and this manuscript. For\nexample: Co-morbidities (anxiety disorders, depression, psychosis), somatic disorders (hepatic function, malnutrition etc),\ntimelines (AUD as lifetime diagnosis or at death, quantity of alcohol consumed, manner of death, cause of death, medication\ntaken, etc etc.) It is critical to look at the clinical variables in depth as this is not a well-controlled animal study. In line 425 the\nauthors refer to ED Fig 1a for cause and manner of death, which is not there. Smoking is the only comorbidity listed here and\nmay be significantly different between groups (age, PMI as well?).\nThe number of DEGs is high for the sample size and I think the approach to use two different methods is good. However, can\nthe authors provide QQ plots or other data that may support the notion that these results are not inflated? The inclusion of\ncovariates in the MAST approach is basic. I understand that including a myriad of clinical covariates is difficult given the\nsample size, however, I wonder if a surrogate variable analysis could cover those known and unknown factors. Although the\nstudy here is the largest up to date on AUD, an independent replication would be critical due to the limit in sample size.\nGiven that the results presented are primarily associative, some functional manipulation of genes or cell types in animal\nmodels would be critically important to better understand the generalizability of the data."}, "author_response": "Response Letter for “Central amygdala single-nucleus atlas\nreveals chromatin and gene transcription dynamics in human\nalcohol use disorder”\nOverview of our Response:\nWe thank each of the reviewers for their time, effort, and constructive feedback on our revision\nof the manuscript. Below, we respond to each of the reviewer’s specific comments and\nconcerns in a point-by-point fashion. To summarize, we have conducted:\n● Comparative analysis of different DEG models: We compared the DEG results for AUD vs\nCON on our snRNA-seq data using different methods: 1) our original MAST+Wilcox\nintersection; 2) pseudobulk approach with DESeq2 and 3) a MAST approach using additional\nclinical covariates (including additional psychiatric disorders, drug use, etc.). Similar GO\nterms were identified between DEGs from MAST+Wilcox and those from DESeq2. In\naddition, we found that >80% of DEGs overlapped between our original model and those\nincluding additional clinical covariates.\n● Confirmation of Gene Expression Changes: We now include a series of fluorescent in situ\n(FISH) experiments for our top AUD GWAS DEGs. Small molecule FISH was used for GRM8,\nCALN1, TCF4, and NCAM1 to confirm expression changes within INH subtypes.\n● Additional Differential Gene Validation and Functional Analysis: As suggested, we have now\nadded an animal model of alcohol use and performed the same molecular analyses as the\nhuman tissue including (snMultiome including RNA and ATAC) after mice were treated with\na 2 g/kg injection of 100% EtOH.\n● Expansion of sex-specific analyses: We have expanded our sex-specific analyses. We\nperformed three sets of sex-specific DEG analyses and included them in new Extended\nData Fig. 4. We now report Male AUD DEGs (Supplementary Data Table 8), Female AUD\nDEGs (Supplementary Data Table 9), Male vs Female Control differences\n(Supplementary Data Table 10) and expand on these findings in the text.\n● Expansion of cell type analyses: We expanded our results and discussion to include more\nmaterial about other cell types with genomic changes. We specifically look at astrocytes\nGiven the extensive number of DEGs and CREs we identified, we also highlighted\nexpression changes, P2G linkages, and cell-type-specificity in 2 genes of astrocytes in\nExtended Data Fig. 7g):\n● Changes to title and main text: As suggested, we have changed our current title to “Central\namygdala single-nucleus atlas reveals chromatin and gene transcription dynamics in human\nalcohol use disorder”.\n● Code update: We have uploaded our code on https://github.com/mjgirgenti/AUDsnCEA.\nRef 1.0 General Comments\nReviewer’s The manuscript by Lee et al., reports an interesting and important\ncomment multi-omic single nucleus profiling of postmortem human amygdala\nfrom donors with alcohol use disorder (AUD). The study aims to\nidentify cell-type specific dysregulated gene regulatory programs.\nThe approach is novel, and the study could improve our\nunderstanding of the molecular underpinnings of AUD; however, the\ndesign, the analysis and the reporting contain significant\nshortcomings that should be addressed.\nAuthor’s We thank the reviewer for their review and have made a point-by-\nresponse point response addressing their comments. Specifically, we have\n1. Added new forms of experimental validations (ie. Mouse\nalcohol addiction model and fluorescent in situ of top DEGs)\nthat validates many of our original conclusions.\n2. Added new DEGs analyses that account for sex, non-INH\ncell types, pseudobulk analysis, and new clinical covariates.\n3. Included new plots and tables clarifying our original results\n4. Shared our GitHub link and edited the manuscript\nPlease see the detailed answers to your question below.\nRef 1.1 DEG Covariate\nReviewer’s Information on the metadata of the donors are significantly lacking,\ncomment the authors do not report on the cause of death, lifetime\nconsumption of alcohol, co-use of other substances of abuse (other\nthan tobacco), presence or absence of psychiatric disorders, and\nrelated medications. All of which should have been included as\ncovariates in the differential expression analyses.\nAuthor’s We thank the reviewer for the opportunity to clarify, and we agree\nresponse that these covariates are important for our analysis. As suggested,\nwe have added all the above metadata to new Supplementary\nData Table 1. We have also included these in a new DEG\nanalysis:\n1. Included the new metadata in Supplementary Data Table 1\n2. Performed a new DEG analysis that includes Manner of\nDeath, Lifetime Alcohol Consumption, Illicit Drug Use, Other\nPsychiatric Disorders, and Meds ATOD as our new\ncovariates. These new DEGs are included in\nSupplementary Data Table 6.\n3. We found that there was 80% overlap with DEGs discovered\nin our original model (Extended Data Fig. 4a)\nExtended Data Fig. 4a. Overlap of DEG analysis with new\ncovariates.\nExcerpt from To assess the impact of additional drug use and comorbidities, we\nthe manuscript performed additional analyses that included covariates for manner\nof death, lifetime alcohol consumption, illicit drug use, comorbid\npsychiatric disorders, and medications at time of death\n(Supplementary Data Table 6). Across cell types, ~80% of DEGs\noverlapped with our model (Extended Data Fig. 4a).\nRef 1.2 Validation\nReviewer’s A major limitation of the study is the lack of validation either in bulk\ncomment tissues or alternative methods confirming cell specificity. Single cell\ntranscriptome approaches are strengthened by including some\nexperimental validations since many of the findings are\ncomputationally based.\nAuthor’s We thank the reviewer for this suggestion and agree that further\nresponse experimental validations will strengthen these findings. As\nsuggested.\nFirst, we performed fluorescent in situ hybridization (FISH). FISH\nallows for confirmation of expression changes within particularly cell\n3\ntypes and compartments. We chose to focus on the GWAS DEGs\n(Fig. 3e) as these seem to be most pertinent to AUD biology\nthrough their link with genetic risk. We designed probes for 4\ntranscripts: GRM8, TCF4, CALN1, and NCAM1. The probe\ndesigned for GABRA2 failed, unfortunately. We were able to\nconfirm expression changes for each of these genes. These are\nnow included in Extended Data Fig. 5c-j and are also included\nbelow for the reviewers.\nc, Example images of FISH from control (left) and AUD (right) tissue\nprobed for PENK (red) and NCAM1 (white). d, Average NCAM1 pixel\nnumber per PENK+ nucleus. Mann-Whitney test, Control: n = 24 nuclei,\nAUD: n = 29 nuclei. e, Example images of FISH from control (left) and\nAUD (right) tissue probed for PENK (red) and GRM8 (white). f, Average\nGRM8 pixel number per PENK+ nucleus. Unpaired t test, Control: n =10\nnuclei, AUD: n = 9 nuclei. g, Example images of FISH from control (left)\nand AUD (right) tissue probed for CALN1 (white). h, Average CALN1 pixel\nnumber per nucleus. Unpaired t test, Control: n = 103 nuclei, AUD: n = 91\nnuclei. i, Example images of FISH from control (left) and AUD (right)\ntissue probed for TCF4 (white). j, Average TCF4 pixel number per\nnucleus. Mann-Whitney test, Control: n = 189 nuclei, AUD: n = 119 nuclei.\n3 control donors and 3 AUD donors were used for these analyses. 3 ROIs\nwere analyzed per donor. P values indicated on violin plots. Green:\nLipofusion. Blue: DAPI.\nperformed the same molecular analysis as the human tissue\n(snMultiome including RNA and ATAC) after mice were treated with\nan 2 g/kg injection of 100% EtOH.\n1. Mouse alcohol model validates CeA cell types: INH, EXC,\nOLI, OPC, END, AST, MIC on snRNA, snATAC, and\nsnMultiome embeddings, respectively (Extended Data Fig.\n9a-c).\nExtended Data Figure 9: Mouse snMultiome data. UMAP\nvisualization across seven cell types of a, snRNA-seq b, snATAC-\nseq across seven cell types and c, snMultiome\n2. Mouse alcohol abuse model snMultiome-seq validates\nDEGs: We identified 1814 total DEGs across all major cell\ntypes with most falling within INHa (Extended Data Fig. 9d).\nWhen compared with the human AUD DEG dataset we found\noverlap for 678 DEGs (Extended Data Fig. 9e).\nExtended Data Figure 9: Mouse snMultiome data. d, Significant\nDEG counts in both directions for each major cell type. e, Intersection\nof alcohol DEGs between human and mouse. f, Top 25 enrichR GO\nterms of the 1,814 unique AUD DEG set.\n3. Mouse alcohol model snMultiome-seq to validates Chromatin\nPeaks: We performed peak-to-gene linkage analysis on our\nmouse snMultiome data and identified 7485 total in the\n5\nmouse AUD model. Of these, 4710 overlapped with the\nhuman AUD CRE genes (62.9%) (Extended Data Fig. 9h).\nExtended Data Figure 9: Mouse snMultiome data. h, Intersection\nof alcohol-focused Cis-Regulatory Element–linked genes between\nhuman and mouse.\n4. Mouse alcohol model snMultiome-seq to validate Gene\nRegulatory Networks. We found that combined KLF16, KLF 7\nand KLF6 form a GRN that regulates 549 DEGs. We found\nincluding DEG-GWAS gene CALN1 confirming its role in\nAUD.\nExtended Data Figure 9: Mouse snMultiome data. i, TF regulatory\nnetwork showing the target genes for TFs KLF6, KLF7, and KLF16 in\nINH. Peak-to-gene correlation > 0.9 was employed in visualizing the\nnetwork.\nthe manuscript mouse model of alcohol abuse. We isolated frozen amygdala from\nmice treated with an 2 g/kg injection of 100% EtOH. We performed\nsingle nuclei isolation and snMultiome in the same manner as our\nhuman donor cohorts. We identified all major cell types across all\nmodalities: RNA (Extended Data Fig. 9a), ATAC (Extended Data\nFig. 9b) and combined Multiome (Extended Data Fig. 9c) for 19,072\nindividual nuclei. We identified 1814 total DEGs across all major cell\ntypes with most falling in the INH (Extended Data Fig. 9d). When\ncompared with the human AUD DEG dataset we found overlap for\nin neuronal cell types and this was reelected in gene set enrichment\nanalysis which identified terms related to neuron projection, synaptic\ntransmission, and glutamatergic signaling (Extended Data Fig. 9f).\nIn addition, we intersected our INH DEGs with a previous central\namygdala mouse model of alcohol withdrawal46 and found 38\noverlapping DEGs in INH and EXC neurons and MIC (Extended\nData Fig. 9g). We performed peak-to-gene linkage analysis on our\nmouse snMultiome data and identified 7485 total in the mouse AUD\nmodel. Of these, 4710 overlapped with the human AUD CRE genes\n(62.9%) (Extended Data Fig. 9h). To confirm the function of our KLF\n7\ngene regulatory network in human AUD (Fig. 5g), we constructed a\nGRN for these same TFs in our mouse model (Extended Data Fig.\n9i). We found that combined KLF16, KLF 7 and KLF6 form a GRN\nthat regulates 549 DEGs. We found 50 DEGs overlapped between\nthe human and mouse GRN including DEG-GWAS gene CALN1\nconfirming its likely role in alcohol abuse. Taken together, these data\nsupport the functional significance of KLF TFs in regulating calcium\nand glutamatergic signaling pathways identified in our dataset.\nEthanol Treatment (from Methods)\nMice were delivered at seven weeks of age and allowed two weeks\nto acclimate to the vivarium, where they were housed in a reverse\nlight-dark cycle housing room. Two hours into the dark cycle, body\nweights were measured. Three hours into the dark cycle, mice\nreceived two sequential intraperitoneal injections each with a 10\nml/kg injection volume: the first was a solution of 5% DMSO in\nsaline directly followed by 2 g/kg ethanol dissolved in water (n = 2M,\n2F) or saline (n = 2M, 2F). Animals were randomly assigned to\ncontrol or ethanol treated groups. Three hours following these\ninjections, mice were euthanized by rapid decapitation and brains\nremoved and immediately snap frozen on dry ice and then stored at\n-80ºC. For tissue analyses, frozen brains were sliced into 2-mm\ncoronal sections on a metal brain matrix (Plastics One, Roanoke\nVA) and amygdala samples were collected bilaterally with a 1-mm\ntissue punch (Fine Science Tools, Foster City, CA). Bilateral\npunches were collected from a coronal section targeted at -1.43 mm\nfrom bregma, according to76. Nuclei were isolated from rodent\namygdala tissue as described in Brain Nuclei Isolation for Multiome\nmethods section.\nRef 1.3 Protocol Clarify\nReviewer’s There is a discrepancy between the reported method and the nuclei\ncomment isolation protocol, 10x Genomics Chromium Next GEM Single Cell\nMultiome ATAC + Gene Expression profiling which simultaneously\nprofiling ATAC and RNA sequencing in the same nucleus and the\nreported nuclei isolation method implies producing two aliquots from\neach sample one of which is processed for single nucleus ATAC\nand one for gene expression profiling.\nAuthor’s We thank the Reviewer for the opportunity to clarify our isolation\nresponse protocol. We have updated the method details in the manuscript.\nThere was only one aliquot of nuclei split for 2 assays of\nsnMultiome - RNA and ATAC.\nExcerpt from Regions of interest were dissected on cryotome (leaflets of ~100-\nthe manuscript 300 micrometers) from frozen Central Amygdala and stored at -80\nºC. Cell nuclei isolation from brain sections were treated similarly to\nalready established protocols with some modifications needed to\nhave nuclei suitable for 10x Genomics Chromium Next GEM Single\nCell Multiome ATAC + Gene Expression. To avoid experimental\nbias, nuclei isolation was done by the same person blinded for the\nmetadata for the particular sample. All reagents were molecular\nbiology grade and sourced from Sigma unless stated otherwise…\nRef 1.4 Cell Count 1\nReviewer’s The profiling yielded 183,391 cells from 50 samples (average 3600\ncomment nuclei per sample), the process as described targets 10,000 cells.\nduring quality control, and report on quality control metrics for\nindividual samples, or cells were lost during the library preparation\nprocess. How many cells were recovered from each sample?\nAuthor’s We thank the reviewer for the opportunity to clarify. Some cells are\nresponse lost due to our rigorous QC. We have now included a new\nSupplementary Figure 3 showing the total number of nuclei\n(before and after QC) by sample. We outline the number of nuclei\nremoved in each step below:\n1. Starting with 618,861 cells, we removed the samples which\nhad less than 2500 cells (612,963 cells).\n2. We manually checked each sample on three different\nUMAPs, RNA-ATAC-Multiome specifically, and filtered for\n“nCount_ATAC > 1000, nCount_RNA > 200, and\nTSS.enrichment > 2” via Signac (533,337 cells).\n3. We removed doublets – 63,268 cells were taken out with\nScrublet, and 23,070 cells were taken out with\nDoubletDetection (446,999 cells).\n4. After pooling all samples together, we filtered for cells that\nhave “at most 10% mitochondrial genes, at least 200 genes,\nand at least 500 UMIs” via Pegasus (357,551 cells).\n5. We then iteratively clustered on the RNA UMAP embedding\nfor 10+ iterations and manually removed streaks of doublets\nthat were not detected by Scrublet and DoubletDetection\n(233,591 cells).\n6. Moving on the ATAC modality, we filtered for cells that had a\ntranscription start site (TSS) enrichment score of at least four\nand at least ten mapped ATAC-seq fragments “filterTSS=4\nand filterFrags=1000 in ArchR” (212,435 cells).\n7. Finally, we removed five samples that most likely had “white\nmatter” and other amygdala subregion contamination\n(174,188 cells).\nThe recovered cell count for the final set of samples can be found in\nthe new Supplementary Data Table 37 and Supplementary Data\nFig. 3a. The excerpted 2 lines below are from the methods.\nafter data preprocessing.\nExcerpt from The final snRNA-seq data included 183,391 nuclei and 26,919\nthe manuscript genes.\nWe produced a final embedding of 174,188 high-quality nuclei in\nboth RNA and ATAC modalities (Supplementary Data Figure 3,\nSupplementary Data Table 37).\nRef 1.5 Barcode Rank\nReviewer’s How many cells were taken out with Scrublet and\ncomment DoubletDetection? ~183,000 cells is quite low for 55 samples. What\ndoes the barcode rank plot look like across all samples, authors\ndenote high and low quality but don’t show QC metrics for these\nsamples.\nAuthor’s We thank the Reviewer for the opportunity to clarify:\nresponse ● 63,268 cells were taken out with Scrublet, and 23,070 cells\nwere taken out with DoubletDetection. While this number is\nlower than might be expected by our n, the number of nuclei\n10\nanalyzed is consistent with and in some cases larger than\nother postmortem brain genomic studies\n● The added barcode rank plots (before and after) demonstrate\nthat only high-quality nuclei were retained (Supplementary\nData Fig. 3b). Looking at the barcode rank plots, it contains\nonly the high quality ‘Cliff’ portion and not the low quality\n‘Knee’ background.\nSupplementary Data Fig. 3b. Barcode rank plot across all\nsamples, before and after data preprocessing.\n● We determined high-quality samples from low-quality\nsamples by inspecting all 55 UMAPs individually for each\ngenomic modality, eliminating samples which did not cluster,\nresulting in a low cell count. Fig. R1 has examples of\nsamples which did not cluster.\nUMAPs across all modalities\nExcerpt from We merged high-quality cells and high-quality samples in Pegasus\nthe manuscript (v1.5.0)38, a python tool for analyzing transcriptomes of single cells,\nresulting in a total of 55 samples (Supplementary Data Figure 3).\nWe produced a final embedding of 174,188 high-quality nuclei in\nboth RNA and ATAC modalities (Supplementary Data Figure 3,\nSupplementary Data Table 37).\nRef 1.6 Method Clarify\nReviewer’s There is a discrepancy in reporting of the Wilcoxon test\ncomment (de_analysis function from Pegasus in the methods and Seurat in\nthe discussion section). The scripts used for performing MAST was\nnot available at the time of review since link for the Github is not\nworking. Please confirm active links to provide access to the scripts\nused in the study.\nAuthor’s We thank the Reviewer for bringing this to our attention.\nresponse ● We have corrected the writing in the discussion because we\nused the Wilcoxon test from Pegasus throughout.\n12\n● Additionally, the Github link\n(https://github.com/mjgirgenti/AUDsnCEA) is now active (Fig.\nR2).\nFigure R2. Screenshot of the working Github page.\n● The Github link\n(https://github.com/mjgirgenti/AUDsnCEA/blob/main/scRNA/\nDEG_MAST.R) also has the MAST DEG script (Fig. R3).\nFigure R3. Screenshot of MAST DEG script on the Github page.\nthe manuscript both the Pegasus based Wilcoxon Test and the more conservative\nMAST method (including covariates).\nRef 1.7 DEG PseudoBulk\nReviewer’s MAST and wilcoxon rank sum testing are quite lenient differential\ncomment methods in single cell studies. This is evident in the higher\nfold change and FDR cutoffs used in the study. Were these directly\ncompared to pseudo bulk methods? Sample size can have\nunexpected impacts on most current DEG methods. Using two\nmethods and then overlapping does help ensure that what DEGs\nare being identified are consistent but it would be good see pseudo\nbulking as well to determine the overlap between single nuclei\nbased methods and a separate, more conservative approach for\nDE.\nAuthor’s We thank the reviewer for this suggestion. Pseudobulk DEGs are a\nresponse more conservative approach and are becoming more popular in\ngenomic studies. We were able to find additional DEGs with the\npseudobulk DESeq239 (adjusted p-value < 0.05) method: 14 in\ninhibitory neurons, 26 in excitatory neurons, 1 in oligodendrocyte,\n10 in astrocytes, and 1 in microglia.\n● We show the overlap of pseudobulk and single-cell DEG lists\nbelow and they are included in Extended Data Fig. 4b.\nExtended Data Fig. 4b. Single-cell and PseudoBulk DEGs overlap\n● We have also included the Pseudobulk DEG list in new\nSupplementary Data Table 7.\nExcerpt from Additionally, to ensure rigor, we tested a sample-based DEG\nthe manuscript method (i.e. pseudobulk) using DESeq239 (Supplementary Data\nTable 7). We identified additional DEGs: 14 in inhibitory neurons,\n26 in excitatory neurons, 1 in oligodendrocyte, 10 in astrocytes, and\n1 in microglia (Extended Data Fig. 4b).\nReviewer’s P values should be reported for proportional change showing no\ncomment difference in cell types between AUD and Control.\nAuthor’s We thank the reviewer for the opportunity to clarify, and we have\nresponse since:\n1. Added p-values into our paper as suggested.\n● (INH) 0.718;\n● (EXC) 0.324;\n● (OLI) 0.977;\n● (OPC) 0.406;\n● (END) 0.551;\n● (AST) 0.945; and\n● (MIC) 0.945.\n2. We include Fig. R4 for the reviewer to visually inspect the\ncell type proportion per cohort for each cell type.\nFigure R4. Cell count per sample for each cell type\nExcerpt from We calculated cell type proportions based on the number of nuclei\nthe manuscript isolated (Fig. 1e) and found no significant proportion shifts due to\ndisease (AUD versus CON): INH P-value=0.718; EXC P-\nvalue=0.324; OLI P-value=0.977; OPC P-value=0.406; END P-\nvalue=0.551; AST P-value=0.945; and MIC P-value=0.945.\nReviewer’s It is stated that microglia and macrophages were found but do not\ncomment show markers to differentiate and collapse them further in\nthe paper. F13a1 is not exclusive to macrophages and can be\nrepresentative of monocyte-derived infiltrating macrophages and in\nsome cases microglia. Cd44, CD169 and Siglech have been\nreported as good marker for infiltrating macrophages rather than\nresident microglia, in addition to higher expression of CD 9. What is\nthe expression of a more robust microglia marker like TMEM119 or\nP2RY12? The identification of a true macrophage population is not\npossible without extensive markers and reporting the differences in\nkey microglia homeostatic genes.\nAuthor’s We thank the reviewer for the suggestion. We have now checked\nresponse our cell subtypes and integrated our original cell-subtyping with the\nreviewer’s markers.\nF13A1 was established as a macrophage marker by Ma et al. in the\n202277, and it emerged as a subpopulation marker de novo via the\nNSForest analysis in our data after we identified microglia using the\nestablished marker APBB1IP (Methods; Figure R6). The high co-\nexpression of F13A1 and SIGLEC1 would suggest that this is a\nmacrophage cluster.\n● We found robust expression of microglia-specific markers\nP2RY12, TMEM119 (now included in Extended Data Fig.\n1a).\n● We also found robust expression of infiltrating macrophage-\nspecific marker SIGLEC1 (aka CD169) in our putative\nmacrophage cluster that we previously identified with F13A1\nand have included in Extended Data Fig. 1a.\nExpression\n● CD44 and CD9 have low expression in both clusters. (Fig.\nR5)\nFigure R5. CD44 and CD9 gene expression in microglia and\nmacrophage.\nFigure R6. Table of de novo identification of macrophage markers.\nThe f_score, a harmonic mean of precision and recall, estimates the\nmarker gene’s importance. In this case, F13A1 was able to\ndelineate 82.96% of the macrophages from microglia.\nExcerpt from We identified two vascular subtypes (endothelial cells, 39.8% and\nthe manuscript mural cells 60.2%) (Fig 2e); and two immune subtypes (microglia,\n97.7% and macrophages, 2.3%) (Fig 2f). To confirm our immune\ncell subtypes, we examined expression of additional microglia-\nspecific markers P2RY12, TMEM119, and found high co-expression\nwith APBB1IP. Additionally, we examined the expression of\ninfiltrating macrophage-specific marker SIGLEC1 (aka CD169) and\nfound high co-expression with our F13A1 macrophage cluster\n(Extended Data Fig. 1a).\nRef 1.10 DEG Gender\nReviewer’s The authors report sex specific AUD changes and provide\ncomment supplementary tables for male AUDvsCON and female AUDvsCON;\nhowever, there was no discussion of how the analysis were\nperformed and specifics were not provided. Given the sample size\nthis analysis is likely underpowered but potential sex differences\nshould be examined and noted in the text.\nresponse of sex-specific DEG analyses and included them in new Extended\nData Fig. 4:\n1. Male AUD DEGs (Supplementary Data Table 8)\n2. Female AUD DEGs (Supplementary Data Table 9)\n3. Male vs Female Control differences (Supplementary Data\nTable 10)\nSpecifically, we intersected our original single-cell DEG list\n(Methods), with their respective sex exclusion/inclusion criteria. We\nhave expanded our results section (excerpt below) and also added\nnew Extended Data Fig. 4c-f:\n● Significant DEG counts for males Extended Data Fig. 4c, in\nboth directions for each major cell type (left) and INH sub cell\ntypes (right)\n● Significant DEG counts for males Extended Data Fig. 4d,\nwith their top enrichR Reactome (Pathway) terms.\n18\n● Significant DEG counts for females Extended Data Fig. 4e,\nin both directions for each major cell type (left) and INH sub\ncell types (right).\n● Significant DEG counts for females Extended Data Fig. 4f,\nwith their top enrichR Reactome (Pathway) terms.\nthe manuscript Data Table 8-10), consistent with previous observations in AUD.\nSpecifically, males have the most DEGs in the INH cell type,\nwhereas female DEGs are more evenly distributed across cell types\n(Extended Data Figure 4c, e). We also found pathway enrichment\ndifferences, specifically in WNT5A signaling in females and ERBB4\nsignaling in males (Extended Data Fig. 4d, f). In addition, male\npathways tended to be more cell type-specific, while females’\ntended to be enriched across cell types. 710 DEGs were common\nacross the male and female specific comparison or approximately\n55% or 25% overlap, respectively.\nReviewer’s Authors discuss the changes in chromatin accessibility that were\ncomment linked into differentially expressed genes, these are potentially\nimportant findings, and additional details should be provided\ndiscussing instances where the changes of chromatin accessibility\ndid not match observed gene expression changes and vice versa.\nAdditional discussion of these findings is warranted.\nAuthor’s We thank the reviewer for the suggestion and have warranted\nresponse additional discussion on this topic:\n1. DEGs that did not have an associated CRE accounted for\nless than 10%. In terms of genes that did not have an\nassociated CRE, only a few DEGs were discovered to not\nhave an associated peak: 4.03% (INH), 2.38% (EXC), 10.5%\n(OLI), 0% (OPC), 9.35% (END), 8.18% (AST), 4.40% (MIC).\n2. CREs that did not have an associated DEG can be\nexplained by 3D genomic structure. There are more CRE\npeaks that do not have associated DEGs. Based on the 3-\ndimensional structure of the genome, many of the CREs that\nare not linked to DEGs likely extend to distal parts of the\ngenome via chromatin-chromatin interactions. At present, we\ndo not have Hi-C data for the amygdala to directly determine\nthis. However, it stands to reason that around 30% of peaks\nare distal to DEGs, which makes long-range chromatin\ncontacts the most likely scenario.\n3. CREs with downregulated DEGs (63.95% of all CREs) –\nopen chromatin region doesn’t necessarily mean active\ntranscription (Fig. R7). Of the DEG-linked CREs, there are\n63.95% CREs associated with downregulated DEGs (Fig.\nR7). We posit several reasons for this:\na. the region may be accessible but not actively used,\ne.g., a poised enhancer or promoter that is open but\nlacks activating transcription factors;\nb. this could also indicate binding repressive\ntranscription factors including silencers or repressive\nenhancers;\nc. An AUD downregulated gene might still have\naccessible regulatory regions due to feedback loops\ntrying to restore expression.\nDEGs. Open chromatin region doesn’t necessarily mean active\ntranscription\nExcerpt from To identify disease gene regulatory mechanisms, we intersected\nthe manuscript our CRE linked genes (CLGs) with DEGs for each cell type. 36.05%\nof CREs were linked to an upregulated DEG.\nRef 1.12 Astro EDFig7\nReviewer’s Additional detail and discussion about the changes in cell types\ncomment other than inhibitory neurons should be provided and any\noverlaps between the identified differentially regulated elements in these\ndifferent cell types. This could be presented in the form of a\npotential supplemental table.\nAuthor’s We thank the reviewer for this excellent suggestion. Given the\nresponse importance of Astrocytes in AUD and the extensive number of\nDEGs and CREs we identified, we also highlighted expression\nchanges, P2G linkages, and cell-type-specificity for 2 genes present\nin astrocytes (Extended Data Fig. 7g):\n● NKAIN3 (logFC=-0.409 and FDR=2.06e-128)\n● FGFR3 (logFC=-1.298 and FDR=1.70e-96)\nExtended Data Figure 7: Astrocytes multi-omic profile. g, Violin\nplot showing AUD versus CON expression of FGFR3, NKAIN3 in\neach cell type with significant up (red asterisk) and down (blue\nasterisk) expression.\nmanuscript.\n● NKAIN3: UMAP plots of snRNA-seq gene expression on the\nExtended Data Fig. 7a, snRNA (left) and Extended Data\nFig. 7b, snATAC-seq (right) embeddings. We see increased\nexpression in astrocytes.\n● NKAIN3: Extended Data Fig. 7c, Chromatin accessibility\nsignal tracks highlighting NKAIN3 peak-to-gene links across\ncell types. Differences across the seven cell types are\nindicated by the dashed boxes. We see increased\naccessibility in astrocytes.\n● FGFR3: UMAP plots of snRNA-seq gene expression on the\nExtended Data Fig. 7d, snRNA (left) and Extended Data\nFig. 7e, snATAC-seq (right) embeddings. We see increased\nexpression in astrocytes.\n23\n● FGFR3: Extended Data Fig. 7f, Chromatin accessibility\nsignal tracks highlighting FGFR3 peak-to-gene links across\ncell types. Differences across the seven cell types are\nindicated by the dashed boxes. We see increased\naccessibility in astrocytes.\nExcerpt from We also identified peak-to-gene links (strong CREs, correlation>0.8)\nthe manuscript for genes previously implicated in addiction in other cell cell types\n(e.g. astrocytes). For example, we find an extensive number of\nATAC peak-to-gene links to NKAIN3 and FGFR3. NKAIN3 is a\nsodium/potassium pump gene that is highly expressed in astrocytes\n(Extended Data Fig. 7a-c) and changes in FGFR3 transcript have\nbeen linked to alcohol exposure of offspring of binge drinking\nbehavior48 (Extended Data Fig. 7d-f); both TSS’s are cell type-\nspecific to astrocytes (Extended Data Fig. 7c, f dashed box).\nAdditionally, both genes were significantly down regulated in AST\n(Extended Data Fig. 7g).\nReviewer’s The authors note a nominally significant up regulation of OPRM1 in\ncomment INH neurons. This finding is somewhat surprising since OPRM1 is\noften not dysregulated in transcriptome studies specifically\nexamining OUD. The donor metadata should be examined for\ncomorbid opioid use or reported OUD/MMD medication usage that\ncould potentially explain this interesting finding.\nAuthor’s We thank the reviewer for this feedback. We have removed this\nresponse sentence due to there only being nominal upregulation. However,\nfor the reviewers benefit, we have also verified that the nominal\nOPRM1 dysregulation is attributed to Alcohol Use Disorder and not\ncomorbid opioid use:\n● There are six AUD cases with reported opioid use (colored\nred; Fig. R8). We found that OPRM1 expression does not\ndiffer among these six cases compared with all other AUD\ncases (P-value=0.23, non-significant).\nFigure R8. OPRM1 Gene Expression by AUD-only Donor\nExcerpt from Interestingly, we found a nominally significant up regulation of the\nthe manuscript mu opioid receptor (OPRM1) in INH neurons cell type and in both\nthe DRD1 and LHX9 INH subtypes (Supplementary Data Table 2).\nReviewer’s Line 96: (mean reads/nuclei), how it is written it suggests that this\ncomment was the target sequencing depth.\nAuthor’s We have since changed the language to clarify the difference.\nresponse\nExcerpt from Here we present a single cell-type multi-omic analysis of 174,188\nthe manuscript high-quality nuclei (average 6028 reads/nuclei) from postmortem\nhuman brain tissue of AUD and neurotypical controls (N =50).\nRef 1.15 (minor) Correct\nReviewer’s Figure 4b is not referenced in the text\ncomment\nAuthor’s We thank the reviewer for pointing this out. Fig. 4b is now\nresponse referenced in the second paragraph of the Results section “Cis-\ngene regulation of calcium activity and glutamatergic signaling\ngenes is disrupted in AUD”.\nExcerpt from We then found 1,526,912 peak-to-gene links (correlation>0.45,\nthe manuscript FDR<1x10-4) using the union peaks and plotted side by side\nheatmaps of linked ATAC regions and gene expression (Fig. 4b).\nRef 1.16 (minor) Code\nReviewer’s Links to GitHub scripts were not working for this reviewer.\ncomment\nAuthor’s We thank the reviewer for the feedback. The GitHub link should be\nresponse working now: https://github.com/mjgirgenti/AUDsnCEA (Fig. R2).\nExcerpt from All code used in this study is freely available online and can be\nthe manuscript found at https://github.com/mjgirgenti/AUDsnCEA.\nRef 2.0 General Comments\nReviewer’s The manuscript by Lee et al provides insight into alcohol use\ncomment disorder in humans by investigating single cell chromatin and gene\nexpression data from post-mortem brain donors. The paper is well\nwritten and provides new data to the field. I’m not fully comfortable\nto provide an in-depth review of the bioinformatic analyses and\ntherefore focus on clinical aspects primarily.\nAuthor’s We thank the reviewer for their constructive comments and have\nresponse made a point-by-point response addressing them. Specifically, we\nhave:\n1. Clarified our writing in title, abstract, and introduction.\n2. Incorporated a more comprehensive list of clinical covariates.\n3. Added a functional validation using mouse amygdala in a\nmouse alcohol use model.\nRef 2.1 Title Clarify\nReviewer’s The title is slightly incorrect because the technique used is single\ncomment nucleus multiome.\nAuthor’s We thank the reviewer’s opportunity to clarify our sequencing\nresponse protocol. We have since changed the title.\nExcerpt from (Title) Central amygdala single-nucleus atlas reveals chromatin and\nthe manuscript gene transcription dynamics in human alcohol use disorder\nRef 2.2 Abstract Clarify\nReviewer’s I’d suggest including the sample size in the abstract.\ncomment\nAuthor’s We thank the reviewer’s suggestion and have since added the\nresponse sample size in our abstract.\nExcerpt from Here we present a multi-omic single nucleus study of ~175,000\nthe manuscript nuclei from 50 donors with alcohol use disorder (AUD) and non-\nchromatin accessibility in the human central amygdala.\nRef 2.3 Intro Clarify\nReviewer’s Line 89: A quick google search identifies snRNAseq studies in\ncomment human CeA, I may be wrong.\nAuthor’s We thank the reviewer for the opportunity to clarify. While there are\nresponse several existing snRNA-seq studies of the human amygdala,\n1. There are currently none specific to Alcohol Use Disorder or\nto subregion-specific Central Amygdala in Humans.\n2. However, we have toned down the language\nExcerpt from To date there have been no single-cell multi-omic studies of the\nthe manuscript human CeA in AUD.\nRef 2.4 DSM-IV combine\nReviewer’s The clinical spectrum includes individuals with alcohol abuse and\ncomment alcohol dependence, two very different phenotypes. What is the\nrationale and justification to combine both beyond the notion that\nDSM-V collapsed those two (can they be combined in this dataset\ngiven that the clinical heterogeneity within AUD is large)? The\ndonor’s age at death is young, what led to their demise? Do you\nhave information on the severity of AUD (mild, moderate, severe)\nand the timeline across life? Did donors consume alcohol at death\n(neurotox data?)? Do you have information on the\nneuropathological assessments of the donors (cases and controls),\nespecially for those >60yrs [neurodegen panel]?\nAuthor’s We thank the reviewer for the opportunity to clarify this point of our\nresponse study subjects.\n● First, the reviewer is correct. In the latest DSM, alcohol\nabuse and dependence have been combined into AUD.\nFrom the point of view of alcohol exposure, there is quite a\nhigh degree of overlap between the level of drinking\nassociated with alcohol abuse and alcohol dependence.\nGiven the importance of these studies and the scarcity of\nhuman postmortem tissue and the current diagnostic criteria,\n29\nwe felt that combining AUD cases with dependence and\nabuse was warranted. Each case is diagnosed by two\npsychiatrists and if there is disagreement in the diagnosis a\nthird is consulted.\n● Nevertheless, in an endeavor to alleviate the reviewers'\nconcerns we tested whether there were any biological or\nstatistical reasons to separate the cohorts. We calculated\nPearson correlations between the transcriptomic profiles of\nall pairwise sample comparisons (Fig. R9,10). While we do\nfind significant differences between AUD and CON samples,\nAlcohol Abuse and Dependence groups do not differ\nstatistically.\nFigure R9. Boxplot of Pearson correlations between each pair of\nsamples, and R10, Boxplot of Pearson correlations between each\npair of samples, after Log Normalization\n● We also projected all samples using PCA for all phenotypic\ndata (Fig. R11). We were able to discern clear distinctions\nbetween the AUDs (blue and green) and the CONs (orange).\nHowever the AUD subtypes (Dependence and Abuse)\nclustered together.\n● Taken together, we feel that it is both biologically and\nstatistically valid to combine Alcohol dependence and abuse\nsamples into an AUD sample as diagnosed by the brain bank\npsychiatrists.\n● We now include new meta data for:\n1. MOD (28 Natural, 15 Accidental, 6 Suicides, 1\nUndetermined) in Supplementary Data Table 1.\n2. We have since included AUD Severity (6 Mild, 9 Moderate,\nand 7 Severe) in Supplementary Data Table 1.\n3. We have since included the 5 AUD samples that had AUD\nNeurotox in addition to many other covariates in\nSupplementary Data Table 1.\nSupplementary Data Table 1. Samples Demographic Metadata\n● Finally, all donated brains underwent physical examination\nby a board-certified neuropathologist at both the macroscopic\nand microscopic level, please see our excerpt below on how\nthese examinations are conducted.\nthe manuscript neuropathologist at both the macroscopic and microscopic level. All\ntissue was screened for confounding neuropathologies, including\namyloid plaques, Lewy bodies (Parkinson’s disease and dementia\nwith Lewy bodies), and transactive response DNA binding protein-\n43 (TDP-43) (ALS, frontotemporal dementia, and limbic\npredominant age-related TDP-43 encephalopathy). All cases were\nalso evaluated at a gross (macroscopic) level to detect evidence of\natrophy, trauma, or infarction. Importantly, cases were excluded if\nmassive trauma occurred (from head injury) that severely damaged\nthe tissue, if the tissue was damaged from a stroke involving a large\nportion of the brain, if the decedent was on a respirator for an\nextended period of time, if the donor had brain cancer, or if the\ndonor had a history of HIV, AIDS, COVID-19, or other\ncommunicable disease.\nRef 2.5 Look Covariate\nReviewer’s Looking at the Extended Data Fig 1a, the details given regarding\ncomment the post-mortem samples is very limited. The U Pitt Brain Bank has\na large amount of data on their donors that is critical to include and\nto evaluate the data and this manuscript. For example: Co-\nmorbidities (anxiety disorders, depression, psychosis), somatic\ndisorders (hepatic function, malnutrition etc), timelines (AUD as\nlifetime diagnosis or at death, quantity of alcohol consumed,\nmanner of death, cause of death, medication taken, etc etc.) It is\ncritical to look at the clinical variables in depth as this is not a well-\ncontrolled animal study. In line 425 the authors refer to ED Fig 1a\nfor cause and manner of death, which is not there. Smoking is the\nonly comorbidity listed here and may be significantly different\nbetween groups (age, PMI as well?).\nAuthor’s We thank the reviewer for the opportunity to explain more about the\nresponse donors. We have now included additional details (Supplementary\nData Table 1) including:\n● Other Psychiatric Disorders\n● Lifetime Alcohol Consumption (Range 0~46)\n● AUD Severity (6 Mild, 9 Moderate, and 7 Severe)\n● AUD Neurotoxicity ATOD (for 5 AUD samples)\n32\n● MOD (28 Natural, 15 Accidental, 6 Suicides, 1\nUndetermined)\n● Medications ATOD\n● Illicit Drug Use (Range 0~4)\nin addition to our original covariates Sex, Race, Age, PMI, RIN.\nWe did find that while 80% of DEGs overlap (Extended Data Fig.\n4a) they are nevertheless significantly different datasets. The\nalcohol attributes (e.g. Lifetime Alcohol consumption), Smoking,\nMedications, Other Psych, and Illicit Drug Use were significantly\ndifferent between groups. Therefore, we have included a new DEG\nlist to reflect this (Supplementary Data Table 6).\nExtended Data Fig. 4a. Overlap of DEG analysis with new\ncovariates.\nSupplementary Data Table 1. Samples Demographic Metadata\nExcerpt from Sociodemographic and clinical details are listed in Supplementary\nthe manuscript Data Table 1 and include the manner of death, alcohol-related\nphenotypes, tobacco use, additional psychiatric disorders, and illicit\ndrug use.\nRef 2.6 Inflated DEGs\nReviewer’s The number of DEGs is high for the sample size and I think the\ncomment approach to use two different methods is good. However, can\nthe authors provide QQ plots or other data that may support the\nnotion that these results are not inflated? The inclusion of\ncovariates in the MAST approach is basic. I understand that\nincluding a myriad of clinical covariates is difficult given the\nsample size, however, I wonder if a surrogate variable analysis could\ncover those known and unknown factors. Although the study here is\nthe largest up to date on AUD, an independent replication would be\ncritical due to the limit in sample size.\nAuthor’s We agree with the reviewer that multiple models should be\nresponse employed in the detection of DEGs. Since our last submission, we\nhave provided the four additional models. Taken together, our\nresults show that our DEG analysis is likely not inflated:\n1. New DEGs which include additional clinical covariates and\npseudobulk analysis (Supplementary Data Table 6, 7). We\nhave found that, with the inclusion of new clinical variates\n(more conservative model) overlapped ~80% with previously\nfound DEGs (Extended Data Fig. 4a).\nExtended Data Fig. 4a. Overlap of DEG analysis with new\ncovariates.\n2. Batch-corrected DEGs via surrogate variable analysis. We\nhave performed surrogate variable analysis via “Harmony”.\nThe data now reflects batch-corrected DEGs to look for\nhidden batch effects (Fig. R12).\nFig. R12. UMAP Plots before and after surrogate variable analysis\n34\n3. Validation using small molecule FISH of GWAS DEGs.\nPlease see R1.4 Validation on page 3.\n4. Replication and validation using a mouse model of alcohol\nuse. We identified 1814 total DEGs across all major cell\ntypes with most falling in the INH (Extended Data Fig. 9d).\nWhen compared with the human AUD DEG dataset we found\na very high overlap of 678 DEGs (~40%) (Extended Data\nFig. 9e).\nExtended Data Figure 9: Mouse snMultiome data. d, Significant\nDEG counts in both directions for each major cell type. e, Intersection\nof alcohol-focused DEGs between human and mouse.\n5. QQ Plot to visualize the distribution of DEG p-values. We\nalso visualize our Wilcoxon-intersected MAST DEGs with a\nQQ Plot (Fig. R13). These results seem minorly inflated, but\nthis is a generally accepted consequence of single-cell\ndifferential expressed genes. Other methods which dealt with\nsingle-cell gene expression have also shown similar results.\nExcerpt from Ethanol Treatment (from Methods)\nthe manuscript Mice were delivered at seven weeks of age and allowed two weeks\nto acclimate to the vivarium, where they were housed in a reverse\nlight-dark cycle housing room. Two hours into the dark cycle, body\nweights were measured. Three hours into the dark cycle, mice\nreceived two sequential intraperitoneal injections each with a 10\nml/kg injection volume: the first was a solution of 5% DMSO in\nsaline directly followed by 2 g/kg ethanol dissolved in water (n = 2M,\n2F) or saline (n = 2M, 2F). Animals were randomly assigned to\ncontrol or ethanol treated groups. Three hours following these\ninjections, mice were euthanized by rapid decapitation and brains\nremoved and immediately snap frozen on dry ice and then stored at\n-80ºC. For tissue analyses, frozen brains were sliced into 2-mm\ncoronal sections on a metal brain matrix (Plastics One, Roanoke\nVA) and amygdala samples were collected bilaterally with a 1-mm\ntissue punch (Fine Science Tools, Foster City, CA). Bilateral\npunches were collected from a coronal section targeted at -1.43 mm\nfrom bregma, according to76. Nuclei were isolated from rodent\namygdala tissue as described in Brain Nuclei Isolation for Multiome\nmethods section.\nReviewer’s Given that the results presented are primarily associative, some\ncomment functional manipulation of genes or cell types in animal models\nwould be critically important to better understand the generalizability\nof the data.\nAuthor’s We thank the reviewer for this suggestion and agree that further\nresponse experimental validations will strengthen these findings. As\nsuggested, we have now added an animal model of alcohol use and\nperformed the same molecular analyses as the human tissue\nincluding(snMultiome including RNA and ATAC) after mice were\ntreated with an 2 g/kg injection of 100% EtOH.\n1. Mouse alcohol model validates CeA cell types: INH, EXC,\nOLI, OPC, END, AST, MIC on snRNA, snATAC, and\nsnMultiome embeddings, respectively (Extended Data Fig.\n9a-c).\nExtended Data Figure 9: Mouse snMultiome data. UMAP\nvisualization across seven cell types of a, snRNA-seq b, snATAC-\nseq across seven cell types and c, snMultiome\n2. Mouse alcohol model snMultiome-seq validates DEGs: We\nidentified 1814 total DEGs across all major cell types with\nmost falling in the INH (Extended Data Fig. 9d). When\ncompared with the human AUD DEG dataset we found\noverlap for 678 DEGs (~40%) (Extended Data Fig. 9e).\nDEG counts in both directions for each major cell type. e, Intersection\nof alcohol-focused DEGs between human and mouse. f, Top 25\nenrichR GO terms of the 1,814 unique AUD DEG set.\n3. Mouse alcohol model snMultiome-seq validates Chromatin\nPeaks:We performed peak-to-gene linkage analysis on our\nmouse snMultiome data and identified 7485 total in the\nmouse AUD model. Of these, 4710 overlapped with the\nhuman AUD CRE genes (62.9%) (Extended Data Fig. 9h).\nExtended Data Figure 9: Mouse snMultiome data. h, Intersection\nof alcohol-focused Cis-Regulatory Element–linked genes between\nhuman and mouse.\n4. Mouse alcohol model snMultiome-seq validates AUD GRN.\nWe found that combined KLF16, KLF 7 and KLF6 form a\nGRN that regulates 549 DEGs (Extended Data Fig. 9i). We\nfound 50 DEGs overlapped between the human and mouse\nGRN including DEG-GWAS gene CALN1 confirming its likely\nrole in alcohol abuse.\nnetwork showing the target genes for TFs KLF6, KLF7, and KLF16 in\nINH. Peak-to-gene correlation > 0.9 was employed in visualizing the\nnetwork.\nthe manuscript mouse model of alcohol abuse. We isolated frozen amygdala from\nmice treated with an 2 g/kg injection of 100% EtOH. We performed\nsingle nuclei isolation and snMultiome in the same manner as our\nhuman donor cohorts. We identified all major cell types across all\nmodalities: RNA (Extended Data Fig. 9a), ATAC (Extended Data\nFig. 9b) and combined Multiome (Extended Data Fig. 9c) for 19,072\nindividual nuclei. We identified 1814 total DEGs across all major cell\ntypes with most falling in the INH (Extended Data Fig. 9d). When\ncompared with the human AUD DEG dataset we found overlap for\nin neuronal cell types and this was reelected in gene set enrichment\nanalysis which identified terms related to neuron projection, synaptic\ntransmission, and glutamatergic signaling (Extended Data Fig. 9f).\nIn addition, we intersected our INH DEGs with a previous central\namygdala mouse model of alcohol withdrawal46 and found 38\noverlapping DEGs in INH and EXC neurons and MIC (Extended\nData Fig. 9g). We performed peak-to-gene linkage analysis on our\nmouse snMultiome data and identified 7485 total in the mouse AUD\nmodel. Of these, 4710 overlapped with the human AUD CRE genes\n(62.9%) (Extended Data Fig. 9h). To confirm the function of our KLF\n40\ngene regulatory network in human AUD (Fig. 5g), we constructed a\nGRN for these same TFs in our mouse model (Extended Data Fig.\n9i). We found that combined KLF16, KLF 7 and KLF6 form a GRN\nthat regulates 549 DEGs. We found 50 DEGs overlapped between\nthe human and mouse GRN including DEG-GWAS gene CALN1\nconfirming its likely role in alcohol abuse. Taken together, these data\nsupport the functional significance of KLF TFs in regulating calcium\nand glutamatergic signaling pathways identified in our dataset.\n40"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors provided a very thorough and thoughtful revision. My questions have been\naddressed adequately. The one remaining detail is that the new mouse data is short-\nterm forced administration and thus limits interpretation with human AUD. This caveat\nshould be addressed in the Discussion."}, "author_response": "We thank the reviewer for this suggestion. We have added the following text to th\nDiscussion:\nWe were also able to confirm the presence of a GRN of KLF6, 7 and 16 in our animal\nmodel of alcohol use (Extended Data Figure 9). However, this was a short-term administration\nof ethanol and thus limits some of the interpretation of this finding."}]

Decision letter