Downregulated transcription in chromosomal domains of midbrain dopamine neurons linked to schizophrenia

Domain: chemistry
Published: 2026-04-06
Editorial outcome: major_revision

Abstract
Ventral midbrain dopaminergic neurons are a key cell type for schizophrenia pathophysiology but information about cell type-specific genomic dysregulation in diseased brains is missing. We generated a unique midbrain functional genomics resource with 97 RNA-seq and 34 Hi-C chromosomal contact libraries for Nurr1+/NeuN+ dopaminergic and their surrounding Nurr1-/NeuN- nuclei, collected from donors diagnosed with schizophrenia (SCZ), bipolar disorder (BD) and compared to neurotypical controls. Among the 954 dopamine neuron genes specifically dysregulated in SCZ, 331 were downregulated, with selective enrichment for risk-associated synaptic plasticity and neuronal connectivity pathways and embedded within dopamine neuron-specific topologically associated chromosomal domains (TAD). Transcript-resolved analysis revealed 2,350 transcripts with altered expression in SCZ dopamine neurons, affecting key susceptibility genes such as the FOXP1, MAPK10, PCM1 and NRNX1. Therefore, genomic dysregulation in the ventral midbrain of subjects diagnosed with SCZ selectively affects dopaminergic neurons, and includes a unilateral association of genetic risk with down-, but not upregulated transcription at the sites of highly organized chromosomal domains harboring neuron-specific genes with complex transcriptional architectures.

Introduction
Schizophrenia (SCZ) and bipolar disorder (BD) are major psychiatric disorders each affecting approximately 1% of the population worldwide. No unifying neuropathology exists, but there is a complex genetic risk architecture that includes, for SCZ, non-coding variants within 287 distinct loci each contributing a very small but tractable liability 1, together with a few dozen very rare mutations and variants of much higher disease penetrance 2. Therefore, exploration of genome organization and function in the brains of subjects diagnosed with SCZ and BD, and comparison with neurotypical brains, carries significant potential to uncover disease-sensitive transcriptomic and epigenomic alterations. To date, however, at base pair resolution, very few genome-scale case-control brain reference sets for gene expression, DNA methylation, chromatin accessibility, and histone marks exist and moreover, these are exclusively limited to the frontal lobe including its higher association cortex 3–7 and portions of the rostral striatum 8,9, and the medial temporal lobe including the hippocampus 10,11.

Unfortunately, to date, there are no genomic studies on critical nodes of the psychosis circuitry outside of the aforementioned frontal and temporal lobe areas, despite the numerous times they have been implicated. For example, functional neuroimaging studies identified the dysregulation of dopaminergic neurons in the ventral midbrain as the critical presynaptic component responsible for altered dopaminergic innervation and signaling in SCZ and BD forebrain 12. Moreover, a pilot study on brain tissue homogenate from 28 diseased brains, measuring expression levels of 12 markers commonly associated with presynaptic dopaminergic signaling, reported alterations in a small subset of subjects 13. However, perhaps the most compelling rationale for studying dopaminergic systems in schizophrenia and related diseases is provided by the undisputed fact that to date, antipsychotic drugs (e.g. haloperidol , aripiprazole, cariprazine) targeting D and other dopamine receptor systems are still the mainstay of treatment worldwide 14.

This lack of genomic reference sets from SCZ and BD midbrain could be a reflection of technical challenges presented by the very miniscule share of dopaminergic neurons in brain tissue. Indeed, even in some of the most compact anatomical structures harboring dopaminergic neurons such as the substantia nigra (SN) and bordering ventral tegmental area (VTA), also known as dopaminergic cell clusters A9 and A10 respectively, less than 2% of the total cell population express presynaptic dopaminergic marker genes 15–17 and thereby are likely to evade functional genomic profiling if bulk tissue without cell-targeted enrichment would serve as the input material.

Here, we bypass these limitations by applying immunotagging combined with fluorescence-activated nuclei sorting (FANS) of ventral midbrain tissue to separately collect dopaminergic neuron nuclei and their surrounding non-neuronal nuclei, followed by deep cell-type specific profiling of the nuclear transcriptome (nucRNA-seq) together with ‘HiC’ 3D genome chromosomal conformation mapping to differentiate between active and inactive chromosomal compartments and domains. Our functional genomics resource encompasses 97 subject and cell-type-specific nuclear transcriptomes and 34 subject and cell-type-specific 3D genomes. Each dataset is deeply sequenced, with average read counts of 59 and 649 million reads per RNA-seq and Hi-C library respectively. We thereby provide reference maps for the ventral midbrain of a SCZ and BD cohort compared to neurotypical controls and report disease-specific alterations in dopamine neuron-specific gene expression in the context of the underlying genetic risk architecture and 3D genome organization.

Results
Dopaminergic neurons in the midbrain are enriched with SCZ risk genes .To understand the role of dopaminergic neurons in SCZ in the midbrain region, we utilized a postmortem brain cohort comprised of 73 donors (1) from the CommonMind Consortium collection18. As illustrated in Fig. 1A, we first profiled a total of 97 RNA-seq libraries by immunotagging nuclei extracts from ventral midbrain SN (area A9) for previously validated 1 FANS sorting of dually labeled nuclei that were NeuN+ (pan-neuronal marker) and Nurr1+ (nuclear receptor subfamily 4 group A member 2 (NR4A2), a transcription factor crucial for maintenance and development of dopaminergic neuron) 19,20 (Fig. S1A). From the same tissues, we sorted and collected the surrounding non-neuronal (predominantly glial) Nurr1-/NeuN- nuclei (AFig. S1B). Altogether, we generated 41 Nurr1+/NeuN+ dopamine neuron-specific (‘dopa neurons’) RNA-seq libraries, comprising 24 control, 8 SCZ, and 9 BD ventral midbrain samples. Additionally, 56 Nurr1-/NeuN- non-neuron specific (‘non-neurons’) were collected from 26 controls, 22 SCZ and 8 BD midbrains. The samples were sequenced to an average depth of 59 ± 5.3 million reads (Fig. S2).

All processed data went through a series of rigorous quality control steps as shown in Fig. S2 and explained in detail in Methods (Supplementary Data 1). Metadata for all QC-passed samples, including sex, ancestry, age, postmortem interval, and RIN, stratified by cell type and diagnosis are shown in Fig. S3C, illustrating demographic and quality metrics across groups. To assess overall transcriptomic structure, we performed principal component analysis (PCA) on log CPM-normalized expression, which clearly separated dopaminergic and non-neuronal samples (Fig. 1B). We next validated the cell-type specificity of the sorted populations through deconvolution 21 against a published single-cell ventral midbrain reference22. This analysis confirmed that FANS-sorted dopaminergic nuclei were highly enriched for dopaminergic signatures (~80.9% ± 0.0206), while non-neuronal nuclei were enriched for non-neuronal signatures (~97.9% ± 0.00320) (Fig. 1C and S2E).

Our initial analysis focused on exploring the statistical differences in the transcriptomes of Nurr1+NeuN+ dopa neurons and Nurr1-NeuN- non-neurons in controls. As expected, the majority (~25%) of variation in transcriptome expression came from cell types (dopa and non-neurons) (Fig S2B). Furthermore more than two-thirds of the transcriptome (20,549 out of 30,687, 67%) showed significant changes at FDR < 0.05 including N=10,383 genes (33.84%) in dopa neurons and N=10,166 (33.13%) in non-neurons (Fig. 1D and Supplementary Data 2A). Inspection of known cell type marker genes further confirmed the specificity of our sorted populations. For example, dopaminergic neuron samples exhibited strong expression of Tyrosine Hydroxylase (TH) (logFC = 2.8;FDR = 7.7 × 10⁻⁷), the rate-limiting enzyme in dopamine synthesis, while non-neuronal samples showed robust expression of Olig2 (logFC = –5.04; FDR = 3.8 × 10-19), a transcription factor specific to oligodendrocytes, which constitute the majority (>65%) of cells in the SN 15 (Fig. 1E).

Notably, dopa-specific transcriptomes showed strong enrichments for psychiatric and neurological disease risk genes including SCZ, BD, MDD (major depressive disorder), and PD (Parkinson’s disease) using MAGMA 23 (Method and Supplementary Data 2B). In contrast, the non-neuron-specific transcriptome showed enrichment for various immunological and other non-brain related traits such as height (Fig. 1F). Functional pathway analysis of the dopa-specific transcriptome was predominantly associated with neuron-specific pathways like synaptic signaling, sharply contrasting with the non-neuron specific transcriptome which was primarily enriched for immune response and cytokine production pathways (Fig. S2D and Supplementary Data 2C). Together, these results establish that our sorted populations capture highly cell type–specific transcriptomic programs, providing a robust framework to investigate the role of dopaminergic neurons in SCZ.

Downregulated dopa neuronal transcriptome in SCZ is associated with brain-related traits. We next assessed disease-related transcriptional changes in 8 SCZ and 9 BD donors compared to 24 controls for dopa neuronal samples, and in 22 SCZ and 8 BD donors compared to 26 controls for non-neuronal samples. To account for biological and technical variability, we modeled variance contributions in gene expression across covariates (Fig. 2A, Fig. S3A) and included cell-type proportion estimates i.e. dopa neuron proportion in dopa neuronal samples and non-neuronal proportion (oligodendrocytes) in non-neuronal samples (Fig. 1C) as covariates. This approach ensured that disease-associated changes reflect true gene dysregulation rather than differences in cellular composition across donors.

Using this model, we identified 954 differentially expressed genes (623 up, 331 down; FDR < 0.05) in SCZ dopa neurons (Fig. 2B, Supplementary Data 3A). No FDR-significant changes were observed in the non-neuronal samples despite the higher statistical power to detect differences in SCZ (Fig. S5A, Supplementary Data 3B). Functional enrichment revealed that downregulated genes in SCZ dopaminergic neurons were enriched for pathways related to neuronal connectivity and signaling, including regulation of neurotransmitter levels, neuronal system, ensheathment of neurons, and positive regulation of nervous system development (Fig. 2C). By contrast, upregulated genes mapped largely to structural and organizational pathways including one neuro related pathway, synaptic signalling. (Fig. S3B, Supplementary Data 3C).

MAGMA analysis confirmed that downregulated genes, but not upregulated genes, were significantly enriched for SCZ genetic risk (Fig. 2D, Supplementary Data 3D). Among the top disease-associated downregulated genes were SLC6A3 (dopamine transporter, DAT) and SLC18A2 (vesicular monoamine transporter, VMAT2), key regulators of dopamine reuptake and storage (Fig. S4A). Their repression is consistent with long-standing evidence of dopaminergic dysfunction in SCZ, including the dopamine hypothesis and prior postmortem findings 24,25. In addition, top upregulated genes in dopaminergic neurons from SCZ donors include MTND1P23, CTD-2336H13.1, and RP11-149A16.12 (Fig. S4A). While their precise roles in neuronal function are less well characterized compared to DAT and VMAT2, their consistent upregulation highlights disease-associated transcriptional shifts that may reflect compensatory processes, mitochondrial activity, or structural reorganization in SCZ neurons. Together, these findings emphasize that SCZ pathophysiology involves both suppression of key dopaminergic regulators and activation of other molecular programs, reinforcing the complexity of transcriptional dysregulation in the disorder.

To evaluate potential medication effects, we compared antipsychotic-exposed and unexposed donors within SCZ (1 unexposed, 8 exposed) and BD groups (7 unexposed, 2 exposed), separately for neuronal and non-neuronal samples. Across all comparisons, no significant genes were detected at FDR < 0.05 (Fig. S4B-E, Tables S3E-H).These findings are consistent with prior reports suggesting that antipsychotic treatment may have only a limited impact on postmortem transcriptional profiles 26,27; however, this cannot be stated with certainty given the limited sample size.

For BD donors, no differentially expressed genes were detected relative to controls in either neuronal or non-neuronal fractions (FDR < 0.05; Fig. S5A, Tables S3I-J). Nevertheless, comparison of disease signatures revealed a negative concordance between SCZ and BD in dopaminergic neurons (ρ = -0.34, p < 2.2 × 10⁻¹⁶) and a positive concordance in non-neurons (ρ = 0.087, p < 2.2 × 10⁻¹⁶) (Fig. S5B-C).

To assess subtype specificity, we leveraged single-cell reference data from Wang et al. 22 and tested cluster-level enrichment of SCZ-associated genes. Downregulated genes were most enriched in the C7 dopaminergic subtype (TH/SLC6A3/ALDH1A1-positive), followed by the C6 dopaminergic subtype (TH/SLC6A3-positive) and the C9 GABAergic subtype (Fig. S7A), indicating a selective repression of dopaminergic programs in SCZ. Of note, modest enrichment was also observed in the non-neuronal C0 cluster, although this signal was weaker and less specific than that in dopaminergic subtypes.

Finally, we compared midbrain findings with published SCZ-associated genes from homogenate dorsolateral prefrontal cortex (DLPFC) and anterior cingulate cortex (ACC) in the CommonMind Consortium dataset 18. In contrast to midbrain dopa neurons, DLPFC and ACC homogenates showed enrichment of upregulated but not downregulated genes for SCZ risk (Fig. S6A-C, Supplementary Data 3K). Correspondingly, SCZ vs control and BD vs control gene expression changes in DLPFC and ACC were positively correlated (DLPFC ρ = 0.57; ACC ρ = 0.52; both p < 2.2 × 10⁻¹⁶), in sharp contrast to the negative concordance observed in midbrain dopa neurons (ρ = -0.34, p < 2.2 × 10⁻¹⁶; Fig. S6D-F). Consistent with this regional divergence, comparison of schizophrenia-associated differential expression between midbrain and DLPFC revealed a modest but significant negative correlation (ρ = –0.21, p = 3.1e–07), indicating regionally opposing transcriptional changes (Fig. 2E). Together, these results reveal a region- and cell-type-specific pattern of dysregulation: midbrain dopa neurons in SCZ exhibit broad gene expression repression of neuronal genes tied to genetic risk, whereas cortical regions show the opposite trend, with risk associated primarily with upregulated expression.

SCZ-associated genes exhibit concordant or discordant transcript-level effect sizes. Alternative splicing and variation in transcription start/stop sites have been implicated in SCZ pathogenesis 26,28. Because gene-level analyses may mask such isoform-specific effects, we examined transcript-level regulation in midbrain sorted populations: dopa neurons and non-neurons. We first tested 90,479 transcripts from 38,656 genes using limma but, given the burden of multiple testing corrections, this yielded only 8 significant transcripts in dopa neurons from SCZ versus controls (FDR < 0.05), while no significant differences were observed in BD versus controls or in non-neuronal populations. (Method, Supplementary Data 4A-D).

Given this result, we implemented the random multivariate meta-analysis method, remaCorr 29 30, to examine the relationship between disease-associated changes in Tx expression. remaCorr integrates the estimated SCZ-specific effect sizes of each transcript isoform as a fixed effect, incorporating the inter-correlation matrix between these effect sizes 31. We further categorized remaCorr-significant genes into concordant and discordant based on whether their Tx showed uniform or opposing effect sizes (Methods). This framework clarifies Tx regulatory interpretation, as discordant genes often reflect more complex isoform-specific mechanisms that may be obscured in gene-level analyses 32–34 (Fig. 3A).

In SCZ dopa neurons, remaCorr identified 294 significant genes comprising 2,350 Tx (FDR < 0.05) (Fig. S8B, Supplementary Data 5A), including 164 concordant (73 concordant downregulated and 91 concordant upregulated) and 130 discordant cases (Fig. 3B). Only ~34% overlapped with limma-based results, reflecting that standard gene-level analysis captures concordant genes well but largely misses discordant ones, whose opposing isoform changes cancel at the aggregate level. Consistent with this, concordant genes showed strong π1 (pi-one) enrichment (~0.9), whereas discordant genes displayed weaker enrichment (~0.65). π1 is a statistical estimate of the proportion of true positives among the tested associations, with higher values indicating stronger and more reproducible signals. Discordant genes were also structurally more complex, with a higher number of transcripts (median: 18 ± 9.5 vs. 5 ± 2.9 transcripts) and longer gene bodies (median: 113,352 ± 274,736 vs. 28,750 ± 102,257 bp) compared to concordant genes (Fig. S8C-D). Also, discordant genes exhibited higher variation in disease-associated changes, as reflected by greater standard deviation of log₂ fold-change (Fig. S9A).

Although functional enrichment analyses did not identify neuro-related pathways (Fig. S9C, Supplementary Data 5B), both concordant downregulated and discordant sets were significantly associated with SCZ risk in MAGMA (Fig. 3D, Supplementary Data 5C). For example, WWOX, a neurodevelopmental gene linked to epilepsy and neurodegeneration 35, showed discordant transcript regulation (Fig. 3E). Finally, LDESC revealed that discordant downregulated genes were most strongly enriched for SCZ heritability (Fig. S9B).

In BD dopa neurons, remaCorr identified 186 significant genes comprising 1715 Tx (FDR < 0.05) (Fig. S10A) including 89 concordant (51 concordant up and 38 concordant downregulated) and 97 discordant cases (Supplementary Data 5D). Discordant genes, again, were structurally more complex, with a higher number of transcripts (11 ± 9.13 vs. 3 ± 4.42 transcripts) and longer gene bodies (median: 117,819 ± 401,078 bp vs. ~27,700 ± ~195,000 bp) compared to concordant genes. Furthermore, both concordant and discordant Tx genes did not show enrichment for BD risk variants for FDR < 0.05 (Fig. S10B, Supplementary Data 5E). Lastly, in the non-neuronal fraction neither SCZ nor BD comparisons yielded significant genes at FDR < 0.05 (Fig. S10C-D, Supplementary Data 5F-G).

Together, these results demonstrate that isoform-level dysregulation, especially discordant downregulation, represents an important axis of genetic risk in SCZ, revealing regulatory complexity not captured by gene-level analysis.

3D genome structures associated with altered gene expression in SCZ dopa neurons. Non-randomly occurring ‘3D’ chromosomal conformations play an important role in genome organization and function. Examples include the kilo-to megabase spanning ‘A’ and ‘B’ chromosomal compartments associated with chromatin more permissive (‘A’) or restrictive (‘B’) transcriptional environment, in addition to ‘self-folded’ topologically associated domains (TADs) and chromosomal loop formations. Previous work conducted in DLPFC of SCZ and BD samples has resolved some of these 3D genome elements at chromosomal sites subject to epigenomic and transcriptional dysregulation 6. To investigate whether SCZ-associated gene dysregulation in midbrain dopa neurons is influenced by 3D genome architecture, we generated 34 deeply sequenced Hi-C libraries from sorted nuclei (17 dopa, 17 non-neuronal; Fig. 1A, Supplementary Data 6). Each set yielded >1.5 billion cis contacts, enabling high-resolution 3D genome maps of the ventral midbrain (Fig. S11C). Clustering of Hi-C profiles using HiCRep 36 separated dopa and non-neuronal samples (Fig. S11D-E), confirming cell type specific chromosomal conformations consistent with their distinct transcriptomes (Fig. 1D).

We next defined topologically associated domains (TADs) per fraction (N= 5,966 dopa and N=6,348 non-neuronal TADs) using TOPDOM 37. TAD boundaries in dopaminergic neurons (Fig. S12 A-D) were strongly enriched for CTCF 6 (3.283 fold enrichment, P = 1.19e–34) and the active promoter mark H3K4me3 (2.56 fold enrichment, P < 1e–4), indicating structural insulation at sites of transcriptional activity. Similarly, boundaries showed high enrichment for the enhancer-associated histone mark H3K27ac (2.91 fold enrichment, P = 1.19e–34), while exhibiting marked depletion of the repressive histone mark H3K27me3 (1.964 fold enrichment, P = 1.14e–34) 38. Consistent patterns were observed in non-neuronal cells (Fig. S12E-H), where TAD boundaries were enriched for CTCF (2.772 fold enrichment, P = 1.19e–34) and H3K4me3 (2.765 enrichment, P = 1.19e–34), as well as H3K27ac (1524 fold enrichment, P < 1e–4). Like in dopaminergic neurons, boundaries were significantly depleted for H3K27me3 (2.485 fold enrichment, P = 1.15e–34) 38. These findings indicate that TAD boundaries in both cell types are flanked by transcriptionally active and structurally defined chromatin, consistent with their role as organizing units of the 3D genome. The high concordance of boundary-associated chromatin marks across dopaminergic and non-neuronal cells supports a conserved architectural framework, despite broader transcriptomic and chromatin state differences between these populations. Similarly, eigenvector analysis identified a total of 39,023 compartments in dopa and 38,991 compartments in non-neurons, using Cscoretools39 (Supplementary Data 7A), revealing extensive cell type-specific switching between active (A) and inactive (B) chromatin states (Fig. S14B) indicating that risk variants are concentrated within transcriptionally active, regulatory-rich chromatin domains.

Specifically, 10,844 bins (27.8%) transitioned from A to B (active in dopa, inactive in non-neurons) and 6,545 bins (16.8%) transitioned from B to A (inactive in dopa, active in non-neurons), while the remainder remained stable as either active or inactive in both cell types (A→A or B→B) (Supplementary Data 7B). Furthermore, genes located in stably active A compartments shared between both cell types (A→A) showed a significant positive association with SCZ (Fig. S14C). Genes that gained activity in non-neuronal cells (B→A) exhibited nominal enrichment (P < 0.05) for SCZ, suggesting that these regions may represent dynamic chromatin domains that become active during neurodevelopmental transitions. In contrast, genes located in consistently inactive (B→B) or deactivated (A→B) compartments showed no enrichment, suggesting that SCZ risk is primarily associated with loci in or transitioning toward active chromatin environments.

Fig. 4A shows a representative dopaminergic Hi-C contact matrix at a chromosome 1 locus spanning 12 genes, organized into A/B compartments, TADs, and loops at 50 kb resolution. We then examined the spatial organization of SCZ-sensitive genes (from limma and remaCorr; Tables S3A, S5A) within the 3D genome to assess whether they preferentially reside within TADs. We then asked whether SCZ-associated genes identified from gene- and transcript-level analyses (Tables S3A, S5A) preferentially localize to specific 3D features. Downregulated genes within TADs showed the strongest enrichment for SCZ GWAS risk compared to upregulated genes or downregulated genes outside TADs (Fig. 4 B, Tables S7B, S7D). For genes with concordant Tx, only genes with downregulated Tx outside TADs were enriched (Fig. 4B, Tables S7B, S7E). By contrast, genes with discordant Tx showed robust enrichment both inside and outside TADs, with the strongest effect inside TADs (Fig. 4B, Supplementary Data 7C).

To better understand the structural complexity of dysregulated genes, we next assessed the extent to which dysregulated SCZ dopa neuron genes spanned multiple TADs. A gene was considered to extend across multiple TADs if its full annotated genomic span (from the outermost TSS to the TES across all isoforms) overlapped more than one TAD, while the TSS remained anchored within a single domain. A significantly larger fraction of genes with discordant Tx (34/130) crossed multiple TADs compared to genes with concordant upregulated (10/73) or downregulated (9/91) Tx (Fig. S14A; Supplementary Data 7E,F). These discordant, multi-TAD genes included large, highly spliced, SCZ-implicated loci such as NRXN140 41, FOXP142, MAPK1011, and PCM143 (Fig. S14A, Supplementary Data 7D). This is in line with prior studies on long neuronal genes with complex isoform structures that span multiple regulatory domains 44,45

Together, these results demonstrate that SCZ-associated dysregulation in midbrain dopaminergic neurons is shaped by 3D genome architecture (Supplementary Data 7G). In particular, downregulation within TADs and isoform-level discordance in multi-TAD genes represent key axes linking chromatin topology to SCZ genetic risk.

Methods
All brain specimens were obtained with informed consent and/or through brain donation programs from the NIMH Human Brain Collection Core (HBCC). NIMH HBCC specimens were collected under NIH Institutional Review Board protocols 90-M-0142 and 17-M-0173 until 2018, after which HBCC research was determined to be nonhuman subjects research, with oversight by the HBCC Oversight Committee. All procedures followed institutional guidelines from Institutional Review Boards of the brain banks for the use of postmortem human brain tissue.

RNA-Seq Library Preparation and Nuclei Isolation: An aliquot of tissue was dissected from each brain sample from the substantia nigra (A10) area of the ventral midbrain, and dounced around 50x in 1ml of lysis buffer (0.32M sucrose, 5mM CaCl , 3mM Mg(Ace) , 0.1mM EDTA, 10mM Tris pH8, 0.5mM DTT, 0.1% Triton X-100) with 400U Recombinant RNase Inhibitor (Cat. 2313; Takara Bio USA, San Jose, CA). An additional 4ml of lysis buffer (without RNase inhibitor) was then added, and the sample solution was dounced an additional 20x until completely homogenous.

The sample homogenate in lysis buffer was transferred to an ultracentrifuge tube (Polypropylene Centrifuge Tubes ⅝ x 3 ¾ in., Ref. 361707; Beckman Coulter, Indianapolis, IN) and underlayed with 9ml of sucrose buffer (1.8M sucrose, 3mM Mg(Ace) , 0.5mM DTT, 10mM Tris-HCl pH8), which was then ultracentrifuged in a SureSpin 630 (17mL) Rotor (Cat. 79363; Thermo Fisher Scientific, Asheville, NC) with 24,000rpm (106,803x g) for 1hr at 4ºC. After removing the supernatant post-centrifugation, the nuclei pellet was resuspended in 1 ml of 0.1% BSA and transferred to a sterile 1.5 ml tube. 2 µl of Anti-Nurr1 antibody (Anti-Nurr1 antibody produced in rabbit, Cat. N4663 or Anti-NR4A2 (N-terminal) antibody produced in rabbit, Cat. N6413, lot no. 118K0446V; Sigma-Aldrich, St. Louis, MO) conjugated in a 2:1 ratio to secondary antibody Alexa Fluor 647 (Goat anti-Rabbit IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor™ 647, Cat. A21245, lot no. 2098544; Invitrogen, Waltham, MA) and 1.5 µl of NeuN neuronal antibody pre-conjugated to Alexa Fluor 488 (IAnti-NeuN Antibody, clone A60, Alexa Fluor™ 488 conjugated, Cat. MAB377X, lot nos. 3481408, 3193849, and 3072355; MilliporeSigma, Burlington, MA) were then added to the sample, and the sample was incubated on a rotator at 4ºC for 2 hrs while protected from light. After 2 hrs, 1 µl of nucleophilic dye DAPI (4',6-Diamidino-2-Phenylindole, Dihydrochloride, Cat. D1306; Invitrogen, Molecular Probes, Eugene, OR) was added to the sample. A BD FACSAria Cell Sorter (BD Life Sciences, San Jose, CA) was used to sort approximately 50,000 nuclei each of NeuN+/Nurr1+ and NeuN-/Nurr1- nuclei populations into collection tubes pre-filled with Trizol LS Reagent (Cat. 10296028; Invitrogen, Life Technologies, Carlsbad, CA). In our hands, 50,000 nuclei as input material was optimal for generating high-quality RNA-seq libraries with the method and procedures described here. The remaining types of nuclei were not sorted and discarded. After sorting, sorted sample volume was measured and compared to the volume of Trizol LS Reagent pre-filled in the tube, and if necessary additional Trizol LS was added to bring the sample to a 3:1 volume ratio of Trizol LS to sample. Sorted samples in Trizol LS were then vortexed, flash-frozen, and stored at −80ºC (Fig. S1A,B).

NeuN/Nurr1 double positive and double negative samples were removed from -80ºC and thawed on ice in an RNase-free environment. An additional 20% volume of Trizol LS was added to each sample, then samples were vortexed for 30 seconds and incubated at room temperature for 10 minutes. As per the Trizol LS Reagent User Guide (Pub. No. MAN0000806; Invitrogen, Life Technologies, Carlsbad, CA), chloroform was added to each sample at 0.2ml chloroform per 0.75ml of total Trizol LS, and the samples were vortexed and incubated at room temperature for 3min. Samples were then centrifuged for 15 minutes at 12,000x g at 4ºC, after which the aqueous phase containing RNA was transferred to a sterile 1.5ml tube. Next, the RNA samples were purified using the Zymo RNA Clean & ConcentratorTM -5 kit (Cat. R1014; Zymo Research, Irvine, CA), including an in-column DNase I treatment which was incubated for 30min instead of the recommended 15min. During the final step of that kit, RNA was eluted into 8.5µl of nuclease-free water. Lastly, the ArcticZymes Heat & Run gDNA Removal Kit (Cat. 80200; ArcticZymes Technologies, Tromsø, Norway) was used as a secondary DNase treatment to remove genomic DNA from each sample.

Library Preparation and Sequencing: Before performing the ArcticZymes DNase treatment, the RNA concentration and RNA integrity number (RIN) were obtained for each sample using the Agilent RNA 6000 Pico kit (Cat. 5067-1513; Agilent, Santa Clara, CA) on the Agilent 2100 Bioanalyzer system. Those values were used to determine the sample input into the SMARTer® Stranded Total RNA-Seq Kit v2 - Pico Input Mammalian (Cat. 634413; Takara Bio USA, San Jose, CA) which was used to generate cDNA libraries for sequencing. The complete libraries were sequenced by the New York Genome Center on the Illumina NovaSeq-6000 platform, with 100bp paired-end reads at a depth of approximately 60 million read pairs per sample. We have assembled a dataset comprising 111 RNA-seq libraries, originating from FANS-sorted Nurr1+/NeuN+ immuno-tagged midbrain dopaminergic neurons, as well as their neighboring glial and other non-neuronal nuclei (Nurr1-/NeuN-) (Fig. S1A,B). These 111 RNA-seq libraries are categorized into 53 samples of dopaminergic neurons and 58 samples of non-neuronal cells.

Within the dopaminergic neuron subset, 28 are designated as control samples, 14 are from subjects diagnosed with schizophrenia, and 11 from subjects diagnosed with BD. Among the non-neuronal samples, 28 represent controls, 22 are from subjects diagnosed with SCZ, and 8 are from subjects diagnosed with BD (Fig. 1A, Supplementary Data1).

RNA-seq data processing, alignment and quality control: For the analysis of RNA sequencing data, the initial steps involved the alignment and pre-processing of raw reads. Trimmomatic (v0.36) 49 was employed to trim the raw reads, followed by the mapping of the processed reads to the human reference genome GRCh38.v24 50 using STAR (v2.5.3a) 51 (Fig. S2A). The resulting BAM files contained paired-end reads that were successfully mapped, inclusive of those spanning splice junctions. After alignment, expression quantification was performed and summarized at both the gene and transcript levels using GENCODE v27 annotations 52,53. Quality control measures, including checks for sample mislabeling, genotype discrepancies, and sex mislabeling, were conducted on 97 RNAseq libraries using the mbv function from QTLtools (v1.3) 54 (Fig. S2B-E). A scatterplot of log₂(CPM) expression for the XIST (chrX) and UTY (chrY) genes confirmed accurate annotation of reported biological sex for all donors (Fig. S2F).

Deconvolution of dopa neuron and non-neuronal samples: To estimate dopaminergic and non-neuronal cell-type proportions per samples for both dopa and non-neuronal fractions, we applied the MuSiC21 tool using a reference single-cell RNA-seq dataset from Wang et.al 22 of the human midbrain. Top 50 DEGs per annotated cluster were used as reference cell type specific markers. The reference comprised of all 11 clusters, annotated as microglia (1), astrocytes (2), oligodendrocytes (0, 3), endothelial cells (4), oligodendrocyte progenitor cells (5), neurons (6, 7, 9), fibroblast-like cells (10), pericytes (8), and T cells (11). For downstream modeling, we grouped clusters into three categories: dopa neurons (clusters 6, 7), non-dopa neurons (cluster 9), and non-neuronal populations (all other clusters) (see Figs. 1C and S2E) . Deconvolution was performed using the music_prop function with cluster identities and sample metadata derived from the reference Seurat object.

Data normalization, selection of covariates, and statistical analysis of differences in gene expression of cases and control: Raw gene expression counts were normalized using the voomWithDreamWeights 55. Filtering was applied to retain only expressed genes, defined as those with counts-per-million (CPM) values greater than 5 in at least 10% of the samples. This step was done separately on three datasets: 1) Dopa and non-neuronal 2) Dopa neuronal, and 3) non-neuronal samples

Covariates were identified across the above mentioned three datasets through forward stepwise model selection using the mvForwardStepwise function from the mvIC R package (V1.6.3) 56, employing the Bayesian Information Criterion (BIC) for model comparison. Specifically, each additional covariate was tested to ensure the BIC difference between models was at least 5. After this, our final models for all comparisons are following:

1) Dopa vs non-neurons in controls:
Dopa and Non-neurons expression: ~ celltype + RIN + (1|Sex) + Age + Age2 + PMI + GC content

2) SCZ vs controls and BD vs controls comparisons in dopa neuronal samples:
Dopa neurons expression: ~ Diagnosis + RIN + (1 | Sex) + Age + Age2 + PMI + (1 | adapters) + GC content + dopa proportion (cluster 7)

3) SCZ vs controls and BD vs controls comparisons in non neuronal samples:
Non-neurons expression: ~ Diagnosis + RIN + (1 | Sex) + Age + Age2 + PMI+ GC content + Non-neuron proportion (cluster 0)

Briefly, celltype in (1) is a categorical variable defined as dopa or non-neuronal, RIN is RNA Integrity Number which is a measure of RNA quality and Dx is the Diagnosis (SCZ, BD and Controls). Fig. 2A, Fig. S2.BA, and Fig.S3A illustrate the variance explained by identified covariates in gene expression for three datasets utilizing the variancePartition 55,57. Statistical analysis to measure the significant gene expression changes across all comparisons were performed using the dream function from the variancePartition 55,57.

To assess potential confounding by antipsychotic exposure (e.g. Haldol, chlorpromazine), we performed a secondary differential expression analysis restricted to SCZ cases. Toxicology data were used to classify individuals as antipsychotic-exposed or -unexposed. Our models were following

1) Antipsychotic exposed vs unexposed in dopaminergic neurons from SCZ/BD:
Dopa neurons expression: ~ Antipsychotic status + RIN + (1 | Sex) + scale(age) + scale(age^2) + scale(PMI) + (1 | adapters) + GC_NC_20_39 + dopa proportions (cluster 7)

2) Antipsychotic exposed vs unexposed in non-neurons from SCZ/BD:
Non-neurons expression: ~ Antipsychotic status + RIN + (1 | Sex) + scale(age) + scale(age^2) + scale(PMI) + GC_NC_20_39 + Non-neuron proportions (cluster 0)

remaCorr analysis of differences in Tx expression. To reduce the impact of multiple testing in transcript-level analyses, we utilized omnibus testing frameworks grounded in meta-analysis. Specifically, we implemented fixed-effects (FE) meta-analysis for correlated test statistics. To capture transcript-level dysregulation beyond overall gene expression changes, we used the R-based tool remaCorr (Random Effects Regression for Transcript-Level Expression 58; specifically the mvTest() function. remaCorr applies a fixed-effects meta-analytic model to aggregate transcript-level effect sizes (log2FC for SCZ vs. controls or BD vs. controls in this study) for each gene, and tests whether at least one transcript deviates from the overall gene-level pattern. Crucially, remaCorr does not test transcripts individually. Instead, it generates a single p-value per gene, reflecting whether there is significant heterogeneity in transcript usage between groups. This design enables detection of relative isoform abundance changes, even when total gene expression remains stable, capturing regulatory events such as isoform switching, alternative promoter usage, or splicing differences. A similar methodological framework has been applied in Kosoy, R. et al. 29,30

Input to remaCorr, we used the transcript-level fit object derived from the voom transformation and linear modeling using the limma framework with the DREAM method from the variancePartition package, which accounts for complex random effects.

The remaCorr test was run using the following command:
variancePartition::mvTest(fit = fit_Tx, vobj = vobjDream_Tx, features = gene_features, coef = "DxSCZ - DxControl", method = "remaCorr ")

HiC library preparation and sequencing; Hi-C comprehensively detects genome-wide chromatin interactions in the cell nucleus. Considering the intricate packaging of two meters of DNA within a nucleus with a diameter of 2 micrometers, the nucleus is compartmentalized into territories occupied by single chromosomes. Zooming into these territories reveals two compartments: active A-compartments and inactive B-compartments. These compartments consist of numerous interaction clusters known as TADs, distributed throughout the nucleus, where chromatin of the same compartment type is associated. Additionally, loops, which are specific points of interaction between distant chromatin regions, further contribute to the three-dimensional organization of the genome.

Hi-C Nuclei Isolation and Fixation: An aliquot of tissue was dissected from each brain sample from the substantia nigra (A10) area of the ventral midbrain, and dounced 50-70x in 5ml of lysis buffer (0.32M sucrose, 5mM CaCl , 3mM Mg(Ace) , 0.1mM EDTA, 10mM Tris pH8, 0.5mM DTT, 0.1% Triton X-100) until completely homogenous. The sample was transferred to a 15ml conical tube and 5ml of sucrose buffer (1.8M sucrose, 3mM Mg(Ace) , 0.5mM DTT, 10mM Tris-HCl pH8) was added. The sample was inverted and mixed to incorporate the sucrose, then centrifuged for 10 minutes at 2300x g at 4ºC. After removing the supernatant, the sample pellet was resuspended in 4ml of lysis buffer, formaldehyde was added to a final concentration of 1%, and the sample was incubated on a rotator for 5 minutes at room temperature. After that, glycine was added to a final concentration of 217mM, and the sample was again incubated on a rotator for 5min at room temperature. Next, the sample was centrifuged for 5min at 2,300x g at 4ºC, and the supernatant was removed. The sample was then resuspended in 5ml of lysis buffer, and transferred to an ultracentrifuge tube (Polypropylene Centrifuge Tubes ⅝ x 3 ¾ in., Ref. 361707; Beckman Coulter, Indianapolis, IN), underlayed with 9ml of sucrose buffer (1.8M sucrose, 3mM Mg(Ace) , 0.5mM DTT, 10mM Tris-HCl pH8), and ultracentrifuged in a SureSpin 630 (17mL) Rotor (Cat. 79363; Thermo Fisher Scientific, Asheville, NC) with 24,000rpm (106,803x g) for 1hr at 4ºC. After removing the supernatant post-centrifugation, the nuclei pellet was resuspended in 1ml of 0.1% BSA, passed through a 40µm cell strainer, and transferred to a sterile 1.5ml tube. 2µl of Anti-Nurr1 antibody (Anti-Nurr1 antibody produced in rabbit, Cat. N4663 or Anti-NR4A2 (N-terminal) antibody produced in rabbit, Cat. N6413; Sigma-Aldrich, St. Louis, MO) conjugated in a 2:1 ratio to secondary antibody Alexa Fluor 647 (Goat anti-Rabbit IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor™ 647, Cat. A21245; Invitrogen, Waltham, MA) and 1.5µl of NeuN neuronal antibody pre-conjugated to Alexa Fluor 488 (Anti-NeuN Antibody, clone A60, Alexa Fluor™ 488 conjugated, Cat. MAB377X; MilliporeSigma, Burlington, MA) were then added to the sample, and the sample was incubated on a rotator at 4ºC for 2hrs while protected from light. After 2hrs, 1µl of nucleophilic dye DAPI (4',6-Diamidino-2-Phenylindole, Dihydrochloride, Cat. D1306; Invitrogen, Molecular Probes, Inc., Eugene, OR) was added to the sample. A BD FACSAria Cell Sorter (BD Life Sciences, San Jose, CA) was used to sort approximately 6,000 nuclei each of NeuN+/Nurr1+ and 50,000 NeuN-/Nurr1- nuclei populations into collection tubes coated with 5% BSA and pre-filled with 300µl of DPBS. After sorting, samples were centrifuged for 5min at 2300x g at 4ºC. The supernatant was removed, and the pellet was stored at -80ºC. In our hands, 5000-50,000 nuclei were ideal as input for high quality Hi-C libraries with the protocol listed below.

Arima Hi-C library prep and sequencing: Post-FANS sorting, approximately 6,000 dopaminergic neuronal nuclei (NeuN+/Nurr1+) and 50,000 sorted glia (NeuN-/Nurr1-) were processed through the Hi-C method using the Arima-HiC Kit User Guide for Mammalian Cell Lines (Cat. A510008; Arima Genomics, San Diego, CA). Subsequently, they were purified using the Beckman Coulter SPRIselect Beads (Cat. B23318; Indianapolis, IN). Then, the samples were sonicated utilizing the Covaris S220 (Woburn, MA), targeting 300-500 base pairs fragment sizes, then size selected and purified using Beckman Coulter SPRIselect Beads to a target of 300–500 base pairs. The samples were then biotin enriched using the Arima-HiC Kit for Library Preparation, alongside the Swift Biosciences® Accel-NGS® 2S Plus DNA Library Kit (Cat. A101020; Arima Genomics, San Diego, CA). Afterward, the Swift Biosciences Accel-NGS 2S Plus DNA library kit (Cat. 21024; Ann Arbor, MI) was utilized for end repair and adapter ligation. Unique indices were ligated to each sample from the Swift Biosciences 2S Indexing Kit (Cat. 26148; Ann Arbor, MI). DNA libraries were amplified and purified using the Kapa Hyper Prep Kit (Cat. NC0709851; Kapa Biosystems, Wilmington, MA) and Beckman Coulter AMPure® SPRIselect Beads. Final libraries were sequenced through Illumina NovaSeq 6000 technology at the New York Genome Center at 150bp paired-end sequencing, at 400 million reads per sample. We have assembled a dataset comprising 34 HiC libraries, originating from FANS-sorted Nurr1+/NeuN+ immuno-tagged midbrain dopaminergic neurons, as well as their neighboring glial nuclei (Nurr1-/NeuN-). These 34 HiC libraries are categorized into 17 samples of dopa neurons and 17 samples of non-neuronal cells (Fig. 1A, Supplementary Data 6).

Alignment and preprocessing: In our exploration of the spatial organization of differentially regulated genes, we utilized the HiCpro (v3.1.0) pipeline59. We performed alignment using the UCSC’s hg38 genome and generated non-duplicated valid interaction pairs, and we recorded the genomic interactions reported by ligated reads. The quality of analysis results in each step was evaluated including inter-/intra-chromosomal pairs, chimeric pairs, duplicates, intra-fragment, intra-long distance ranges, and ligations. The interaction matrices were generated at different bin sizes, like 20 kb, 40kb, 150kb, 500kb, and 1mb. We used HicRep (v1.12.2)36, QTLtools (v1.3) 60, and PCA for quality check, outlier detection and to find mislabeled samples (Fig. S11).

Following the initial processing steps, we transitioned to the transformation of HiCpro output into a more suitable format for downstream analysis. Utilizing the hicprojuicebox tool (a hicpro utils), we converted the HiCpro output into .hic format, which was then further processed into .cool format using hicConvertFormat (v2.2.1.1) 61. This allowed for seamless integration with the HiCExplorer toolkit, where we proceeded to normalize interactions using hicNormalize and apply correction methods with hicCoIrrectMatrix (v2.2.1.1) 61, including the Knight–Ruiz balancing algorithm with default parameters. These steps were essential for preparing the data for subsequent analyses of chromatin organization and gene regulation in dopaminergic neurons.

For identifying TADs, we utilized TopDom 37, while compartment calling was performed using CscoreTools (v1.1). Additionally, loop calling was conducted using HICCUPS (juicer_tools_1.22.01) 37,39,62.

Loop discovery: Chromatin loops were detected using the juicer tool HICCUPS at resolutions of 5, 10, and 25 kb, employing default window sizes. Only loops that were reproducible across two resolutions were kept and advanced to the next step of analysis. It's important to note that our analysis specifically targeted chromatin loops within chromosomes, disregarding interchromosomal interactions to focus solely on intrachromosomal interactions.

Compartment discovery: To identify chromatin compartments, we used CscoreTool 39, which performs eigenvector decomposition of normalized Hi-C contact matrices to assign compartment scores (PC1 values). Positive PC1 values indicate active (A) compartments, while negative values represent inactive (B) compartments. Because eigenvector-based compartment scores can arbitrarily flip sign between chromosomes or runs, we implemented a sign-alignment step based on transcriptional activity. For each chromosome, we calculated the mean expression of genes overlapping positive and negative compartments (CPM > 0.5 in ≥ 10% of samples). If the negative-score regions showed higher mean expression, we inverted the sign of that chromosome’s compartment scores, ensuring that positive values consistently represent transcriptionally A compartments and negative values represent B compartments. This correction was applied separately to dopaminergic and non-neuronal datasets.

Valid interaction pairs were generated using HiC-Pro, and compartment scores were computed separately for each cell type at multiple resolutions (10 kb, 30 kb, 50 kb, 60 kb, 70 kb, 80 kb, and 100 kb). For all downstream analyses, we used the 50 kb resolution matrices, which provided an optimal balance between resolution and signal robustness. Each genomic bin was assigned to compartment A or B based on the (sign-aligned) PC1 score, and the absolute value of the score was retained to represent compartment strength. To assess cell-type-specific compartment switching, we compared compartment calls between Dopa and non-neuronal cells and classified bins into four categories:
(1) A in both cell types,
(2) B in both cell types,
(3) A in Dopa and B in Non-neurons,
(4) B in Dopa and A in Non-neurons.

TAD discovery: For TAD identification, we implemented a custom R-based workflow built around TopDom v0.0.2 37. First, .cool format HiC matrices were converted to full chromosome-wise contact matrices at 10kb, 50kb, and 100kb resolution using the HiCcompare v 1.20.0 (63) and cooler2bedpe 64 tools. For each chromosome (chr1–22), contact matrices from all samples were summed bin-by-bin to generate deep coverage consensus matrices for each cell type. We then ran TopDom separately on these merged matrices using a window size of 60, 70, 80 , 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 200 (optimized empirically), identifying domain boundaries based on local minima in contact frequency (“binSignals”) Our analysis revealed that a contact matrices at 10kb with window size of 50 provided the best results (Fig. S13A,B). This process resulted in 5,966 TADs in dopaminergic neurons and 6,348 TADs in non-neuronal cells.

Disease trait enrichment analysis using MAGMA: We utilized MAGMA (v1.08b) 23 to perform trait enrichment analysis for both brain-related and non-brain-related traits. The gene lists used for this analysis were derived from limma and remaCorr analyses, as described in the previous sections. Our goal was to assess the enrichment of these GWAS traits in cell-region- and disease-specific genes. MAGMA calculates gene-level p-values for each gene and trait by jointly analyzing the association of all SNPs from GWAS summary statistics that map to the gene region, while accounting for linkage disequilibrium (LD) between SNPs. Gene regions were defined as 35 kb upstream and 10 kb downstream of the gene body. LD was estimated using data from the European panel of the 1000 Genomes Project phase 3 65. Finally, using a linear regression framework, MAGMA tested whether differentially expressed genes showed stronger associations with GWAS traits compared to the rest of the genome.

Functional pathway analysis: Functional pathway analysis was done using gene set enrichment analysis (GSEA) 66 , specifically using the fisher's exact test as the method for the test. The gene lists used for this analysis were derived from limma and remaCorr analyses, as described in the previous sections. The resulting gene set enrichment analysis identifies key biological processes, pathways, or functional gene groups in which the input differentially expressed genes are significantly enriched.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Singh et al provide a comprehensive transcriptomic and chromatin conformational analysis of midbrain dopaminergic neurons in schizophrenia (SCZ). The dataset, including RNA-seq and Hi-C libraries, is robust and represents a valuable resource for understanding the molecular underpinnings of SCZ. The study is well-structured and addresses an important gap in our knowledge of SCZ pathophysiology beyond cortical regions. The Hi-C data is an exciting aspect of this study. However, there are notable issues in the statistical rigor and data presentation that must be addressed before publication.\n\nMajor Comments\n1. Differential Expression Analysis and Confounders\no The differential expression (DE) analysis is concerning due to the relatively small sample sizes for case-control comparisons and the lack of includiont of potential impact of biological confounders (e.g., age, and ancestry) that were not included in the confounders for analysis. Include a detailed breakdown of the sample demographics for each comparison. Supplementary Table 1 does not provide sufficient clarity on the characteristics of included samples. For each comparison, present a table summarizing:\n  Sample size (n)\n  Ancestry composition (e.g., % European)\n  Sex distribution (% Male)\n  Age at death (mean ± SD)\n  PMI (mean ± SD)\n  RIN (mean ± SD) for both original tissue and sorted nuclei.\no Age and genotype principal component analysis (PCA; as a surrogate for ancestry) were not included in the DE models. These are known confounders, especially given the significant age disparity between SCZ donors (mean age: 56.3) and controls (mean age: 43.25). Reanalyze the data including these variables in the models. This is needed for all analysis in the manuscript.\n2. Impact of Antipsychotic Treatment\no Antipsychotic medications, which target dopaminergic systems, are konwn to significantly alter the transcriptome and likely the chromosomal architecture. This is a critical biological effect that was not discussed, addressed or appropriately analyzed in the manuscript. The authors must conduct a stratified analysis comparing antipsychotic-treated and untreated SCZ cases to assess the influence of treatment on the findings. Discuss these results in the context of prior studies on antipsychotic-induced transcriptional changes.\n3. Expression Plots\no The manuscript would benefit greatly from expression plots for the top differentially expressed genes (e.g., CNTNAP2, NRXN1, SHANK2). These should show individual data points and group distributions to allow readers to visually assess the robustness of the findings. Include violin or scatter plots for the most significant DE genes for each comparison (e.g., SCZ vs. controls in dopaminergic neurons).", "reviewer_2": "This manuscript by Singh and colleagues generates RNA-seq and Hi-C datasets from midbrain, dopaminergic positive neurons in human brain (control, SCZ and BD). The data just from the control brains could serve as a resource alone because as the authors point out snRNA-seq from this area currently does not have the same depth as the profiling they did here. The comparisons with the SCZ and BD tissues are somewhat of a separate story and provide a list of differentially expressed genes in these neurons in SCZ.\nOverall, the authors have provided important datasets for the community. I have some comments that should be addressed to increase the impact of the study.\n1) Interestingly, the authors do not find any genes changed in BD. One wonders if the BD dataset suffers from a lack of power since there are half as many BD donors. Could the authors comment on this or carry out some sort of downsampling of the SCZ datasets to see if this is a power issue.\n2) The authors do not find any differentially expressed genes among the non-neurons in SCZ. This is interesting and surprising. Is it possible this is because the non-neurons encompass a variety of cell types and thus any differential signal gets washed out by the cell type heterogeneity? Could this also be due to region of interest? The authors later compare their results with prior results in other cortical areas in SCZ. Were there non-neuronal genes changed in those other datasets from other brain regions?\n3) The rationale to reduce the FDR to 0.1 when there are hundreds of differential genes uncovered using 0.05 seems baffling. Could the authors please provide a more cogent rationale for this choice for the downstream analyses?\n4) In the differential expression model, why are the cell type and sex covariates treated as random effect factors in the model?\n5) For the FANS sorting, why are the number of dopaminergic neurons RNAseq libraries n = 53, and non-neurons RNAseq libraries n = 58 if the same donors were used? Were some libraries dropped for QC reasons?\n6) What is the cell proportion of non-neurons and dopaminergic neurons according to the FANS? Does the proportion change in SCZ and BD?\n7) In line 654, approximately 50,000 nuclei each of NeuN+/Nurr1+ and NeuN-/Nurr1- nuclei populations…; In line 785, Post-FANS sorting, approximately 6,000 dopaminergic neuronal nuclei (NeuN+/Nurr1+) and 50,000 sorted glia (NeuN-/Nurr1-) were processed…, why is the dopaminergic neuronal nuclei (NeuN+/Nurr1+) number so different between RNAseq libraries and Hi-C library prep?\n8) In Fig S6 A, is this for Non-neurons or Dopaminergic neurons in BD?", "reviewer_3": "In this manuscript, Singh and colleagues present a dataset of transcriptomic and chromatin conformation data of midbrain dopaminergic neurons from control individuals and from patients affected with schizophrenia (SCZ) and bipolar disorder (BD). The authors exploited fluorescence-activated nuclei sorting (FANS) to focus the study on a specific and highly relevant population of neurons, bypassing the limitations of bulk tissue-level analysis. In this population, the authors revealed a downregulation of genes associated with synaptic functions in SCZ. Downregulated genes were also associated with genetic risk factor of SCZ from GWAS studies. Transcript-level analysis also identified genes with complex isoform level dysregulation in SCZ. Exploiting their Hi-C dataset, the authors also characterized the topologically associated domains (TADs) of dopamine neurons. The authors observed that transcripts dysregulated in SCZ and present in TADs are also enriched in SCZ genetic risks elements, suggesting a link between chromatin conformation and transcript regulations in genes with high genetic risk association with SCZ.\nThe cell type level multi-Omic dataset presented by Singh and colleagues is an important contribution to the field of psychiatric disorders. It explores the understudied potential alterations associated with SCZ in the midbrain dopaminergic system, with less focus on BD. Several analyses are performed to exploit the dataset and link genetic risk, transcripts and chromatin conformation together. This allowed the authors to make several interesting observations.\nThere are however several limitations and weaknesses to the study and to the manuscript. Authors should consider further demonstrating their purification of dopamine neuron central to the dataset. Further efforts need to be made to clarify the analytical approaches and more systematically report all results across all analyses. The manuscript itself suffers also from several discrepancies between text, legends and figures, raising questions on some observations. Few statements need additional analyses, notably for the TAD-related section. These issues overall undercut the strength and impact of the observations made.\nTherefore, in my opinion, this study requires additional analysis and clarification before it can be published in Nature Communication. Comments and suggestions are detailed below.\n1. Regarding FANS validation\n1.1 The FANS approach to obtain dopaminergic neurons relies on antibodies of the known markers Nurr1 and NeuN. This FANS is essential to the study presented and the dataset generated. Therefore, more efforts need to be done or presented to validate this sorting approach, notably to ensure that the sorted fractions consist of pure dopaminergic neuron. FANS is sensible to contamination due to aspecific binding and the specificity of the antibodies is not demonstrated formally in this study. The authors need to either evaluate it in their report or refer to other published works which estimated used antibodies' specificity in human brain cells. Both Nurr1 and NeuN+ fluorescence profiles do not seem clearly bimodal (based on example in Fig. S1), which could create a certain subjectivity in gating and variations across samples or operators. Currently the demonstration that the sorting was efficient principally relies on 2 genes presented (TH and Olig2, Fig. 2C) that while important to show do not measure efficiency, purity or consistency of the fractionation.\n1.2. It is unfortunate that a comparison with Nurr1-negative/NeuN-positive neurons was not performed. Notably this weakens any statement related to a dopaminergic-specific results, as other neurons could have yielded similar SCZ-alterations, but were not evaluated.\n1.3. Please check reference [19] in \"(NR4A2), a transcription factor crucial for maintenance and development of dopaminergic neuron [19]\" which corresponds in the document received to \"Sudden death in Schizophrenia\" (Li et al 2018) a review that does not mention the gene or the FANS method. Similarly, reference [1] cited after NeuN sort don't seem to mention NeuN or Nurr1 FANS.\n1.4. Two different antibodies seemed to have been used for Nurr1 (Methods) which could have an impact. Did the author evaluate potential effects of chosen antibody. It should be reported in sample metadata.\n1.5. The authors could have performed a single-nuclei sequencing on a few samples after FANS to demonstrate the purity of the Nurr1+ fraction studied and evaluate the composition of the NeuN-negative fraction.\n1.6. Alternative would be to perform single cell-based deconvolution of bulk RNA-seq. This could rely on available single nuclei dataset of human midbrain, midbrain substantia nigra or midbrain dopamine neurons (Adams et al 2024, Agarwal et al 2020, Kamath et al 2022, among others). Sample-level cell type deconvolution results might reveal potential contamination of cells in certain samples or fractions (dissection errors, FANS issue, blood cell contamination in NeuN-, other type of neurons in NeuN+, etc.). At the very least, authors could evaluate and present several specific markers of midbrain cell types expression across samples and diagnosis to further demonstrate the robust dopaminergic nature of the fraction.\n1.7. If not already tested, cell type-specific markers enrichment could be systematically done in parallel to GO pathways enrichment for DEGs. Excluding non-neuronal cell contamination is notably relevant as the authors note that upregulated traits in SCZ relate to immune disease and coronary artery (Table S3H).\n1.8. Subtypes of dopaminergic neurons have been reported (Kamath 2022) and differential capture of those subtypes could explain variations between samples or be associated with diagnosis. Have the authors check potential dopa neurons' subtype markers?\n1.9. The authors made this discussion point: \"Furthermore, our datasets on control brains provide an important complement to other reference sets in the field such as those that profiled human midbrain transcriptomes with single nucleus resolution at the expense of much shallower profiling with typically < 2000 expressed genes called per cell type, or 1-2 orders of magnitude below the depth accomplished here\" This statement warrants a focused comparison between gene expression profiled by FANS+RNA-seq vs snRNA-seq results demonstrating the advantage of the dataset. Could DEGs described in the study be evaluated by snRNA-seq assays?\nGeneral comments and questions\n2. More generally, can the authors assess or discuss if there was any evidence of clinical or transcriptomic heterogeneity between SCZ/BD patients studied that could have confounded the differential expression results? The authors also evaluated the variations dependent on Sex but didn't discuss potential sex-specific effects on the disorders profiled.\n3. The leniency cutoff of 0.1 applied for Fig.S5 seemed unnecessary and should be avoided. There were already hundreds of DEGs and MAGMA enrichment at FDR < 0.05. It seemed only included in the manuscript to obtain MAGMA enrichment for SCZ and BD GWAS. The authors should refrain from mentioning associated results in the abstract or main conclusions as this could be misinterpreted as a strong result. Larger cohort might be needed to reach those conclusions.\n4. Importantly, more effort needs to be done on the reporting of results with statistical leniency in all figures. Both cutoffs should be presented in the same figure systematically for easier comparison. Since different cutoff are applied, the cutoff applied and to which test (DEG, RE2C, enrichment, MAGMA, etc.) need to be systematically visible in every panels.\n5. The authors note a MAGMA enrichment for SCZ genetic risk in cortical tissue upregulated DEG, which contrasts to its nominal enrichment in downregulated DEG in Midbrain neurons. Is there an overlap in the genes/genetic risk association between the regions or do they represent different subsets of the SCZ genetic risk?\n6. The comparison of DEG results from homogenates cortical tissue (i.e. multiple cell types) and sorted cell type (Nurr1+ neurons) is difficult to interpret. The first one represents an aggregated result of both glial and neuronal cells DEGs and their relative proportions. A counter argument would be that if correctly sorted, certain subtype of cortical neurons could exhibit the similar phenotype identified in dopaminergic neurons. Therefore, a comparison with a NeuN+ FANS done in the SCZ cortex would reinforce the statement that dopa neurons exhibit a regional-specific downregulation in SCZ. Notably, the statement in Discussion: \" However, it is important to note that according to our analyses, the disease-sensitive transcriptome of midbrain dopa neurons is significantly anti-correlated with its two target sites (PFC and striatum (Fig. S7 D-F)\" seems uncorroborated at this point due to scale difference (tissue- vs cell-level DEG).\n8. Points regarding RE2C transcript analysis\n- RE2C analysis is not trivial to understand and interpret. It seems to return statistics at the gene level only in Fig.4E, but Methods states: \"For each transcript, RE2C provides a statistical test and outputs a p-value reflecting the significance of the observed differences\". The authors should state clearly what RE2C actually test. Does it test for a relative change of any transcript in the overall diversity of transcripts present for a given gene regardless of the gene expression change? Relative transcript abundance of a gene could change across conditions with or without affecting the overall expression level of the gene itself.\n- The principle of the RE2C analysis seems to have been proposed and applied for GWAS and cohorts overlap according to the reference given by the authors. If this is the first application of RE2C to transcript isoform analysis, the authors should expand the description in the methods for reproducibility, explain why this approach can be used for isoform/transcript measurement, and how it compares with current existing approach.\n- Related, have the authors tried other approach used for differential transcript usage, such as LeafCutter, exon-level and percent spliced-in (PSI) analyses/reporting?\n- The term \"discordant genes\" does not seem adapted since it's the transcripts that are discordant, not the gene.\n- Could the author report changes in transcripts relative abundances alongside overall fold change of the gene? Would it be possible for instance to classify genes by overall DEG significance and report differential transcript usage between conditions significant by RE2C in Fig.4E? How is the distribution of log2FC computed in Fig.4E?\n- The authors note higher number of Tx and longer size for genes with discordant transcripts but do not comment on it being a statistical effect or a biological result. Should the analysis somehow correct for gene's length or number of exons to limit this inflation?\n- The discrepancy between RE2C and limma significance are puzzling. One would expect that genes with concordant transcripts would meet significance at the gene level. Fig. S8B should be separated into concordant vs discordant genes.\n- Did the authors notice changes in relative transcript abundances across individuals?\n- The author performed RE2C analysis in non-neuronal fraction, but the results are not presented or commented. Same remarks for BD results? Are the distribution of concordant/discordant genes obtained for SCZ similar in BD?\n- There are issues with the reported gene numbers for concordant up down/genes related to panel A and B of Figure S9. The reported gene numbers in x axis labels (104, 28, 97) are different between Figure S9 panels A and B axes, but also in the legend and the main text!! Panel A legend reads \"GWAS traits in RE2C detected concordant upregulated (n=28), concordant downregulated (n=104), and discordant SCZ-specific genes (n=97) at FDR<.05\" and Main Text reads: \"132 concordant genes (104 concordant upregulated and 28 concordant downregulated) and 97 discordant genes\". The authors should double-check numbers and correct label throughout manuscript and verify the follow-up analyses (MAGMA, functional pathway enrichment, TADs, etc.) have been performed on the correct sets and the related statements.\n- \"MAGMA analysis at strict (FDR<.05) and lenient threshold (FDR<.10) revealed that both concordant downregulated Tx genes and discordant Tx genes were significantly associated with brain-related traits compared to other traits (Fig. S9A, Fig. S10B)\". Results with those 2 FDR thresholds should be presented side by side in a unique figure with clear indications of FDR and number of genes obtained. There seems to be no significant enrichment for concordant downregulated Tx in either panel. Please check statement.\n9. Points on Hi-C and TAD analyses\n- Half of the Hi-C samples come from BD donors. This is not mentioned or properly evaluated. BD effects should be assessed and affected TADs removed for SCZ interpretation. It's unfortunate that no SCZ Hi-C was collected, making SCZ vs. TAD observations weaker.\n- The TAD methods section is unclear. Notably it mentions evaluating consensus TAD across different fractions of the samples without precising the final selected number or a final strategy. Was this consensus TAD process affected by BD samples?\n- Statement from previous work from the same group in (Girdhar et al. 2022): \"Interestingly, acetylated CRDs consistently showed, both in H3K27ac NeuN+ (Study-1) and H3K37ac Tissue (Study-2), maximum density at the center of sub-TADs and TADs (Supplementary Fig.11). In contrast, methylated CRDs were enriched at TAD boundaries (Supplementary Fig.12), resonating with earlier reports on H3K4me3 enrichment at TAD boundaries.\"\nStatement from this study: \" Consistent with studies in other brain regions and cell types, TAD boundaries in dopa neurons were enriched for CTCF chromosomal scaffolding protein binding sites and for transcription (promoter-enhancer) histone marks H3K27ac, while showing significant depletion for some of the repressive histone marks including H3K27me3 (Fig. S12).\" There seems to be a discrepancy between the 2 studies' statements here. An enrichment of H3K27ac at TAD boundaries would be counterintuitive. Please comment, clarify or correct the statement.\n- What is the origin of the histone marks peaks used in Fig.S12? Was it obtained from ChIP-seq of the same FANS fractionation, or midbrain bulk tissue?\n- Figure S12B: Discrepancy between the legend and left panels title mentioning \"H3K27ac\" for dopa neurons whereas non-neurons (right) and legend report h3K4me2. The text mention H3K27ac. Please show all marks for both cell type. Grey line is not described in legend.\n- In example region shown in Fig 5 A, the TADs and compartment boundaries seem to not align and not align with the contact map matrix triangular domains. Please explain or comment the discrepancies. One TAD boundary (marked by space in the green bar?) seem to cut the PDE48 gene in half despite its strong expression by RNA-seq signal? Similarly, the 2nd boundary corresponds to a region with high contact signal on the top heatmap. Similarly loops overlap multiple TADS. Wouldn't these discrepancies be unexpected from the TAD analysis and definition as an insulating unit?\n- \"Consistent with the previous studies, 98% of TADs showed overlap across cell types.\" This statement should be clarified and is not corroborated by any panel. The \"studies\" should be cited. Notably the authors already obtained different numbers of TADs from the 2 cell types evaluated (N=5,966 and N=6,348), making 98% unclear.\n- \" we also identified cell-specific ‘A’ active and ‘B’ inactive compartments (a total of 65,639 compartments) using Cscoretools -based eigenvector analyses\" Similarly, this is not described in a panel or measured? How many and which regions are differentially active between cell types?\n- Figure 5B evaluates genes inside or outside TAD across \"All genes\" and across \"DEG\". Does \"all genes\" include both expressed and not expressed genes? Is there otherwise a link between expression level and inside/outside TAD status? The result also shows significant enrichment of SCZ GWAS traits across all genes? Only a minority of genes seems to have an associated TAD, is this number in range with previous studies?\n- \"and asked whether these genes reside predominantly within the TADs. Indeed, downregulated SCZ-specific genes within TADs showed significant enrichment for SCZ risk genes in dopa neurons.\" The 2 sentences are disconnected since the first sentence's question is not answered. Are the SCZ-DEG genes enriched in TAD compared to other expressed genes in dopaminergic neurons?\n- If available from previous work and not described in other studies, I would encourage the authors to compare HiC/TAD results in dopa neurons to those obtained for other neurons such as cortical ones (rather than to glial cells).\n- \" Importantly, the proportion of dysregulated SCZ dopa neuron genes that extended across multiple TADs (Fig. S9.E) was significantly higher for the group of discordant (28/97) as compared to concordant (10/132) genes (Tables S5H-I).\". Related to above comments it is highly unclear how the gene body of an expressed gene in a purified cell type can overlap several TADs identified in the same cell type? It would suggest that the TAD definition used in this study is too lenient to truly represent large strongly independent regulatory domain? This could depend on a different way of interpreting TADs on my part, but this needs further clarification on how TAD are defined and their sensitivity and relationship with gene expression. This weakens the TAD interpretation.\n- \"This speaks to the complex regulome for this group of transcripts, which included a number of large, high penetrance schizophrenia risk genes including the extensively spliced synaptic adhesion molecules, NRXN1 [36,37], CNTNAP2 [38], SHANK2 [39,40] and CTNND2 delta-catenin [41].\" This is indeed highly interesting results, and I encourage the authors to present a detailed view of the genomic region of those genes as in Fig.5A, alongside alternative transcript measured and their differential transcript measurement. Did the author check if downregulated vs upregulated alternative transcripts of a same gene spanned a TAD boundary?\n- Discussion point: \" In the present study, we specifically focused on Mb-scaling TADs, and the smaller sub-TADs hierarchically nested into them, which represent a type of conformation and cis-aligned transcriptional regulome defined by dynamic chromosomal loop extrusions of individual chromatin fibers, constrained by strong boundary elements at TAD peripheries and weaker in-TAD boundaries demarcating sub-TADs\" At no point do the author mention or do a distinction between TAD and sub-TAD in their Results. This needs to be more clearly presented throughout the results.\nOverall, the TAD analysis needs to be more systematic with a clear definition of TAD/sub-TAD and their relationship with gene locations (most genes have no identified TAD) and gene expression and dysregulation. If loops, compartments and contact measurement can span multiple TADs, it needs to be clarified. Enrichment in TAD/subTAD vs outside TAD need to be systematically reported (before MAGMA analysis).\nMinor comments\n- How do the authors explain variation across Dopa neuron samples spread across PC2 in Fig.S2B?\n- Verify the boxplots IQR rule for outlier definition mentioned in legends. My understanding of default boxplot plotting parameters is that outlier lies further than 1.5x IQR from the 25/75 quantiles, not from the median.\n- Indicate actual p-value in S4A as the color scale is not visually evident from the lines color.\n- Disorders Abbreviations in S4B and all subsequent figures should be given at least once in a panel or legend for the reader.\n- Similarly, ACC /PFC/DLPFC abbreviations are not spelled out nor in the cited reference.\n- Figure S6A name the BD significant DEGs. Overall, top significant genes could be indicated on many plots.\n- Figure S11 lacks x axis, or at least separation of samples into clinical groups for this QC report to be appreciated. (panel D has a typo in \"correlation\")\n- Fig 5B upregulation / downregulation indications near x axis label should be made clearer than arrows to avoid any possible confusion on which is which.\n- Fig 4D typo in second panel title: \"Cocornordant\"\n- \"We then asked if the directionality of SCZ transcriptome is similar in other cortical regions.\" This sentence should be modified as it suggests that midbrain is a cortical region.\n- The axis labels (logCPM, Eigenvector, enrichment, FDR?) of each part of the panel Fig.5A are lacking, A and B compartment should be clearly colored and TAD boundaries made more evident.\n-\" much more strongly associated with neuronal signaling (GO) pathways compared to those positioned outside TADs (Fig. 5B, C).\" The figures do not report the outside TADs result. The panel C enrichment should name the 2 to 4 genes overlapping each GO term.", "reviewer_author_response": "General Comments\nWe thank the editor and reviewers for their constructive feedback on our manuscript. We have\ncarefully considered each comment and have undertaken substantial revisions to address the key\nconcerns raised. In the revised manuscript, we have conducted following major revisions:\na) Cell-type deconvolution for purity validation: To directly address concerns about potential\ncontamination and cellular heterogeneity, we performed MuSiC-based 1 deconvolution using a\nhigh-resolution single-cell midbrain reference dataset 2.\nb) Updated covariates in differential expression models: In response to reviewer concerns about\npotential confounds, we have reanalyzed all transcriptomic data with updated linear mixed-effects\nmodels that now include age as a covariate, in addition to RIN, sex, PMI, GC content, and\ndeconvolution-derived cell type proportions.\nc) Enhanced quality control: We applied stringent QC criteria to remove outlier samples in\nprincipal component analysis (PCA) space of dopaminergic and non-neuronal samples,resulting in 98\nhigh-quality samples (previously n=111).\nWe believe these revisions have significantly enhanced the robustness of our findings.\nReviewer comments are written in black, whereas authors answers are written in blue. Below,\nwe provide detailed point-by-point responses to all reviewer comments.\nReviewer #1\nSingh et al provide a comprehensive transcriptomic and chromatin conformational analysis of\nmidbrain dopaminergic neurons in schizophrenia (SCZ). The dataset, including RNA-seq and\nHi-C libraries, is robust and represents a valuable resource for understanding the molecular\nunderpinnings of SCZ. The study is well-structured and addresses an important gap in our\nknowledge of SCZ pathophysiology beyond cortical regions. The Hi-C data is an exciting\naspect of this study. However, there are notable issues in the statistical rigor and data\npresentation that must be addressed before publication.\n1.1. Differential Expression Analysis and Confounders\n1.1.1 The differential expression (DE) analysis is concerning due to the relatively small sample sizes\nfor case-control comparisons and the lack of inclusion of potential impact of biological confounders\n(e.g., age, and ancestry) that were not included in the confounders for analysis. Include a detailed\nbreakdown of the sample demographics for each comparison. Supplementary Table does not provide\nsufficient clarity on the characteristics of included samples. For each comparison, present a table\nsummarizing: Sample size (n), Ancestry composition (e.g., % European), Sex distribution (% Male),\nAge at death (mean ± SD), PMI (mean ± SD), RIN (mean ± SD) for both original tissue and sorted\nnuclei.\nReviewer #2 :\nThis manuscript by Singh and colleagues generates RNA-seq and Hi-C datasets from midbrain,\ndopaminergic positive neurons in the human brain (control, SCZ and BD). The data just from the\ncontrol brains could serve as a resource alone because as the authors point out snRNA-seq from this\narea currently does not have the same depth as the profiling they did here. The comparisons with the\nSCZ and BD tissues are somewhat of a separate story and provide a list of differentially expressed\ngenes in these neurons in SCZ.\nOverall, the authors have provided important datasets for the community. I have some comments that\nshould be addressed to increase the impact of the study.\n2.1) Interestingly, the authors do not find any genes changed in BD. One wonders if the BD dataset\nsuffers from a lack of power since there are half as many BD donors. Could the authors comment on\nthis or carry out some sort of downsampling of the SCZ datasets to see if this is a power issue.\nReviewer #3\nIn this manuscript, Singh and colleagues present a dataset of transcriptomic and chromatin\nconformation data of midbrain dopaminergic neurons from control individuals and from\npatients affected with schizophrenia (SCZ) and bipolar disorder (BD). The authors exploited\nfluorescence-activated nuclei sorting (FANS) to focus the study on a specific and highly\nrelevant population of neurons, bypassing the limitations of bulk tissue-level analysis. In this\npopulation, the authors revealed a downregulation of genes associated with synaptic\nfunctions in SCZ. Downregulated genes were also associated with genetic risk factor of SCZ\nfrom GWAS studies. Transcript-level analysis also identified genes with complex isoform level\ndysregulation in SCZ. Exploiting their Hi-C dataset, the authors also characterized the\ntopologically associated domains (TADs) of dopamine neurons. The authors observed that\ntranscripts dysregulated in SCZ and present in TADs are also enriched in SCZ genetic risks\nelements, suggesting a link between chromatin conformation and transcript regulations in\ngenes with high genetic risk association with SCZ.\nThe cell type level multi-Omic dataset presented by Singh and colleagues is an important\ncontribution to the field of psychiatric disorders. It explores the understudied potential\nalterations associated with SCZ in the midbrain dopaminergic system, with less focus on BD.\nSeveral analyses are performed to exploit the dataset and link genetic risk, transcripts and\nchromatin conformation together. This allowed the authors to make several interesting\nobservations.\nThere are however several limitations and weaknesses to the study and to the manuscript.\nAuthors should consider further demonstrating their purification of dopamine neuron central to\nthe dataset. Further efforts need to be made to clarify the analytical approaches and more\nsystematically report all results across all analyses. The manuscript itself suffers also from\nseveral discrepancies between text, legends and figures, raising questions on some\nobservations. Few statements need additional analyses, notably for the TAD-related section.\nThese issues overall undercut the strength and impact of the observations made.\nTherefore, in my opinion, this study requires additional analysis and clarification before it can\nbe published in Nature Communication. Comments and suggestions are detailed below.\n3.1. Regarding FANS validation\n3.1.1 The FANS approach to obtain dopaminergic neurons relies on antibodies of the known markers\nNurr1 and NeuN. This FANS is essential to the study presented and the dataset generated. Therefore,\nmore efforts need to be done or presented to validate this sorting approach, notably to ensure that the\nsorted fractions consist of pure dopaminergic neuron. FANS is sensible to contamination due to a\nspecific binding and the specificity of the antibodies is not demonstrated formally in this study. The\nauthors need to either evaluate it in their report or refer to other published works which estimated used\nantibodies' specificity in human brain cells. Both Nurr1 and NeuN+ fluorescence profiles do not seem\nclearly bimodal (based on example in Fig. S1), which could create a certain subjectivity in gating and\nvariations across samples or operators. Currently the demonstration that the sorting was efficient\nprincipally relies on 2 genes presented (TH and Olig2, Fig. 2C) that while important to show do not\nmeasure efficiency, purity or consistency of the fractionation.\nreviewer Fig.6 already visualizes log2FC for each transcript within select SCZ-associated genes,\nincluding both statistically significant and non-significant features, classified by remaCorr significance.\nThese distributions represent per-transcript effect sizes estimated from the limma model and are\nstratified by gene class (concordant up, concordant down, discordant), illustrating intra-gene transcript\nregulation. The distribution of log2FCs in Fig. 3E is computed directly from limma-derived SCZ-specific\neffect sizes for each transcript (Tx), with error bars showing the standard error of the log2FC estimate.\nRed indicates FDR < 0.05 for individual transcripts, highlighting where specific isoforms are\nsignificantly altered even if the gene-level expression is not. In summary: Transcript-level relative\nabundance changes are already captured via remaCorr, and the results are reported and visualized in\ndetail (Fig. 3, Response to reviewer Fig.6). remaCorrallows identification of genes with differential\ntranscript usage, including those missed by gene-level DEG analyses. Fig.3E/Response to reviewer\nFig.6 presents transcript-level log2FC distributions, contextualized within RE2C significance calls,\naddressing the reviewer’s suggestion to link isoform changes with disease-associated transcript\nremodeling.\nFigure 3. SCZ-specific concordant and discordant transcripts in dopa neurons. (A) Genes X and Y are examples of\ngenes with concordant transcripts, meaning their transcripts are consistently either upregulated or downregulated, as\nindicated by the colors red and yellow. In contrast, Gene Z is a discordant gene, with its transcripts in orange showing both\npositive and negative effect sizes. (B) Barplot showing the number of genes which have concordant (up and down) and\ndiscordant transcript. The crosshatch in each bar represents the number of DEG identified in gene analysis (C) pi1 stats of\ngenes from remaCorr analysis using gene expression analysis. (D) MAGMA gene set enrichment analysis of concordant and\ndiscordant genes in SCZ for enrichment in SCZ genetic risk (FDR < 0.05, Benjamini–Hochberg adjustment across all tests\nshown) (E) Example of significant discordant genes (WWOX) from remaCorr analysis. Summary statistics of a gene from the\nremaCorr model including P Values, remaCorr statistics, and FDR values are added above the plot. The X-axis shows the\ngene and mapped transcripts that are expressed within the gene. The Y-axis shows the SCZ-specific effect sizes from\nlimma-based analyses. Effect sizes of transcripts marked in red have limma-based P values <.05 and non-significant\ntranscripts are shown in black.\nResponse to reviewer Fig.6 | Concordance of transcript-level and gene-level expression changes in SCZ. Scatter plot\nshowing the relationship between gene-level log₂ fold changes (x-axis; FDR < 0.05, limma) and the average log₂ fold change\nof all associated transcripts (y-axis; derived from transcript-level limma model) for SCZ vs. control comparisons. Each point\nrepresents a gene, with transcript-level changes aggregated across isoforms. A strong positive correlation (ρ = 0.81, p <\n2.2×10⁻¹⁶) indicates that transcript-level shifts largely mirror gene-level expression changes, while also capturing additional\nisoform-specific variability.\n3.3.6 The authors note higher number of Tx and longer size for genes with discordant transcripts but\ndo not comment on it being a statistical effect or a biological result. Should the analysis somehow\ncorrect for gene's length or number of exons to limit this inflation?\nreviewer comment 3.4.8). See for example the screenshot of gene NRXN1 and FOXP1 from IGV\ntrack is shown in Response to reviewer Fig.8.\nResponse to reviewer Fig.8 | Genomic organization and dopaminergic neuron TAD structure at NRXN1 and FOXP1 loci\n3.4.13 Discussion point: \" In the present study, we specifically focused on Mb-scaling TADs, and the\nsmaller sub-TADs hierarchically nested into them, which represent a type of conformation and\ncis-aligned transcriptional regulome defined by dynamic chromosomal loop extrusions of individual\nchromatin fibers, constrained by strong boundary elements at TAD peripheries and weaker in-TAD\nboundaries demarcating sub-TADs\" At no point do the author mention or do a distinction between TAD\nand sub-TAD in their Results.\nReviewer #3 (Remarks on code availability):\nThe code is not made available. The authors have relied on existing package and tools, including\nsome developed by them (VariancePartition/Dream, etc.).\nFor reproducibility, the description in Methods could be improved for RE2C and TAD, especially if not\ncode is provided."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have satisfied my concerns. Thank you for the careful response. This is an important contribution to the field.", "reviewer_2": "The authors have been responsive to my prior critiques. I have no further concerns.", "reviewer_3": "In their revision, Singh and colleagues performed a cell type deconvolution that validated the compositions and purity of the sorted dopaminergic neuron fractions. They also improved their statistical models to include additional covariates and corrected all noted reporting issues raised by the reviewers. The authors appropriately modulated their statements following reviewers’ concerns, or in other cases offered convincing counterarguments or evaluations. The effort and depth of the modifications were appreciated, and the revised manuscript now presents methods, results, and discussion with great details.\n\nThe manuscript present an interesting dataset and an important advancement for our understanding of Schizophrenia's impact on midbrain dopaminergic neurons.\n\nI have only a single remaining request to clarify the compartment analysis that was added during revision that, while not central to the report, might warrant authors' verification.\n\nFollowing this, I recommend the publication in Nature Communication.\n\nRegarding Compartment analysis.\n- In Fig.S14, the authors show large compartment switching between cell types.\n- Response 3.4.7 mentions using ChIP-seeker to link them to schizophrenia-relevant genes, but it does not seem presented ? (Since a similar analysis is done at TAD level, this might not be a necessary addition).\n- Seeing the important shift observed, it is important to show that inactive-to-active switching corresponds to transcriptomic change between dopa-neuron and non-neuron.\n- The Method on compartment analysis states: \" To identify chromatin compartments, we used CscoreTool, which performs eigenvector decomposition of normalized Hi-C contact matrices to assign compartment scores (PC1 values). Positive PC1 values indicate active (A) compartments, while negative values represent inactive (B) compartments.\". It seems from looking at CscoreTool methods that it returns a likelihood value (-1 to 1 range), not a PC1 value?\n- More importantly, my understanding from CscoreTool description and github is that obtained sign output of CscoreTool (or Eigenvector-method, for that matter) can flip between run and chromosomes. In consequence for the same chromosome ABBAB and BAABA would be equivalent results simply due to sign flip. In that case CscoreTool attribution to A and B does not necessarily correspond to active and inactive compartment status, which is difficult to infer from Hi-C alone. Additional data are required to infer which of A and B is active compartment, notably transcriptomic activity, as suggested above.\n- If this is a correct caveat, was inactive/active compartment attribution between samples verified without ambiguity? The authors might have considered or verified that aspect already, but the Method description suggest an automatic attribution of negative value to B compartment, and of B compartment as \"inactive\". Importantly, could the author verify that this sign flip is not the source of compartments switching between dopa-neurons and non-neuron?\n- Also, the text mention 65k compartments were identified, but the sum of bars in S14B and Table reports less total bins? Is there a change in definition between compartment and bins, or some compartment excluded for the comparison? This could be indicated in the figure.\n- Finally, it seems that even minor switches between low positive to low negative score were considered as switch in Table S7. Since low compartment score could more easily switch from negative to positive, they could be filtered out as low quality compartments.\n- From a quick look at Table S7, it also seems like there is a chromosome bias in the results?\nAltogether I think this warrants a verification of the Method description and a verification that transcriptomic data match compartment attribution across chromosomes and samples."}, "author_response": ""}]

Decision letter