Disturbances of paraventricular thalamic nucleus neurons in bipolar disorder revealed by single-nucleus analysis

Domain: biology
Published: 2026-01-07
Editorial outcome: minor_revision

Abstract
Bipolar disorder (BD) is a major global health burden, and its treatment challenges highlight the need for pathology-based therapeutic development. Emerging evidence suggests that the thalamus, particularly the paraventricular thalamic nucleus (PVT), is a key region in mood regulation. We performed single-nucleus RNA sequencing on 82 thalamic and cortical samples from 21 patients with BD and 20 controls to compare transcriptional pathology. PVT neurons showed the most striking abnormalities, including the largest number of differentially expressed genes and ~50% fewer cells in BD, whereas cortical alterations were comparatively modest. PVT neurons exhibited marked downregulation of synaptic and ion channel-related genes such as SHISA9, CACNA1C, and KCNQ3, which are linked to BD risk and serve as central nodes in downregulated networks. We also observed disrupted interactions between thalamic excitatory neurons and microglia. Overall, PVT neurons emerge as a central pathological substrate and a promising diagnostic and therapeutic target in BD.

Introduction
BD is a psychiatric condition that affects approximately 1% of the global population. Patients with BD experience recurrent episodes of mania and depression, often accompanied by suicide attempts. Although existing treatments can be effective for certain patients, potential adverse effects and treatment resistance highlight urgent needs for novel therapeutic approaches. Advancements in treatment will hinge on a deeper understanding of the biological mechanisms underlying BD; however, pathological processes within the brain remain elusive. Previous genetic studies have pointed to dysregulation of synaptic and ion-channel-related pathways in BD1–4. Yet, the investigation of detailed mechanisms has been hampered due to a lack of analytical methods to observe cellular and molecular changes in the human brain extensively.

The advancement of single-nucleus RNA sequencing (snRNA-seq) has provided opportunities to comprehensively analyze gene expression in the human brain at single-cell resolution5–7. Emerging studies using snRNA-seq have illuminated cellular and molecular pathologies underlying schizophrenia8–10, autism spectrum disorder11,12, and major depressive disorder13,14. For instance, Ling et al. identified disruptions in the synaptic neuron-astrocyte program in schizophrenia10. However, the pathological mechanisms of BD remain underexplored despite the previous bulk RNA-seq studies15–18. Furthermore, prior snRNA-seq investigations into psychiatric disorders have predominantly focused on cortical regions, particularly the frontal cortex. Subcortical regions warrant further investigation to elucidate the mechanisms underlying psychiatric disorders.

Comprehensive volumetric studies using MRI have revealed significant volume reductions in subcortical regions, particularly the thalamus and hippocampus, in patients with BD19, as well as in cortical regions20. Notably, lithium treatment has been reported to restore thalamic volume, highlighting the thalamus as a potential core pathological region and therapeutic target for BD19. Within the thalamus, the paraventricular thalamic nucleus (PVT) is particularly promising for BD research. PVT serves as a critical hub for emotional processing, with projections to the nucleus accumbens, amygdala, and medial prefrontal cortex (mPFC) in mice21–23, and has been implicated as a key region associated with mood disorders in the Polg transgenic mouse model24 and neural circuit manipulation models24,25. Additionally, the PVT is one of limited regions enriched with serotonin and noradrenaline in mice26, suggesting its relevance to mood modulation. While Schulmann et al. have reported the physiological cytoarchitecture of the medial thalamus by snRNA-seq and spatial transcriptomics27, the thalamic cellular and molecular pathologies in BD warrant further investigation.

Here, we examined 82 postmortem thalamus and cortex samples from 21 patients with BD and 20 controls by snRNA-seq to investigate the basic cytoarchitecture of the human PVT and its alterations in BD. We compared the pathological changes in the thalamus to those observed in the frontal cortex (BA10) to determine which brain region is more closely associated with BD pathology (the analytical workflow is depicted in Fig. 1a). By analyzing both brain regions from the same cohort, we minimized confounding factors related to individual background and specimen conditions, thus enabling a clearer comparison between the thalamus and cortex. Among various cell types in these brain regions, the PVT showed a prominent reduction and dysregulated synaptic/neuronal genes in BD, alongside disturbed synaptic neuron-microglia interactions. In PVT neurons, SHISA9, CACNA1C, and KCNQ3 exhibited notable changes as putative BD risk genes and/or critical hubs in the downregulated synaptic/neuronal gene networks. In contrast, the cortex exhibits less transcriptomic and compositional changes than the thalamus. Our study offers insights into the pathological mechanisms of BD and highlights the PVT as a potential target for future drug development.

Results
A total of 383,130 single-nucleus transcriptomic profiles were obtained from 21 patients with BD and 20 controls, encompassing 183,354 thalamic nuclei and 199,776 cortical nuclei before QC (Supplementary Data 1). The thalamic nuclei were specifically isolated from the medial region of the thalamus (Fig. 1a).

Cytoarchitecture of the medial thalamus
Initially, we focused on the thalamic and cortical snRNA-seq data from 20 controls to delineate the physiological cytoarchitecture of the human medial thalamus. Our QC yielded 32 distinct cell clusters (Fig. 1b). One cluster, designated as NeuCyto, was removed from the subsequent analysis as non-nuclei neuronal cytoplasmic components (see Methods). Accordingly, 147,477 nuclei in 31 cell nuclei clusters underwent subsequent analysis. Excitatory and inhibitory neurons were distinctly segregated between the thalamus and cortex, while astrocytes and oligodendrocyte progenitor cells (OPCs) were continuously separated. This observation aligned with Siletti et al.6.

The thalamic excitatory neurons (ExN_Tha) predominantly expressed SLC17A6 (coding VGLUT2), whereas the cortical excitatory neurons (ExN_FrC) primarily expressed SLC17A7 (coding VGLUT1) (Supplementary Fig. 1a). RBFOX3, encoding a well-established neuronal marker known as NeuN, exhibited robust expression in ExN_FrC, yet was minimally expressed in ExN_Tha. Hierarchical clustering of the 31 cell clusters, using all expressed genes, indicated that thalamic excitatory neurons were more closely aligned with cortical inhibitory neurons than with cortical excitatory neurons (Fig. 1c), corroborating the findings by Mathys et al.28. Furthermore, mapping through cortex reference29 confirmed that thalamic excitatory neurons corresponded more closely with cortical inhibitory neurons than with their excitatory counterparts (Supplementary Fig. 1b). These observed similarities were represented by specific subtype markers of cortical inhibitory neurons, particularly calcium-binding protein genes, including PVALB, CALB1, and CALB2, which encode parvalbumin, calbindin, and calretinin, respectively (Fig. 1d and Supplementary Fig. 1c).

When compared to ExN_FrC, ExN_Tha showed a preferential expression of genes associated with cardiac muscle development, such as SHOX2, PROX1, HCN4, and CACNA1G, as determined by gene set enrichment analysis (GSEA) (Fig. 1e, Supplementary Data 2, 3). CACNA1G, which encodes a subunit of T-type calcium channels involved in electrical pace-making activity, is specifically expressed in the thalamus and the heart30, and has been identified as a core gene associated with schizophrenia31. Among the ExN_Tha-specific genes, CASQ2 was highly characteristic to ExN_Tha, with virtually no expression in other cell clusters (Supplementary Fig. 2a). CASQ2 encodes calsequestrin-2, a high-capacity calcium-binding protein particularly expressed in cardiac muscle cells32, further underscoring the commonality between thalamic excitatory neurons and cardiac muscle cells.

We also assessed the expression of monoamine receptor genes as primary targets for antipsychotics (Supplementary Fig. 2b). Notably, DRD2, HTR7, and HRH3 exhibited pronounced expressions in thalamic excitatory neurons compared to other cell types in the thalamus and cortex. DRD2 is a key target of antipsychotic drugs and is implicated in schizophrenia, a genetically related disorder to BD33. HTR7 represents the primary target of lurasidone, an antipsychotic used in the treatment of both BD and schizophrenia. The notable expression of DRD2 in the thalamus, coupled with its minimal expression in the frontal cortex, underscores the thalamus’s relative significance in psychotic disorders. Figure 1f illustrates the expression of representative genes characteristic to thalamic excitatory neurons.

Inhibitory neurons were classified into three distinct classes based on their source tissues and ganglionic eminence markers compiled from Schmitz et al.34. Thalamic inhibitory neurons expressed lateral ganglionic eminence (LGE) markers including ZFHX4, FOXP2, MEIS2, and TSHZ134, as well as midbrain-derived inhibitory markers such as SOX14 and OTX2, consistent with the findings from Kim et al.35 (Supplementary Fig. 1c). These neurons were annotated as InN_Tha without reference to their developmental origin, which, while intriguing, is beyond the scope of the present study. Compared to cortical inhibitory neurons, thalamic inhibitory neurons also demonstrated a preferential expression of genes associated with cardiac muscle, specifically in atrial cardiac muscle tissue development (Supplementary Data 2, 4). This enrichment, analogous to ExN_Tha Gene Ontology (GO) enrichment, was also driven by the expression of CACNA1G and HCN4.

Cytoarchitecture of PVT neurons
Among the thalamic excitatory neurons, ExN_Tha_CALB cells in Fig. 1b were putatively supposed to be PVT neurons, as characterized by the expression of CALB1 and CALB2. We used these genes as guidance for locating PVT neuron clusters because they have been employed in previous studies as indicators of PVT neurons in both humans36 and rodents37,38. However, note that CALB1 and CALB2 expressions are not exclusive to PVT neurons, and additional heterogeneity likely exists within the human PVT beyond what can be defined by these genes. These putative PVT neurons were also robustly identified in the Human Brain Cell Atlas (HBCA) dataset6, as demonstrated by the integration of our and HBCA datasets (Fig. 2a). ExN_CALB in Fig. 2a corresponded to ExN_Tha_CALB in Fig. 1b. The putative PVT neurons (ExN_CALB cells) were detected from all 20 our samples and 20 HBCA samples (see Methods). This integrated dataset was used for subsequent analyses to enhance generalizability.

As CALB1 and CALB2 expression alone were insufficient to confirm that the ExN_CALB clusters represented PVT neurons, we validated this classification using two approaches: (1) spatial transcriptomic deconvolution and (2) anchor-based label transfer using mouse PVT neuron referenced data. We obtained spatial transcriptomic data of the medial thalamus from one control sample as a representative specimen (Fig. 2b). The PVT was characterized through unbiased clustering of spatial transcriptomic data and partially corresponded to the areas immunostained by calretinin, a protein coded by CALB236 (Fig. 2c). Deconvolution of thalamic spatial transcriptomic data using the integrated cell clusters by CytoSPACE39 indicated that ExN_CALB cells were predominantly localized in the paraventricular regions of the thalamus with a clear demarcation from other excitatory neuron populations (Fig. 2d, e), aligning with the anatomical definition of the PVT by Uroz et al.36. Oligodendrocytes were concentrated in the lateral side of the thalamus, while astrocytes were primarily found in the medial side. Oligodendrocyte precursor cells (OPCs) and microglia displayed a diffuse distribution throughout the thalamus. Expectedly, ExN_CALB corresponded most closely with mouse PVT neurons when employing anchor transfer with published single-cell datasets from mice23,40 (Fig. 2f). These data confirmed the validity of ExN_CALB as PVT neurons. Upon closer examination of PVT neuron subtypes, the majority of human PVT neurons, particularly those from our dataset, corresponded to mouse posterior PVT neurons (Fig. 2g).

Compared to other thalamic excitatory neurons, PVT neurons exhibited a distinctive expression profile for genes (Supplementary Data 5) associated with dopamine signaling (e.g., MOXD1, SNCA, PRKN) and the regulation of calcium and ion channels (e.g., MCUB, CACNG3, PDE4B, SCN9A), as determined by GSEA (Fig. 2h, Supplementary Data 2). These PVT neuron markers were particularly expressed in the PVT regions identified through spatial transcriptomics (Fig. 2i). Notably, DRD2, highly expressed in the medial thalamus relative to the frontal cortex, was prominent in the PVT compared to other thalamic areas. This observation corroborates the enrichment of dopamine-signaling-related genes in PVT neurons and aligns with Rieck et al.’s report of high D2 receptor density within paraventricular thalamic regions41. The higher expression of SCNA in PVT is consistent with Schulmann et al.27. Collectively, these findings define PVT neurons as a subset of excitatory neurons uniquely characterized by the expression of calcium-binding protein genes, with a specific emphasis on dopamine signaling and calcium ion regulation.

Cellular compositional changes in the medial thalamus of BD
We subsequently examined the cellular compositional and transcriptional alterations in the medial thalamus and frontal cortex between 21 patients with BD and 20 controls. The thalamus and frontal cortex were analyzed separately to enhance precision. Cellular compositions refer to the relative proportions of cell clusters within the total population of cell nuclei.

We obtained 146,643 cell nuclei in 18 unique cell clusters and 2139 droplets in two NeuCyto clusters after QC (Fig. 3a and Supplementary Fig. 3a). Two NeuCyto clusters were removed as non-nuclei from the subsequent analysis. Tha_ExN_CALB in Fig. 3a corresponded to PVT neurons. All major cell classes were consistently detected across the 41 samples (Supplementary Fig. 3b). The Tha_ExN_CALB cluster was also robustly identified in all samples among thalamic excitatory neurons with two outliers, BD13 and CT10 (Supplementary Fig. 3c). Significant compositional reductions were observed exclusively in thalamic excitatory neurons (FDR<0.05, FC=0.54–0.59), including PVT neurons (Tha_ExN_CALB, FDR=0.00488, FC=0.55), as determined by sccomp42 (Fig. 3b, Supplementary Data 6). In contrast, five oligodendrocytes and astrocytes clusters showed compositional increases (FDR<0.05, FC=1.22–1.55). The reduction of Tha_ExN_CALB remained significant even after excluding the two outlier samples, BD13 and CT10 (FDR=0.011, FC=0.592). In examining six major cell classes, we observed a significant compositional reduction in excitatory neurons (FDR=0.00525, FC=0.58) and compositional increases in oligodendrocytes and astrocytes in BD (FDR=0.0168 and 0.004, FC=1.29 and 1.44, respectively) (Supplementary Fig. 3d). Complementary analysis using scCODA43 also supported a consistent compositional reduction in thalamic excitatory neurons, including PVT neurons (Supplementary Data 7). In the scCODA analysis, Tha_ExN_CALB and Tha_ExN_PVALB_1 exhibited comparable compositional decreases (FC=0.686 and 0.663).

Variability in tissue dissection may also influence cell-type composition; however, such variability is difficult to control during sampling and to explicitly model statistically. To address this potential confound, we randomly selected 10 BD and 10 control samples from the full dataset (21 BD vs. 20 controls) and used sccomp to evaluate BD-associated FCs across 1000 random subsamples. Notably, over 97.5% of these subsamples yielded diagnosis-associated FCs < 0.79 for both Tha_ExN_CALB and the broader thalamic excitatory neuron population (Supplementary Fig. 3e), suggesting that the observed compositional reductions are unlikely to be driven by a small number of anatomically atypical outlier samples.

To validate this finding, we compared the cell density of thalamic excitatory neurons and oligodendrocytes by immunohistochemical staining of VGLUT2 and SOX10, calculating cell density from 548,423 excitatory neurons and 1,819,912 oligodendrocytes in 36 thalamic samples (Supplementary Data 8). This immunohistochemical assessment revealed significant reductions of the thalamic excitatory neurons in BD (P=0.00781, FC=0.53) and no significant differences in the oligodendrocytes (P=0.252, FC=0.91) (Fig. 3c, d). If incorporating PMI and/or pH to the GLM models, the results remained consistent (Supplementary Data 9). Note that immunohistochemistry directly assessed cell density, while compositional analysis using snRNA-seq data assessed relative proportional changes. The immunohistochemical assessment indicated that the compositional increase in oligodendrocytes and astrocytes in BD should be mainly derived from prominent reductions of excitatory neurons, not reflecting absolute cell density increases. If we set Tha_Oligo_1 as a reference cluster in scCODA analysis, assuming that the absolute counts of oligodendrocytes were comparable between BD and control from IHC results, we consistently observed compositional decreases in thalamic excitatory neurons (Supplementary Data 7).

Differentially expressed genes in the medial thalamus of BD
Subsequently, we identified differentially expressed genes (DEGs) within each cluster using MAST44, and categorized these DEGs into “probable” (robustness score > 0.5) and “possible” (robustness score > 0) DEGs, according to the confidence level of detection. To evaluate the functional implications, we conducted gene ontology (GO) enrichment analysis for the DEGs, defining GOs with robustness scores above 0.5 as “robust DEG-GOs” for each cluster (see Methods).

This analysis revealed 56 downregulated and 69 upregulated probable DEGs in the medial thalamus (Fig. 3e, Supplementary Data 10). Tha_ExN_CALB (PVT neurons) displayed the highest DEG count, with 22 downregulated and 43 upregulated probable DEGs (64 downregulated and 157 upregulated possible DEGs). Tha_ExN_PVALB_1 exhibited the second-largest number of upregulated DEGs, with 81.4% of these genes overlapping with those upregulated in Tha_ExN_CALB. There were 29 downregulated and 67 upregulated robust DEG-GOs in the thalamus, detected only in Tha_ExN_CALB, Tha_ExN_PVALB_1, and Tha_Micro_1 (Fig. 3f, Supplementary Data 11). Tha_ExN_CALB also demonstrated the highest number of DEG-GOs, 24 downregulated and 58 upregulated. The robust DEG-GOs were predominantly enriched in neuronal and synaptic processes, and several DEG-GOs demonstrated overlap between downregulated and upregulated categories in Tha_ExN_CALB, such as the regulation of trans-synaptic signaling (GO:0099177) and modulation of chemical synaptic transmission (GO:0050804). Given the enrichment of synaptic terms in the DEG-GOs, we also performed SynGO enrichment analysis45. The downregulated DEGs in the Tha_ExN_CALB cluster again showed the strongest SynGO enrichment, including three terms related to synaptic structure: integral component of postsynaptic membrane (GO:0099055), integral component of presynaptic membrane (GO:0099056), and integral component of postsynaptic density membrane (GO:0099061) (Supplementary Data 11).

Complementary GSEA based on fold changes supported the enrichment of neuronal and synaptic processes among downregulated genes in Tha_ExN_CALB (Supplementary Data 12). GSEA identified potassium channel activity, including KCNQ3, KCND2, and KCNC2 as the top GO term for the downregulated genes in Tha_ExN_CALB. The overlap between DEG-GO and GSEA results (GSEA-GO) was more pronounced for downregulated than upregulated genes (9/24 vs. 6/58; two-sided Fisher’s exact test, P=0.00931, odds ratio=5.07) (Supplementary Data 11). Notably, the top three DEG-GO terms among downregulated genes—transporter complex (GO:1990351), monatomic ion channel complex (GO:0034702), and main axon (GO:0044304)—were also identified by GSEA. The downregulated genes in Tha_ExN_CALB exhibit greater biological coherence and interpretability across enrichment methods. Supplementary Data 13 summarizes the downregulated probable DEGs indicated by both DEG-GO and GSEA as biologically interpretable genes for BD. Among them, KCNQ3, CACNA1C, SHISA9, KCNC2, LRRC7, ROBO2, and CNTNAP4 were consistently enriched across all three analytical approaches (DEG-GO, GSEA, and SynGO).

Upon closer inspection of probable DEGs, MAGMA-based BD genomic risk1,46 was most enriched in the downregulated DEGs in Tha_ExN_CALB (P=0.00123). Additional fine-mapping by PolyFun + SuSiE/FINEMAP47–49 further supported CACNA1C and SHISA9 as putative BD risk genes within the downregulated set. Especially, SHISA9 is one of the 36 credible genes for BD in the largest GWAS to date1. While CACNA1C is not listed in these 36 credible genes, CACNA1C has consistently ranked among the top genes by P-value in prior PGC1-4 studies1,2,50,51, and was also implicated by our FUSION-based TWAS integration. In contrast, the upregulated probable DEGs in Tha_ExN_CALB exhibited weaker enrichment for MAGMA-based BD genomic risk (P=0.0916), including no putative BD risk genes despite the greater number of upregulated genes (n=22 for downregulated vs. 43 for upregulated, Supplementary Data 10).

Transcriptional disturbances of PVT neurons in BD
The notable enrichment of DEGs and DEG-GOs in Tha_ExN_CALB (PVT neurons) led us to focus on the DEGs in this cluster (Fig. 4a). When assigning the DEGs to co-expression network modules identified through hdWGCNA52 (Supplementary Data 10), the downregulated DEGs in Tha_ExN_CALB were significantly enriched in Module-4 of excitatory neurons (Excit-M4), while the upregulated DEGs were primarily enriched in Excit-M11 and Excit-M9, with adjusted p-values < 0.05 and FCs > 2.5 for both probable and possible DEGs (Fig. 4b and Supplementary Fig. 4a). The downregulated genes exhibited strong aggregation within a single module (>67% in Excit-M4), whereas the upregulated genes displayed more heterogeneous module enrichment. The clear delineation of modules between downregulated and upregulated genes suggests distinct molecular pathways.

Among the 1540 Excit-M4 genes in Tha_ExN_CALB, the module eigengene-based connectivity measure (defined as kME in the WGCNA package53) for downregulated DEGs were significantly higher than those for other Excit-M4 genes (Supplementary Fig. 4b), indicating the DEGs’ centrality within the Excit-M4 co-expression network. Similarly, the kME for upregulated DEGs was elevated among the 236 genes in Excit-M11 and the 162 genes in Excit-M9. In addition to higher kME, the downregulated DEGs in Excit-M4 were significantly associated with core genes from the top 10 GSEA-derived GOs for this module (P=1.54×10⁻⁵ and 8.59×10⁻⁵ for probable and possible DEGs, respectively, Supplementary Data 14). The upregulated DEGs in Excit-M11 also showed enrichment in core genes from the top 10 GSEA GOs for that module (P=2.32×10⁻⁴ and 1.43×10⁻⁷ for probable and possible DEGs). We did not observe similar enrichment in Excit-M9. Therefore, the downregulated DEGs in Excit-M4 and the upregulated DEGs in Excit-M11 are likely core genes within their respective modules that could exert influence on, or be influenced by, a broad array of genes within the same co-expression network. Excit-M4 genes are enriched in GOs related to axonal function and calcium ion regulation, while Excit-M11 genes are enriched in GOs related to calcium ion dynamics and synaptic activity (Supplementary Data 14). This consistency of module-based GO and DEG-GO in PVT neurons further supported the validity of the DEG-GOs in these neurons.

Molecular insight into PVT neuron pathology
In addition to GO and module enrichment analyses, we examined protein-protein interaction (PPI) networks of DEGs for molecular insights into BD pathology. The downregulated and upregulated probable DEGs within PVT neurons (Tha_ExN_CALB) exhibited significant PPI network enrichment (P=3.43×10⁻⁹ and 1.46×10⁻¹², respectively). When combined, these DEGs formed cohesive PPI networks (PPI enrichment P<1.0×10⁻¹⁶, the lower limit of STRING54, Fig. 4c). The cohesive PPI networks of downregulated and upregulated DEGs indicated their tight interactions at protein levels, despite the distinctive co-expression modules.

In PVT neurons, the downregulated genes showed consistent GO term enrichment in both GO analysis and GSEA, and the MAGMA-based BD genomic risk was relatively prominent for downregulated genes. Thus, we focused on the downregulated genes. When extending the PPI analysis to possible DEGs, the downregulated possible DEGs in PVT neurons also demonstrated robust PPI networks (enrichment P<1.0×10⁻¹⁶), encompassing genes within the top three DEG-GOs as biological processes (BP) or cellular components (CC) (Fig. 4d). These key genes, including CACNA1C and SHISA9 as putative BD risk genes, were identified as core DEGs based on their PPI and GO enrichment. Particularly, KCNQ3, GRIN2A, and GRM5 intersected with the top three BP GOs—modulation of chemical synaptic transmission (GO:0050804), regulation of trans-synaptic signaling (GO:0099177), and regulation of monatomic ion transmembrane transport (GO:0034765). KCNQ3 and KCNC2 intersected with the top three CC GOs—transporter complexes (GO:1990351), monatomic ion channel complexes (GO:0034702), and main axon (GO:0044304). These synaptic and ion-channel-related genes as PPI network hubs aligned with the biological pathways indicated by BD genomic studies1,2,4 and our DEG-GO/GSEA/SynGO results. Particularly, the overlap of KCNQ3 in both BP and CC top GOs underscores its significant relevance to BD, besides putative BD risk genes.

Disturbed interaction between PVT neurons and microglia in BD
In addition to Tha_ExN_CALB (PVT neurons), we identified downregulated robust DEG-GOs within Tha_Micro_1 (microglia). Notably, these two clusters were the only ones displaying downregulated robust DEG-GOs (Fig. 3e, f). Furthermore, all five DEG-GOs identified within Tha_Micro_1 were related to synaptic function (Supplementary Data 11). When we assigned the DEGs in Tha_Micro_1 to co-expression network modules, six of the ten downregulated probable DEGs (60%) and 32 of the 55 downregulated possible DEGs (58.2%) were assigned to Module-11 of microglia (Micro-M11) (P=8.5×10⁻⁴ and 1.20×10⁻¹⁴, FC=4.36 and 4.22, hypergeometric test) (Supplementary Data 10). Among the probable DEGs, SYNDIG1, CYFIP1, and TIAM1 emerged as promising synaptic DEGs (Fig. 5a), each exhibiting robustness scores exceeding 0.95 and included within the five synaptic DEG-GOs. Additional GSEA identified synaptic enrichment among the downregulated DEGs in Tha_Micro_1 (Supplementary Data 12), again highlighting SYNDIG1 and CYFIP1. SynGO enrichment analysis suggested potential enrichment of GO:1905274 (regulation of modification of the postsynaptic actin cytoskeleton) (Supplementary Data 11).

The simultaneous downregulation of synaptic genes in both microglia and PVT neurons implicates a disturbance in microglia-neuron interactions for synapses in BD. To investigate this possibility, we inferred potential cell-to-cell interactions as trans-cellular meta-gene programs using non-negative matrix factorization (NMF) by GeneNMF55. While each NMF-based meta-gene program was generally confined to one major cell class (Supplementary Fig. 5a), we extracted one meta-gene program encompassing excitatory neurons and microglia (Excit-Micro meta-gene program, Fig. 5b), likely reflecting their physiological interactions. This meta-gene program (MP7), mainly driven by PVT neurons (Tha_ExN_CALB), was the only one showing significantly lower scores in BD (Bonferroni-adjusted P=0.0217 by generalized linear mixed model [GLMM], Fig. 5c). This reinforced the possibility that interactions between excitatory neurons and microglia are disrupted in BD. The Excit-Micro meta-gene program with downregulation in BD was robustly detected in another down-sampling condition (Supplementary Fig. 5b), supporting the robustness of our method. Supplementary Data 13 summarizes the downregulated genes in Tha_ExN_CALB and Tha_Micro_1 that were enriched in one or more GO terms in both DEG-GO and GSEA analyses. Notably, these genes belonged to the same co-expression modules in each cluster.

At the gene level, SYNDIG1 regulates excitatory synapse number and SHISA9 encode a regulator of AMPA receptor (AMPAR) complex56,57, with both genes’ robustness scores exceeding 0.7 (Supplementary Data 10). Notably, SYNDIG1 was prominently expressed in microglia, while SHISA9 was primarily expressed in excitatory neurons (Fig. 5d). SHISA9 consisted of core PPI networks of BD-associated downregulated genes in PVT neurons (Fig. 4d) and overlapped in the downregulated genes in Tha_ExN_PVALB_1, as one of the 36 putative BD-causal genes1. Beyond microglia and excitatory neurons, thalamic T cells were included in this meta-gene program. Figure 5e illustrates potential ligand-receptor interactions between neurons and microglia/T cells, as analyzed through LIANA58. Although we did not detect downregulated DEGs as ligand–receptor pairs, KCND2 in Tha_ExN_CALB, LAMB1 in Tha_ExN_PVALB_2, and KCNQ1 in Tha_Micro_1 were identified as downregulated DEGs among the LIANA-based ligand-receptor pairs. These genes encode ion channels or extracellular matrix proteins that may potentially contribute to Excit–Micro interactions. However, as this analysis is an inference, the precise ligand–receptor mechanisms remain to be elucidated in future studies.

Cellular compositional changes in the frontal cortex of BD
We performed a detailed analysis of frontal cortex samples to characterize compositional and transcriptional alterations associated with BD relative to those observed in the thalamic samples. We obtained 165,585 cell nuclei in 22 unique cell clusters and 2369 droplets in one NeuCyto cluster from the same cohort of 21 patients with BD and 20 controls (Fig. 6a and Supplementary Fig. 6a, b). NeuCyto cluster was removed from the subsequent analysis.

We observed a compositional decrease in FrC_ExN_L6ITCar3 (FDR=0.00463, FC=0.65) and compositional increases in FrC_Astro (FDR=0.0416, FC=1.36) and FrC_Oligo (FDR=0.00450, FC=1.41) by sccomp42 (Fig. 6b, Supplementary Data 6). Within the major cell class analysis, we observed compositional reductions in excitatory and inhibitory neurons (FDR=0.0283 and 0.0218, FC=0.848 and 0.863, respectively) and compositional increases in oligodendrocytes (FDR=0.0488, FC=1.39) (Supplementary Fig. 6c). The scCODA analysis was not consistent to the sccomp analysis, exhibiting the compositional decreases of FrC_ExN_L2/3IT_1 (FC=0.797), FrC_ExN_L5IT_4 (FC=0.866), and FrC_InN_CGE_LAMP5 (FC=0.848). In contrast to findings in the thalamic excitatory neurons, the cortical neurons showed inconsistent results and smaller reductions (Fig. 6b vs. Fig. 3b). Note that the absolute cell count in situ would be different from these compositional (proportional) changes. Given the immunohistochemical assessment for the thalamus, compositional increases in oligodendrocytes might be derived from more prominent reductions of neurons, not reflecting the absolute increases of oligodendrocytes.

Differentially expressed genes in the frontal cortex of BD
Our DEG analysis revealed 55 downregulated and 73 upregulated probable DEGs in the frontal cortex (Fig. 6c, Supplementary Data 10). FrC_ExN_L5IT_4 and FrC_ExN_L6Car3 exhibited the highest numbers of probable DEGs (Supplementary Data 10) and 11 robust DEG-GOs (Supplementary Data 11). While the DEG-GOs primarily pertained to synaptic and neuronal functions, their number in the frontal cortex was lower than that observed in the thalamus (Fig. 6d vs. Fig. 3d). Among cortical DEG-GO, only GO:0045211 (postsynaptic membrane), an upregulated GO term in FrC_ExN_L6ITCar3, was identified by FC-based GSEA, and no SynGO enrichment was observed. These findings suggest that cortical DEGs exhibit lower biological interpretability compared to thalamic DEGs.

The enrichment of DEGs in FrC_ExN_L5IT_4 and FrC_ExN_L6Car3 prompted us to focus on the DEGs in these clusters (Fig. 7a). When assigning the DEGs to co-expression network modules identified through hdWGCNA, we observed specific module enrichments in FrC_ExN_L5IT_4 and FrC_ExN_L6Car3 (adjusted P<0.05 and FC>2.5 for both probable and possible DEGs). However, the enrichment for FrC_ExN_L6Car3 was less pronounced (Supplementary Fig. 7a, b). In FrC_ExN_L5IT_4, the upregulated DEGs were primarily enriched in modules ExcitDL-M3 and ExcitDL-M19 (ExcitDL: deep-layer excitatory neuron, Fig. 7b). The module scores of DEGs within FrC_ExN_L5IT_4 exceeded those of other genes in their respective modules, indicating more central roles of these DEGs within their co-expression networks (Supplementary Fig. 7c, d). In contrast to the findings in the thalamus, FrC_ExN_L5IT_4 and FrC_ExN_L6Car3 exhibited weaker protein-protein interaction (PPI) network enrichment, lacking a clear distinction in interactions between downregulated and upregulated DEGs (Fig. 7c vs. Fig. 4c). Among the probable DEGs in these two clusters, only upregulated DEGs in FrC_ExN_L6Car3 were nominally enriched in MAGMA-based BD genomic risk (P=0.00454), including MACROD2 as a putative BD risk gene.

Comparison of the thalamus and cortex
As described above, the medial thalamus exhibited more pronounced compositional and transcriptional alterations than the frontal cortex. Specifically, PVT neurons showed the most significant compositional and transcriptional shifts among the various cell types across both regions. Furthermore, downregulated genes in PVT neurons showed the most prominent enrichment of MAGMA-based BD genomic risk among all the probable DEG sets in the thalamus and cortex analysis (Fig. 8a) as the only one category surviving Bonferroni correction (adjusted P=0.0295, n of tests=14). This enrichment pattern extended beyond BD, as downregulated DEGs in PVT neurons also showed the highest enrichment for MAGMA-based genomic risk associated with schizophrenia33 and major depressive disorder59, while showing weaker or nonsignificant enrichment for other psychiatric conditions: attention-deficit/hyperactivity disorder, post-traumatic stress disorder, obsessive-compulsive disorder, substance use disorder, eating disorder, panic disorder, and autism spectrum disorder (Fig. 8a, details in Supplementary Data 15). These enrichment comparisons should be interpreted with caution, as differences in statistical power across the original GWAS datasets may influence the observed patterns.

We examined which cell types most expressed putative BD risk genes. Thalamic excitatory neurons, including PVT neurons, showed the highest expression of the putative BD risk genes concerning the number of expressed genes and their average expression rate (Fig. 8b). For example, CACNA1C and SHISA9 exhibited relatively higher expression in thalamic excitatory neurons than in cortical excitatory neurons (Fig. 8c). However, the differences were modest, and both genes showed the highest expression in Tha_ExN_PVALB_1, suggesting that additional factors likely contribute to the transcriptional alterations observed in PVT neurons, as discussed later.

We examined the expression patterns of BD-associated genes identified through rare variant analyses, specifically HECTD2 and AKAP113,60. Neither gene was differentially expressed in any cell type clusters in the thalamus or cortex at an FDR threshold of <0.1. The expression levels of HECTD2 and AKAP11 were not notably elevated in thalamic excitatory neurons, and their spatial expression patterns showed no specificity to the PVT region (Supplementary Fig. 8a, b). These findings may reflect distinct mechanisms underlying BD-associated genes identified through rare variants compared to those implicated by common variant analyses.

In parallel with region-specific variations, significant parallels were identified between the thalamus and cortex. First, a recurring trend of compositional neuronal reduction coupled with compositional increases in oligodendrocytes and astrocytes were observed (Figs. 3b, 6b). These findings suggest potential brain-wide neuronal damage in BD. Second, probable DEGs identified in PVT neurons (Tha_ExN_CALB) overlapped with those in deep-layer cortical excitatory neurons (FrC_ExN_L5IT_4 and FrC_ExN_L6Car3), exhibiting similar patterns of upregulation and downregulation (Fig. 8d). Notably, KCNQ3, SLC8A3, and ATXN1 emerged among the commonly downregulated genes. They were highlighted as key candidates in the DEG-GO analyses of both PVT neurons and cortical deep-layer excitatory neurons, underscoring their potential as markers of widespread pathology. KCNQ3 encodes a subunit of voltage-gated potassium channels integral to neuronal excitability regulation. KCNQ3 was a core DEG in PVT neurons, and its broad downregulation extended to the cortex was consistent with Kaminsky et al.’s report of KCNQ3 downregulation in bulk brain tissues from patients with BD61. SLC8A3, a transmembrane sodium/calcium exchanger, is essential for maintaining Ca²⁺ homeostasis in both neuronal and cardiac cells, while ATXN1, a DNA-binding protein, is associated with spinocerebellar ataxia, a neurodegenerative disease. Although these genes were downregulated across the thalamus and cortex, their functional impact on PVT neurons might be more pronounced, given the unique characteristics of PVT neurons identified in previous analyses (see Fig. 2f).

Third, Excit-Micro meta-gene programs were also identified in the cortex (Fig. 8e). Besides, SYNDIG1, which exhibited the largest downregulation in the thalamic microglia, also exhibited the most notable downregulation in cortical microglia (Fig. 8f). Despite a low robustness score and a lack of significant downregulation of cortical Excit-Micro meta-gene programs in BD (Fig. 8e), the observed downregulation of SYNDIG1 in cortical microglia mirrors findings in thalamic microglia, suggesting possible broad disturbances of Excit-Micro meta-gene programs in BD. This hypothesis warrants further validation with more cell nuclei counts and additional samples. We applied GeneNMF on the integrated dataset encompassing both thalamic and cortical regions to examine trans-regional gene programs. While this analysis did not reveal any discernible trans-regional meta-gene programs, the thalamic and cortical Excit-Micro meta-gene programs were still identified in the thalamus-cortex integrated dataset, further validating the robustness of our approach (Supplementary Fig. 8c).

Although we did not identify trans-regional cortico-thalamic meta-gene programs, we extended the LIANA analysis to explore potential interactions between the cortex and thalamus, focusing specifically on connections between cortical deep-layer excitatory neurons (ExN_FrC_DL) and thalamic excitatory neurons (ExN_Tha) (Fig. 1b). Among the top 20 ligand-receptor pairs in the cortex-to-thalamus direction, eight involved DEGs in either ExN_Tha or ExN_FrC_DL (Supplementary Fig. 9), and four DEG-containing pairs were detected in the thalamus-to-cortex direction. These results suggest potential disruptions in cortico-thalamic connectivity mediated by altered ligand-receptor interactions. However, this analysis is exploratory and hypothesis-generating; further studies are required to more definitively characterize these connections.

Medication signature
Multiple medications and substances were detected in most postmortem brain samples at the time of death (Supplementary Fig. 10a). As subgrouping based on substances at death or medication history was not feasible in our dataset, we assessed the potential overlap between medication signatures (previously reported DEGs induced by specific medications/substances) and the DEGs identified in our study. In the four neuronal clusters with the largest numbers of DEGs, there was generally no substantial overlap with medication signatures for alcohol, antipsychotics, antidepressants, lithium, or valproic acid, except in a few instances (Fig. 8g, Supplementary Data 16). The upregulated genes in Tha_ExN_CALB significantly overlapped with genes upregulated by olanzapine in medium spiny neurons (MSNs)62 (Bonferroni-adjusted P=0.00355, n of tests=18 per cell type). A similar, though weaker, overlap was observed in Tha_ExN_PVALB_1 (adjusted P=0.0772). Another finding was the overlap between downregulated genes in FrC_ExN_L6ITCar3 and genes downregulated by valproic acid in forebrain organoids63 (adjusted P=0.0573). Although MSNs represent a distinct neuronal subtype from thalamic excitatory neurons, the consistent overlap between Tha_ExN_CALB/Tha_ExN_PVALB_1 and olanzapine-treated MSNs for upregulated genes suggests that these transcriptional changes could be influenced by medication exposure. However, no notable overlaps were observed for other combinations.

Importantly, the downregulated genes in Tha_ExN_CALB—which were a central focus of our results—were unlikely to be attributable to medication or substance effects in this analysis.

Single-molecule fluorescence in situ hybridization for SYNDIG1 in microglia
The unexpected detection of SYNDIG1 in thalamic and cortical microglia prompted further validation of its expression and spatial distribution using RNAscope. SYNDIG1 transcripts were detected in CX3CR1-positive cells, consistent with a microglial identity (Fig. 9a). CX3CR1/SYNDIG1 double-positive cells—representing SYNDIG1-expressing microglia—were preferentially localized to the gray matter of the frontal cortex (Fig. 9b). In the thalamus, these cells appeared more abundant in medial regions, although this pattern was not definitive. Comparison of representative thalamic male samples from five BD cases and five controls showed relatively lower SYNDIG1 expression in BD microglia (raw FC=0.801, Fig. 9c). While generalized linear mixed model (GLMM) analysis with a negative binomial distribution did not reach statistical significance, it indicated a trend toward reduction (GLMM P=0.153, estimated FC=0.783).

Methods
Post-mortem brain samples
Postmortem brain specimens were sourced from the Douglas–Bell Canada Brain Bank (www.douglasbrainbank.ca). Frozen samples were obtained from 21 patients with BD and 21 matched controls (Supplementary Data 1), with written informed consent provided by the next of kin. No compensation was provided to donor families. Due to the limited availability of postmortem thalamic tissue, no formal sample size calculation was performed; instead, we included all eligible samples available at the time of collection, ensuring sex matching between BD cases and controls. Each group comprised 14 males and 7 females prior to quality control. Sex was determined from coroner files; gender identity information was not available for deceased donors. Analyses were not stratified by sex because of the limited sample size within each diagnostic group. The average age at death was significantly higher in the control group than in the BD group (52.5 vs. 44.3 years; two-sided t test, P=0.0472), reflecting our prioritization of paraventricular thalamic region availability and sex matching over age matching. The cause of death was significantly skewed towards suicide in the BD group (two-sided Fisher’s exact test, P=2.33×10−8), consistent with the elevated suicide risk in BD and the characteristics of this brain bank. Of the 42 donors, 41 were Caucasian and one Asian. Donors exhibited substantial heterogeneity in substances present at death, precluding stratification into subgroups of more than two individuals. Summary statistics of substances present at death are presented in Supplementary Fig. 10a.

Bulk specimens of the medial thalamus and frontal cortex (BA10) were dissected from archived fresh-frozen slabs and then sectioned at a thickness of 200μm using a CM3050S cryostat-microtome (Leica Microsystems, Wetzlar, Germany). First, the whole brain was bisected into left and right hemispheres. The left hemisphere was snap-frozen in 2-methylbutane precooled to −40°C and stored at −80°C, while the right hemisphere was fixed in 10% neutral buffered formalin (phosphate buffer, pH 7.0; Fisher Scientific, Cat. #SF-100-20) for paraffin embedding and sectioning. Within a half hemisphere, the PVT was macroscopically identified, with its anterior boundary defined by the coronal plane at the anterior edge of the mammillary body and its posterior boundary by the coronal plane at the anterior edge of the posterior commissure, spanning approximately 20mm along the anterior-posterior axis. A 10-mm-thick coronal slab encompassing the PVT was then selected and excised. Each slab extended anteriorly to the center of the mammillary body and posteriorly to midbrain structures, including the red nucleus and substantia nigra. From each slab, tissue blocks were prepared to include the PVT together with major adjacent landmarks (putamen, caudate, and insular cortex), ensuring accurate anatomical orientation and precise localization of the PVT (Fig. 1a, Supplementary Fig. 10b-d). Dissections were typically performed from the posterior aspect of the slab, with each block carefully inspected to confirm that the slab contained the appropriate rostrocaudal level of the PVT. Although the procedure was designed to target the mid-PVT, unavoidable interindividual variation in gross neuroanatomy occasionally resulted in slight shifts in sampling position along the anterior-posterior axis. The medial region of the thalamic cryosections and the gray matter part of the cortical cryosections were excised for snRNA-seq. Approximately 50mg of cryosections, corresponding to two to four cryosections, underwent the subsequent cell nuclei isolation. This research was conducted in accordance with the Helsinki Declaration and received approval from the Research Ethics Committee of the Faculty of Medicine, Juntendo University (approval ID: M19-0278).

Cell nuclei isolation
Cell nuclei were isolated from the prepared specimens, with myelin sheath removal achieved through a four-layer density gradient centrifugation. The tissue was homogenized in 900µl of homogenization buffer (HB; 250mM sucrose, 25mM KCl, 5mM MgCl2, 20mM Tris-HCl pH 8.0, 1× cOmpleteTM, EDTA-free Protease Inhibitor Cocktail (Roche, Basel, Switzerland), 0.1mM dithiothreitol (Fujifilm, Tokyo, Japan) and 0.04Uµl−1 Protector RNase Inhibitor (Roche)) by using approximately 20 strokes of a pellet pestle and then by pipetting 10-times through a 1000-µl low retention fine point tip (Rainin, Oakland, USA). The homogenized tissue was collected as a pellet by centrifugation at 500×g for 5min at 4°C. The pellet was suspended in 1000µl of HB supplemented with 0.05% polyoxyethylene sorbitan monolaurate (Tween 20, Nakarai, Kyoto, Japan) by pipetting, and then incubated for 12min with gentle pipetting at times through a 1000-µl low retention fine point tip to dissociate nuclei. The dissociated nuclei were collected as a pellet by centrifugation at 500×g for 5min at 4°C, and then suspended in 900µl of nuclei suspension buffer (NSB; 1× PBS, 1% (w/v) UltraPure BSA (Thermo Fisher, Waltham, USA) and 0.04Uµl−1 Protector RNase Inhibitor). The nuclei suspension was filtered through 70-µm and 40-µm Mini Cell Strainer (Funakoshi, Tokyo, Japan), serially, and then collected as a pellet by centrifugation at 500×g for 5min at 4°C. For density-gradient centrifugation, 50% iodixanol solution was made, which consists of the same composition as HB except where sucrose was replaced with iodixanol using OptiPrep (60% iodixanol, Serumwerk Bernburg AG, Bernburg, Germany), and then mixed with HB to prepare 25%, 29%, and 35% iodixanol solution. The nuclei pellet was suspended in 1000µl of 25% iodixanol solution. A four-layer discontinuous iodixanol density gradient was assembled by sequentially overlaying 200µl of 35% iodixanol solution as a cushion, 400µl of 29% iodixanol solution, all nuclei suspension in 25% iodixanol solution, and 150µl of HB. The density-gradient was centrifuged by 2-steps: at 500×g for 6min at 4°C, and then at 5000×g for 30min at 4°C. After 1400µl of upper layers (partially including 29% iodixanol layer) was removed totally from the top surface of layers, nuclei on 35% iodixanol layer were collected as a suspension in remaining solution. For iodixanol removal, the nuclei suspension was diluted by adding 1000µl of NSB, and then centrifuged at 500×g for 30min at 4°C to collect nuclei as precipitates. The nuclei were washed with 500µl of NSB, and then re-precipitated by centrifugation at 500×g for 10min at 4°C. The nuclei were suspended in 100µl of NSB, 30µl of which was diluted 5-fold with NSB, and then filtered through a 40-µm porosity Flowmi cell strainer (Sigma Aldrich, Burlington, USA). Isolated single nuclei were counted, and if needed, more diluted with NSB to adjust the concentration of nuclei to the recommended concentration for the following library preparation according to the manufacture’s guidelines (10x Genomics). All procedures were performed on ice or at 4°C. Two specimens from one BD patient and one control were processed simultaneously from nuclei isolation to snRNA-seq library preparation to minimize systematic bias due to experimental batches. All procedures were performed by the same experimenter (KI), who was blind to the disease status of the samples, to avoid inter-experimenter differences.

Single-nucleus RNA-sequencing (snRNA-seq)
Nuclei isolated from the medial thalamus and frontal cortex underwent library preparation using the 10x Chromium 3’ RNA v3.1 dual index protocol, following the manufacturer’s instructions (10x Genomics, Pleasanton, USA). The target recovery rate was set at 5000 nuclei, with library amplification achieved through 14 PCR cycles. The snRNA-seq libraries were sequenced on the NextSeq 2000 platform (Illumina, San Diego, USA) in accordance with the protocols provided by Illumina and 10x Genomics. To minimize experimental bias, 10x Genomics library preparation was performed exclusively by the same experimenter (KI), with all case–control pairs processed in parallel. Due to the limitations of the NextSeq 2000 platform, which could accommodate only a few single-nucleus RNA-seq libraries per run, the dataset spans multiple sequencing batches, including several salvage runs needed to achieve sufficient read count per nuclei.

The output files were processed with Cell Ranger (v7.0.0, 10x Genomics) using default parameters (reference: GRCh38-A) to generate count data. Sample BD21 thalamus was excluded due to a low recovery rate (<1000 nuclei, target of 5000). We did not process the BD21 frontal cortex to maintain consistency in paired analysis. The final dataset comprised 82 snRNA-seq count datasets from 41 individuals (21 BD cases and 20 controls) across the thalamus and frontal cortex. All samples achieved read depth per nuclei exceeding 50,000, with a median depth of 74,000 reads per nucleus.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Nishioka et al. submit an interesting manuscript that investigates postmortem brain samples from bipolar disorder (BD) patients using single-nucleus RNA-seq. The authors analyzed a cohort of >80 samples of the prefrontal cortex and thalamus from ~40 subjects to identify new insights into brain dysfunction in BD. The authors showed that gene expression and cell composition changes are predominant in the thalamus, and much lower in the cortex. This led them to focus most of their analyses on the paraventricular thalamic nucleus. Subsequent analyses showed that neurons of the paraventricular thalamic nucleus had large decrease in cell number in BD and exhibited the largest changes in gene expression. In particular, genes involved in synaptic processes were downregulated and few dysregulated hub genes were proposed. Using tools predicting cell-cell interaction changes, the authors proposed that dysfunction of microglia might mediate impairment of synaptic function in BD paraventricular thalamic nucleus neurons. Overall, the study is novel for the field and provides good insights into specialization of thalamic excitatory neurons and their dysfunction in BD. Nevertheless, few major issues, mainly concerning data analysis, should be addressed before accepting the publication.\nMajor issues:\n1. Covariate control needs significant improvement.\na. There have been multiple batches of sequencing, however, it is not described how many batches, how many samples per batch, etc. Given that sequencing batch effect is one of the major confounders in snRNA-seq data, the authors should provide thorough analysis of sequencing batch impact on their clustering, cell type annotation, differential gene expression and cell composition. Supplementary Table 1 should also contain information regarding sequencing batch.\nb. Medication is another large confounder in such analyses, since BD brains are usually exposed to a variety of drugs for a long period of life vs neurotypical brains. Potential impact of medication should be addressed, such as tests whether different medications within their BD cohort have impact on gene expression, e.g. by comparing subgroups of patient samples, comparing to other snRNA-seq studies, or to medication signatures identified in previous animal/patient studies.\nc. Some samples that were analyzed have high PMI, >24hours, and there are many studies showing the impact of high PMI on the mRNA content, including benchmarking studies for RNA-seq experiments. The authors should perform a thorough analysis for potential confounding impact of PMI and if necessary exclude samples that are too different from low PMI samples and interfere with data analysis.\nd. Additionally, the authors should statistically analyze potential impact of all metadata variables, including age, sex, drug exposure etc.\n2. Although authors did include filter of ambient RNA, it should be improved.\na. There is some ambiguity with the analysis of detection and removal of ambient RNA. In Methods Line 565: DecontX was chosen over SoupX and no ambient RNA decontamination\nDo the authors refer to “no ambient RNA decontamination” as a tool or package? The authors should be more specific.\nb. Overall, given potential contamination with ambient RNA, more robust analysis is necessary. We propose to use CellBender and provide the following data in Supplementary figures:\n- Sample-wise expression of ambient RNAs\n- CellBender ELBO (test and train) plots\n- CellBender cell calling curves displaying cell probabilities of droplets\n3. Data integration\nThe authors used iterative integration with Seurat and Harmony. There is no justification provided for why Seurat-only or Harmony-only integration are not enough. Simple comparison plots will help the reader to understand the rationale for choosing iterative integration.\n4. Differential gene expression\na. We could not find information whether p values for differentially expressed genes were adjusted. The authors should outline how they did the adjustment.\nb. While GO analysis of DE genes is often used in snRNA-seq studies, it works better for larger gene lists than the authors obtained in this study. GSEA is better applicable here, since GSEA implements all genes, not only short list of significantly changed DE genes.\n5. Integration with genetics is rather preliminary and more sophisticated tools can be implemented for integration with the latest GWAS BD study. This will help validating the findings in the study for their relevance to BD. The thalamus and cortex should be analyzed separately by integration with GWAS data to show potential regional specialization with respect to BD-linked brain dysfunction.\nFinally, at the moment, it seems that most association of DE genes in this study with BD risk factors is based only on 1 gene from 36 credible genes from the last GWAS study - SHISA9. We would caution the authors for not to overinterpret such association.\n6. Compositional analysis\nWhile it is justified to use sccomp package for compositional analysis, the package does not control for inter-dependency of proportional changes, and the authors admit it when comparing to immunohistochemical assessment, Lines 169-172:\nNote that immunohistochemistry directly assessed cell density, while compositional analysis using snRNA-seq data assessed relative proportional changes. The immunohistochemical assessment indicated that the compositional increase in oligodendrocytes and astrocytes in BD should be mainly derived from prominent reductions of excitatory neurons, not reflecting absolute cell density increases.\nTo avoid such issue, there are packages that allow for independent compositional analysis of cell proportions, such as scCODA. The authors should use one of such packages to cross-compare snRNA-seq with immunohistochemistry.\n7. Cell-cell interaction\na. Based on the text lines 270-283, it sounds as the authors propose that SHISA9 protein on neurons interacts with SYNDIG1 on microglia? If yes, then some literature search should help confirming this potential interaction. Otherwise, we suggest concentrating interpretation on general discussion on how microglia modulate synaptic processes in excitatory neurons via AMPA receptor complex. Additionally, LIANA analysis should help with such discussion, however the data is just shown in Figure 5e, without any interpretation.\nb. One missing opportunity in the study is analysis of changes in cortico-thalamic connectivity in BD. This can be done using similar cell-cell interaction tools.\n8. Such snRNA-seq studies usually require validation of some findings by single-molecular FISH. We propose selecting the most interesting finding in gene expression and validate it by RNAscope or similar technologies.\n9. The cortex does not show large changes in gene expression that is why the authors argue to focus on the thalamus. Nevertheless, this is surprising given previous role of the cortex that has been proposed in brain dysfunction in BD. There is a lot of focus in the Discussion section on the thalamus, but such surprisingly minor changes in the cortex are not discussed. We think these findings are very important, since the study can pave the way for future analyses of postmortem brain tissue in psychiatric disorders. Currently, almost all studies focus on the prefrontal cortex and maybe in the future the thalamus should be the major region for such snRNA-seq analyses.\nMinor\n1. Discussion lacks similarity of BD and schizophrenia, although CACNA1C is one of the major genes in the current study and it has been identified as genetic risk factor by schizophrenia GWAS and not BD.\n2. In general, although Discussion is quite long, some points are extensively discussed, whereas others are missing. For instance, close relation of thalamic excitatory neurons with GABAergic neurons is interesting. Is it related to cell fate or neuronal specialization? Similarly, cardiac tissue similarity can be moved from Results to Discussion section.\n3. Lines 137-141 compare data from the current study to previously published organoid data. This comparison is completely separated from all other analyses in the study and in the current state should be removed from the manuscript.\n4. In Figure 3b – it will help to add direction of \"increased/decreased\" for BD.\n5. We are not sure why the authors decided to include reference to NeuCyto in Lines 158-159. If this “cell type” is ambient RNA contamination, then it is difficult to make conclusions based on such “cell type”.\n6. Please correct Lines 255-256:\nrelated to synaptic function, although microglia are primarily regarded as an immune-related cell type\nMicroglia has well known role in synaptic plasticity and modification of neuronal connectivity, which is also one of their primary roles in the brain.\n7. Figure 8a: show p-adjusted line on the figure\n8. Lines 452-460 repeat lines 463-472.\n9. Correct misprint:\nFig1a Neruon-microglia synaptic interaction in BD (Figure 5)", "reviewer_2": "Nishioka et al. present a snRNA-seq study of PVT in 21 bipolar disorder (BD) patients and 20 controls. BD is a severe mental health condition marked by recurring manic and depressive episodes, requiring new treatments due to resistance, adverse effects, and high suicide rates. Recent research has identified the thalamus, particularly the paraventricular thalamic nucleus (PVT), as a crucial brain region for mood regulation, with significant differences in BD patients compared to controls. Authors report that thalamus samples from the BD patients they studied exhibited a 50% reduction in PVT neurons, which also showed notable downregulation in genes related to synaptic function and ion channel regulation. Additionally, disruptions in interactions between PVT neurons and microglia, particularly in AMPA receptor signaling pathways, were observed. These findings suggest that targeting AMPA signaling in PVT neurons could pave the way for innovative treatments and diagnostics aimed at improving outcomes for BD patients.\nThis paper has several strengths: Well written, addresses important questions, focuses on an understudied brain region, uses mostly state-of-the-art methods.\nThe biggest concern has to do with overinterpretation of results, with too little consideration for confounding and other limitations. There is also insufficient methodological detail in several areas.\nConcerns in order of importance:\n• Anatomical definition of PVT\n• CALB1 and CALB2 were used to define PVT, but they are not exclusive to PVT; CALB1 is broadly expressed in all ‘matrix’ nuclei; CALB2 is expressed in all midline and all intralaminar nuclei\n• Needs more detailed description (ideally with illustration) how the “medial thalamus” was dissected; e.g. more anterior part bordering anterior thalamic nuclei or more posterior part bordering mediodorsal thalamus\n• Gao et al. and Shima et al. have identified major differences within PVT, especially along anterior-posterior axis; the author only included a general scoring of cells relative to these mouse datasets, not specific to PVT subdivisions\n• Consideration of alternative explanations for differences in cell type composition\n• Dissection likely influences which excitatory populations were captured. How are cell populations represented across the 40 individuals?\n• Spatial transciptomics shows Exc_PVALB3, but no Exc_PVALB1/2. Why?\n• Is neuronal loss in BD limited to PVT (as defined by Exc_CALB cluster) or only detectable there due to more consistent dissection of the midline?\n• Overinterpretation of mechanistic link between excitatory neurons and microglia\n• SHISA9 and SYNDIG1 are both AMPAR auxiliary subunits, but an interaction between them is not known; SYNDIG1 encodes a transmembrane protein that is known to interact with AMPAR in neurons. Its function in microglia is not known. More evidence is needed of any claimed interactions in microglia.\n• Discussion would benefit from a more nuanced handling of the data, given its limitations, and more attention to areas where evidence is less clear and multiple conclusions are possible, e.g.,\n• Findings in PVT may not apply to the rest of the thalamus\n• Avoid anthropomorphic statements such as “Mice with recurrent depressive episodes”\n• Claims about high-frequency firing and role of KCNQ3 in potential response to high firing are too speculative; citations would be needed to support such a claim for this gene\n• Lack of relevance of thalamic organoid integration analysis\n• Authors did not sample the entire thalamus to claim that organoids recapitulate properties of PVT\n• This analysis is not really related to the rest of the paper\n• Insufficient methodological detail on Quality Control\n• Authors should describe QC and filtering process in more detail\n• Insufficient methodological detail on anchor-based integration\n• Was CCA or RPCA approach used, and why?\n• Insufficient exploration of confounding variables for DGE (and cellular composition)\n• Why were the covariates selected as described? The authors did not appear to have investigated other quality-related variables and explore variance they explained.\no Defining BD risk genes\no Reliance on MAGMA-based gene level scores from GWAS has significant limitations\no What about genes implicated by rare lof alleles in exome studies?\no Developmental origin of thalamic inhibitory neurons\no Only MGE, CGE, and LGE were considered, but it is not convincing that the observed thalamic inhibitory neurons were derived from LGE\no Existing literature points to midbrain-derived (SOX14+) inhibitory neurons\no Comparison with Siletti et al. is unclear, esp. regarding sample subject variation", "reviewer_4": "This manuscript investigates cellular and molecular alterations in bipolar disorder (BD) through single-nucleus RNA-seq (snRNA-seq) on postmortem brain tissue from BD patients and controls, focusing on the medial thalamus and frontal cortex. The authors first identify a subpopulation of thalamic excitatory neurons that represent neurons of the paraventricular thalamic nucleus (PVT), and provide a detailed molecular and spatial characterization of this cell population. Subsequent analyses reveal changes in cell type compositions and gene expression in BD. Intriguingly, these changes are more substantial in the thalamus compared to the cortex, highlighting the importance of this region and PVT neurons in particular. Overall, this study makes important contributions to our understanding of both the cytoarchitecture of the human PVT and alterations in BD at the cellular compositional and transcriptional levels. In addition, the analyses of the snRNA-seq data are thorough and I commend the authors for their careful assessments of QC metrics and the robustness of their results. There are several considerations that if addressed would strengthen the study.\nMajor comments:\n(1) Dysregulation of circadian rhythms is a feature of BD. Were any associations observed between the transcriptional changes in BD and patients' time of death?\n(2) The gene set enrichment analyses in this study use terms from the Gene Ontology (GO) database. Given that synaptic dysfunction is implicated in BD pathology, I recommend the authors also perform these analyses using the curated gene sets from the SynGO database (Koopmans et al. 2019), which were also used in O'Connell et al. (2025).\n(3) The cohort includes both BDI and BDII subtype cases. Is the dataset sufficiently powered to stratify analyses by subtype? I understand that the sample size of the BDII cases may preclude such analyses; if so, it would be worth mentioning this in the limitations section of the Discussion (together with sex differences in lines 449-451).\n(4) In Fig. 8a, downregulated genes in Tha_ExN_CALB are enriched in SNP-based genomic risk of all 3 phenotypes tested. Is this observed for other psychiatric disorders (such as post-traumatic stress disorder, attention deficit-hyperactivity disorder, borderline personality disorder, and autism spectrum disorder)?\n(5) In addition to the analyses of common variant-based BD genomic risk, it would also be helpful to comment on how these findings relate to the rare loss-of-function variants recently described in Thorgeirsson et al. (2025). For example, which cell types in thalamus and cortex express HECTD2 and AKAP11, and how does their expression in BD patients compare to controls?\n(6) While the data that support potential disruptions in PVT-microglia trans-cellular interactions in BD are compelling (MP7, Figs. 5b-c), it was unclear to me why the authors focused on SHISA9 in PVT neurons as a representative MP7 gene throughout the manuscript. Although SHISA9 is part of the core PPI network of BD-associated downregulated genes in PVT neurons (Fig. 5d), it appears to be expressed in only in a small fraction of Tha_ExN_CALB (PVT) neurons, while is expressed in most Tha_ExN_PVALB_1 neurons (Figs. 5d and 3a). I suggest that the authors either characterize this SHISA9-positive subpopulation of PVT neurons further (for example, are these detected in all thalamic samples?), or reduce/remove references to SHISA9 in PVT neurons throughout the text. These include parts of the section on \"Disturbed interaction between PVT neurons and microglia in BD\" (lines 251-284), the Discussion (lines 415-417), as well as the following sentence in the abstract: \"This interaction was further exemplified by the concurrent downregulation of SYNDIG1 in microglia and SHISA9 in PVT neurons, both of which are involved in AMPA receptor signaling pathways\" (top of pg. 3, line numbers not available).\nMinor comments:\n(1) Please provide details on the cohort in the Methods, including information on how individuals were selected for inclusion and the approval of the study protocol. Details are in the Reporting Summary, but must also be provided in the manuscript.\n(2) The final two paragraphs of the Discussion appear to be redundant."}, "author_response": "Point-by-point response\nWe would like to express our sincere gratitude to the reviewers for their encouraging and\nconstructive comments. Below is our point-by-point response to each reviewer’s specific\nfeedback. The reviewers' comments are presented in blue italics. The revised sections of our\nmanuscript are highlighted in red. The rectangles indicate the revisions made to the\nmanuscript. The references and tables have been renumbered to reflect the changes in the\nrevised submission.\nWe list the Supplementary materials as follows. We have deposited our count data to Gene\nExpression Omnibus (GEO) with accession number GSE306819. At this moment, only reviewers\nand editors can access these data with the following reviewer token and access URL.\nSupplementary Tables\nSupplementary Table 1. Sample information and sequencing metrics\nSupplementary Table 2. GO enrichment analysis by GSEA.\nSupplementary Table 3. Compositional changes in BD.\nSupplementary Table 4. Differentially expressed genes in BD.\nSupplementary Table 5. Robust DEG-GOs in BD (FDR < 0.05).\nSupplementary Table 6. Downregulated genes in Tha_ExN_CALB and Tha_Micro_1 with\nbiological interpretability\nSupplementary Figures\nSupplementary Figure 1. Cell-type annotation of thalamic and cortical cell clusters.\nSupplementary Figure 2. Characterization of thalamic excitatory neurons and PVT neurons.\nSupplementary Figure 3. Cell-type annotation in the thalamus.\nSupplementary Figure 4. Co-expression modules in PVT neurons.\nSupplementary Figure 5. Thalamic meta-gene programs estimated by geneNMF.\nSupplementary Figure 6. Cell-type annotation in the cortex.\nSupplementary Figure 7. Co-expression modules in cortical excitatory neurons.\nSupplementary Figure 8. HECTD2 and AKAP11 expression and thalamus-cortex meta-gene\nprograms estimated by geneNMF.\nSupplementary Figure 9. Thalamus-cortex communication estimated by LIANA.\nSupplementary Figure 10. Sample status and anatomical dissection.\nSupplementary Figure 11. Quality control (QC) for cell clustering.\nSupplementary Figure 12. CellBender-based ambient RNA decontamination.\nSupplementary Figure 13. NeuCyto clusters and additional spatial deconvolution.\nSupplementary Figure 14. Quality control (QC) for cell clustering and comparative analysis.\nSupplementary Figure 15. sccomp analyses incorporating PMI as an additional covariate.\nSupplementary Figure 16. DEG analyses incorporating PMI as an additional covariate.\nSupplementary Data\n01.DataS01.ExN_Tha_specific_genes.txt\n02.DataS02.InN_Tha_specific_genes.txt\n03.DataS03.PVT_Tha_specific_genes.txt\n1\n04.DataS04.scCODA.txt\n05.DataS05.IHC_count_data.txt\n06.DataS06.IHC_GLM.txt\n07.DataS07.DEG_GSEA.txt\n08.DataS08.PVT_module_GSEA_top10_GOs.txt\n09.DataS09.MAGMA_PsychiatricDisorders.txt\n10.DataS10.medication_signatures.txt\n11.DataS11.HBCA.profiles.txt\n12.DataS12.Celltype_markers.txt\n13.DataS13.thalamus.exp_seq_batches.cellcomposition.txt\n14.DataS14.cortex.exp_seq_batches.cellcomposition.txt\n15.DataS15.CompositionalAnalysis.GLM.AIC.txt\n16.DataS16.CompositionalAnalysis.GLM.Anova_ChiSq.txt\n17.DataS17.PutativeBDriskgenes.txt\n18.DataS18.Sccomp_PMI.txt\n19.DataS19.DEG_by_other_covariates.txt\n20.DataS20.DEG_Dx_withPMI.txt\n21.DataS21.DEG-GO_Dx_withPMI.txt\n22.DataS22.DEG_GSEA_withPMI.txt\n23.DataS23.lowPMI.DEG.txt\n24.DataS24.lowPMI.GSEA.txt\n25.DataS25.Mediation_signatures.references.txt\nExternal Data (figshare)\nReviewer access: https://figshare.com/s/dba4d20d4d7d02a67f9f\nBDsnRNAseq.thalamus.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nBDsnRNAseq.cortex.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nCellBender.output_pdf.zip\nCellBender.removedRNA.zip\nRepository data (GEO/SRA)\nGEO: GSE306819 (reviewer token: izenmomqlhihpqd)\nSRA: PRJNA1309835 (reviewer access:\nhttps://dataview.ncbi.nlm.nih.gov/object/PRJNA1309835?reviewer=t8fe31h5nf0gtpj229q7a6\n6d02)\nReview Materials\nThaFrC.CTRL.merge.harmony.QC.IDrenamed.rds (Seurat Object)\nthalamus.HBCAmerge.harmony.QC.IDrenamed.rds (Seurat Object)\nthalamus.merge.harmony.QC.IDrenamed.rds (Seurat Object)\ncortx.merge.harmony.QC.IDrenamed.rds (Seurat Object)\nDissection_PVT.wmv (video)\nNishioka et al. submit an interesting manuscript that investigates postmortem brain samples\nfrom bipolar disorder (BD) patients using single-nucleus RNA-seq. The authors analyzed a\ncohort of >80 samples of the prefrontal cortex and thalamus from ~40 subjects to identify new\ninsights into brain dysfunction in BD. The authors showed that gene expression and cell\ncomposition changes are predominant in the thalamus, and much lower in the cortex. This led\nthem to focus most of their analyses on the paraventricular thalamic nucleus. Subsequent\nanalyses showed that neurons of the paraventricular thalamic nucleus had large decrease in\ncell number in BD and exhibited the largest changes in gene expression. In particular, genes\ninvolved in synaptic processes were downregulated and few dysregulated hub genes were\nproposed. Using tools predicting cell-cell interaction changes, the authors proposed that\ndysfunction of microglia might mediate impairment of synaptic function in BD paraventricular\nthalamic nucleus neurons. Overall, the study is novel for the field and provides good insights\ninto specialization of thalamic excitatory neurons and their dysfunction in BD. Nevertheless,\nfew major issues, mainly concerning data analysis, should be addressed before accepting the\npublication.\nWe thank the reviewer for the thoughtful and encouraging assessment of our manuscript and\nfor recognizing the novelty and significance of our findings on thalamic excitatory neuron\ndysfunction in BD. We appreciate the constructive feedback on data analysis and have\naddressed all major points in the revised manuscript through methodological clarifications,\nadditional analyses, and new supplementary materials. We believe these revisions have\nstrengthened the rigor and clarity of our study.\nMajor issues:\n1. Covariate control needs significant improvement.\nWe appreciate the reviewer’s insightful comment and apologize for omitting the detailed\nrationale for covariate selection. In response, we have now added a comprehensive\ndescription of the covariate selection process and thoroughly revised the relevant section in\nthe Methods.\nCovariate Selection Strategy\nInitially, we considered a wide range of potential covariates: diagnosis (Dx), Sex, Age,\nPostmortem interval (PMI), brain pH (pH), race, cause of death, refrigerator delay, substances\nat death, experimental/sequencing batch. Given our modest sample size (N = 41), covariates\nshould be selected to minimize the risk of overfitting and multicollinearity.\n- Race: The cohort was heavily biased toward Caucasian individuals, with only one Asian\nsample. Thus, race was excluded from statistical modeling due to lack of variability.\n- Cause of death: A significant imbalance was observed, with suicide overrepresented\namong cases (Fisher’s exact test, P = 2.33 × 10-8). To reduce multicollinearity, this covariate\nwas excluded.\n- Refrigerator delay (RD): RD contained three missing values, which is substantial for a\n3\nsample of 41. Furthermore, RD was significantly correlated with PMI (Pearson’s correlation,\nP = 0.0181). We prioritized PMI as a relatively more relevant potential covariate.\n- Medication: Due to the wide heterogeneity of substances at death and missing information\nin several donors, subgrouping was not feasible with the available sample size. The details\nare described later.\n- Experimental and sequencing batches: We had 21 batches for nuclei isolation and 16\nsequencing batches in thalamic analysis. Modeling batch effects failed due to the large\nnumber of batches relative to sample size (N = 41). The details are described later.\nThe remaining candidate covariates were Dx, sex, age, PMI, and pH. Among them, pH was\nparticularly correlated with Dx (GLM P = 0.0074; see revised Supplementary Fig. 14b),\nsuggesting potential collinearity. Nonetheless, we tested models with all five covariates (Dx,\nSex, Age, PMI, and pH), as these are commonly used in postmortem brain studies.\nCompositional Analysis\nFor the compositional analysis, model fit was evaluated using Akaike Information Criterion\n(AIC) from generalized linear models (GLM) employing negative binomial distribution (glm.nb,\nMASS package), given that the sccomp package does not provide scores indicating model fit\nlike AIC. We compared multiple candidate models and ranked them by summed AIC rankings\n(rankSum) across both thalamic and cortical datasets. The thalamic microglia cluster was\nexcluded in this GLM-based calculation because the cell count distribution of this cluster did\nnot fit well with the negative binomial distribution.\nThe model including Dx, age, and sex (Dx + Age + Sex) exhibited the lowest AIC by rankSum,\nindicating the best overall fit (Supplementary Data 15). The contributions of PMI and pH\ndiffered between the thalamic and cortical analyses. We selected PMI as the secondary\ncovariate because pH showed a stronger correlation with diagnosis and represents a less direct\nmeasure of tissue integrity for both the thalamus and cortex (as described in the original\nmanuscript, Line 672). Based on these findings, we selected Dx + Age + Sex as the primary\ncovariate combination for compositional modeling, and Dx + Age + Sex + PMI as the secondary\ncovariate combination.\nIn line with the reviewer’s suggestion, we also assessed the relative contribution of each\ncovariate using Type III likelihood ratio chi-square tests (Anova(type = 3) from the car package).\nThe results are provided in Supplementary Data 16. Dx, Sex, and Age consistently\ndemonstrated significant contributions to compositional variation.\nDifferential Expression Analysis\nIn the differentially expressed gene (DEG) analysis, we adopted the same covariate\ncombination used in the compositional analysis. In response to the reviewer’s comments, we\ncalculated the AIC of GLMs for each gene’s log-transformed expression to specifically evaluate\nmodel fit when incorporating PMI and/or pH, as MAST does not provide scores indicating\nmodel fit like AIC. Gene expression was aggregated at the specimen level and limited to genes\nexpressed in at least 10% of cells within each cluster, consistent with the DEG criteria. This\nGLM-based modeling was performed solely to assess the impact of covariates; the DEG\nanalysis itself was conducted using zlm in the MAST package.\nfollowed by Dx + Age + Sex + PMI. This supports the use of Dx + Age + Sex as the primary\ncovariate combination and Dx + Age + Sex + PMI as a secondary model. In addition, we\nassessed the relative contribution of each covariate using Type III chi-square statistics\n(Anova(type = 3) from the car package) for both models. As these data are too huge to submit\nfor the journal web system, we deposited these data in figshare as External Data.\nWe have revised the Methods section to incorporate these details. The updated description\nenhances the clarity and transparency of our statistical approach.\nLine 897 in the revised manuscript\nCovariate selection in compositional analysis\nInitially, we considered a wide range of potential covariates: diagnosis (Dx), sex, age,\npostmortem interval (PMI), brain pH (pH), race, cause of death, refrigerator delay, substances\nat death, experimental batch, and sequencing batch. Given our modest sample size (N = 41),\ncovariates should be selected to minimize the risk of overfitting and multicollinearity. Race\nwas excluded from statistical modeling because the cohort was heavily biased toward\nCaucasian individuals, with only one Asian sample. Cause of death was excluded to reduce\nmulticollinearity, because it exhibited a significant imbalance with suicide overrepresented\namong BD cases (two-sided Fisher’s exact test, P = 2.33 × 10-8). Refrigerator delay contained\nthree missing values, and RD was correlated with PMI (Pearson’s correlation, P = 0.0181). We\nprioritized PMI as a relatively more relevant potential covariate. Due to the wide\nheterogeneity of substances at death and missing information in several donors, subgrouping\nbased on substances at death was not feasible in our dataset. Regarding experimental and\nsequencing batches, we had 21 batches for nuclei isolation and 16 sequencing batches in\nthalamic analysis. Modeling batch effects was unsuccessful due to the large number of\nbatches relative to the sample size (N = 41). Although experimental and sequencing batches\ncould not be incorporated into the subsequent statistical modeling due to their high\ndimensionality relative to sample size, we confirmed that the identified cell clusters were\nbroadly and consistently represented across batches (Supplementary Fig. 14a, Supplementary\nData 13 and 14).\nThe remaining candidate covariates were Dx, sex, age, PMI, and pH. Among them, pH\nvalues obtained from cerebellum samples, may not directly represent the specimen\nconditions of the thalamus and frontal cortex; given probable regional variability, we treated\npH as an indirect indicator for this study. Notably, pH and PMI were associated with Dx, with\npH showing a stronger signal (Supplementary Fig. 14b-c). These associations likely reflect\ndifferences in cause of death, as BD samples predominantly originated from individuals who\ndied by suicide, unlike control samples. Nonetheless, we tested statistical models with all five\ncovariates (Dx, Sex, Age, PMI, and pH), as commonly used covariates in postmortem brain\nstudies.\nFor the compositional analysis, model fit was evaluated using Akaike Information\nCriterion (AIC) from generalized linear models (GLM) employing negative binomial distribution\n(glm.nb in MASS [v7.3-60] package), because the sccomp and scCODA packages do not provide\nscores indicating model fit like AIC. In the GLM, total nuclei counts were normalized to 10,000,\n5\nwith each cluster’s count rounded to integers to reduce variability in total nuclei counts across\nindividuals. We compared multiple candidate models and ranked them by summed AIC\nrankings across both thalamic and cortical datasets (see rankSum in Supplementary Data 15).\nThe model including Dx, age, and sex (Dx + Age + Sex) showed the lowest overall AIC, indicating\nthe best model fit. While the contributions of PMI and pH differed between the thalamic and\ncortical analyses, we selected PMI as the secondary covariate because pH showed a stronger\ncorrelation with diagnosis and represents a less direct measure of thalamus and cortex tissues.\nBased on these findings, we selected “Dx + Age + Sex” as the primary covariate combination\nfor compositional modeling, and “Dx + Age + Sex + PMI” as the secondary covariate\ncombination. The relative contribution of each covariate was assessed by Type III likelihood\nratio chi-square tests (Anova(type = 3) from the car package, Supplementary Data 16).\nAs support for the above covariate selection, adding PMI did not increase the model\nfitness in the sccomp analysis, without changes of elpd_diff/se_diff > 5. More than five of\nelpd_diff/se_diff change is suggested for significant model improvement according to the\ndevelopers’ guideline42.\nLine 985 in the revised manuscript\nDifferentially expressed gene analysis\nDifferentially expressed genes (DEGs) were analyzed within each cell cluster in BD compared\nto controls using zlm function in MAST44 (v1.28.0), a statistical package for single-cell RNA-seq\nDEG analysis. The model included Dx, age, sex, and cngeneson (a covariate for cell quality) as\nfixed effects, and sample ID as a random effect. Sample-level covariates were selected based\non the prior compositional analysis and further validated by evaluating AIC values from GLMs\nfitted to log-transformed gene expression values. In this GLM-based model to select\ncovariates, gene expression was aggregated at the specimen level and restricted to genes\nexpressed in at least 10% of cells within each cluster, consistent with the downstream DEG\nfiltering criteria. These GLMs were used solely to assess covariate impact. Among the tested\ncovariate combinations, “Dx + Age + Sex” yielded the lowest overall AIC across gene-cluster\npairs, followed by “Dx + Age + Sex + PMI”, as determined by cumulative AIC rank sums\n(External Data, see Data Availability). These results support the use of “Dx + Age + Sex” as the\nprimary model and “Dx + Age + Sex + PMI” as a secondary model. The relative contribution of\neach covariate was further quantified using Type III chi-square statistics (Anova(type = 3) from\nthe car package), as reported in External Data. Following MAST guidelines, we incorporated\ncngeneson, a variable of scaled gene expression counts, to adjust for cell nuclei quality, where\nhigher gene detection rates suggest better nuclei quality within the same cluster.\nSupplementary Data\n15.DataS15.CompositionalAnalysis.GLM.AIC.txt\nAkaike information criterion (AIC) assessment for the GLM in compositional analysis.\n16.DataS16.CompositionalAnalysis.GLM.Anova_ChiSq.txt\nAnalysis of variance (ANOVA) Chi-square assessment for the GLM in compositional analysis.\nExternal Data (figshare)\nBDsnRNAseq.thalamus.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each thalamic cluster.\nBDsnRNAseq.cortex.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each cortical cluster.\n6\na. There have been multiple batches of sequencing, however, it is not described how many\nbatches, how many samples per batch, etc. Given that sequencing batch effect is one of the\nmajor confounders in snRNA-seq data, the authors should provide thorough analysis of\nsequencing batch impact on their clustering, cell type annotation, differential gene expression\nand cell composition. Supplementary Table 1 should also contain information regarding\nsequencing batch.\nWe appreciate the reviewer’s emphasis on addressing sequencing batch effects and apologize\nfor the omission of these details in the original manuscript. While we attempted to evaluate\nthe potential influence of experimental and sequencing batches for differential gene\nexpression and cell composition, the number of batches was too high to allow robust statistical\nmodeling. As shown in the revised Supplementary Table 1, the thalamic analysis involved 21\nbatches for nuclei isolation and 16 sequencing batches, including salvage sequencing. As a\nresult, GLM-based approaches failed to estimate the effects of these batch variables. Similarly,\nsccomp and scCODA did not perform when incorporating these batch covariates, and the zlm\nfunction in MAST either failed to converge or arbitrarily dropped variables, rendering\nestimation of batch effects infeasible.\nTo mitigate experimental variability, we limited nuclei isolation and 10x library preparation to\na single experimenter and consistently processed case-control pairs in parallel. We recognized\nthat nuclei isolation and library preparation were likely sources of bias and therefore\nprioritized strict procedural consistency at these stages. Due to the limited throughput of our\nsequencing platform (Illumina NextSeq 2000), which allowed only at most 3-4 libraries per run\nusing the P3 reagent (as the P4 reagent was not available at the time), the dataset comprises\nmany sequencing batches, including additional salvage runs to obtain sufficient read count per\nnuclei. Although pooling and sequencing all libraries together on a NovaSeq 6000 or others\nwould have been ideal for batch control, we opted against this approach due to the pace of\nour library preparation and concerns about potential degradation during prolonged storage.\nIn hindsight, a NovaSeq-based approach could have improved batch uniformity, although it\nwould not have addressed upstream variability introduced during nuclei isolation, a process\ninherently difficult to standardize due to its low-throughput, manual nature.\nGiven the limitations in incorporating batch variables into our statistical models, we instead\nevaluated potential biases in cell clustering and annotation across experimental and\nsequencing batches. As shown in Supplementary Fig. 14a and Supplementary Data 13/14, all\nidentified cell types were well represented across nearly all batches, indicating that batch\neffects did not drive cell-type clustering.\nIn addition to revising the covariate selection strategy (as described above), we have now\nincluded: Cell-type representation per experimental and sequencing batch, a detailed\nexplanation of our experimental design aimed at minimizing batch bias, and a discussion of\nbatch effects as a limitation of the study.\nIn this study, numerous experimental and sequencing batches were required due to limited\nresources. Future work would benefit from automated library preparation and concurrent\nsequencing of pooled libraries on high-throughput platforms to further control batch effects.\nLine 712 in the revised manuscript\nSingle-nucleus RNA-sequencing (snRNA-seq)\nNuclei isolated from the medial thalamus and frontal cortex underwent library preparation\nusing the 10X Chromium 3’ RNA v3.1 dual index protocol, following the manufacturer’s\ninstructions (10x Genomics, Pleasanton, USA). The target recovery rate was set at 5,000 nuclei,\nwith library amplification achieved through 14 PCR cycles. The snRNA-seq libraries were\nsequenced on the NextSeq2000 platform (Illumina, San Diego, USA) in accordance with the\nprotocols provided by Illumina and 10X Genomics. To minimize experimental bias, 10x\nGenomics library preparation was performed exclusively by the same experimenter (KI), with\nall case–control pairs processed in parallel. Due to the limitations of the NextSeq 2000\nplatform, which could accommodate only a few single-nucleus RNA-seq libraries per run,\nthe\ndataset spans multiple sequencing batches, including several salvage runs needed to achieve\nsufficient read count per nuclei.\nLine 908 in the revised manuscript\nRegarding experimental and sequencing batches, we had 21 batches for nuclei isolation and\n16 sequencing batches in thalamic analysis. Modeling batch effects was unsuccessful due to\nthe large number of batches relative to the sample size (N = 41). Although experimental and\nsequencing batches could not be incorporated into the subsequent statistical modeling due to\ntheir high dimensionality relative to sample size, we confirmed that the identified cell clusters\nwere broadly and consistently represented across batches (Supplementary Fig. 14a,\nSupplementary Data 13 and 14).\nSupplementary Figure 14. Quality control (QC) for cell clustering and comparative analysis.\na) Proportions of thalamic and cortical cell clusters across experimental batches.\nSupplementary Table 1 (abridged)\nSupplementary Data\n13.DataS13.thalamus.exp_seq_batches.cellcomposition.txt\nCell cluster composition by experimental and sequencing batches in the thalamus integration\nset.\n14.DataS14.cortex.exp_seq_batches.cellcomposition.txt\nCell cluster composition by experimental and sequencing batches in the cortex integration set.\nb. Medication is another large confounder in such analyses, since BD brains are usually exposed\nto a variety of drugs for a long period of life vs neurotypical brains. Potential impact of\nmedication should be addressed, such as tests whether different medications within their BD\ncohort have impact on gene expression, e.g. by comparing subgroups of patient samples,\ncomparing to other snRNA-seq studies, or to medication signatures identified in previous\nanimal/patient studies.\nWe agree with the reviewer that medication exposure is an important potential confounder in\npostmortem studies of psychiatric disorders. Before addressing the scientific concerns, we\nwould like to apologize for including individual-level medication data in the original submission.\nIn accordance with Canadian regulations on donor privacy, such individual data should not\nhave been included for publication. As Coroner data in Canada, where the brain tissue\nspecimens were collected are nominal and public, there is a theoretical risk of deanonymizing\nthe information by providing individual level data.\nIn the revised manuscript, we have instead included aggregated summary data, including a\nsummary of medication exposure (new Supplementary Table 1 and Supplementary Fig. 10a).\nGiven that the “Last 3 months medication” information reflects prescribed records rather than\nactual adherence, we have instead used “Substance at death” to represent medication status\nfor biological studies. In fact, the substances at death were not consistent with the record-\nbased last 3 months medication.\nacross donors substantially limited the feasibility of subgroup analyses in our sample (N = 41).\nAmong BD cases, no substance combination category included more than two donors,\nprecluding meaningful group-level comparisons. We also attempted compositional and\ndifferential expression analyses using only donors without any reported medication or\nsubstance exposure; however, only one BD donor met this criterion, rendering such\nstratification infeasible. As a result, detailed analysis of medication effects was not possible\nwithin this dataset.\nAs an alternative approach, we cross-referenced our DEGs with transcriptomic signatures from\npsychotropic medication studies in both human and animal models by examining the overlap\nbetween our DEGs and transcriptional changes induced by target psychotropic\nmedications/substances. The target medications/substances included those detected in more\nthan 20% of BD donors and at least twice as frequently as in controls (antipsychotics,\nantidepressants, and alcohol), as well as medications commonly prescribed for BD (lithium and\nvalproic acid). Because snRNA-seq data derived from the human thalamus were not available\nfor these medications/substances, we instead used bulk RNA-seq data from postmortem brain\nsamples, human induced neurons derived from induced pluripotent stem cells, or animal brain\nsamples, selecting the most appropriate report for each medication/substance. While this\ncomparison is not ideal, it represents the most feasible approach given current data\nconstraints.\nAs a result, we observed minimal overlap in most of the combinations, suggesting that the\nobserved transcriptional alterations are unlikely to be solely attributable to\nmedication/substance effects, except for the effect of olanzapine on the upregulated genes in\nthalamic excitatory neurons. At least, downregulated genes in PVT neurons, our focus in this\nstudy, are unlikely to be caused by psychotropic medication/substances in this analysis.\nWe have added a summary of this analysis in the revised Discussion and Methods sections,\na\nnew Fig.8g, and Supplementary Data 10\nLine 441 in the revised manuscript\nMedication signature\nMultiple medications and substances were detected in most postmortem brain samples at the\ntime of death (Supplementary Fig. 10a). As subgrouping based on substances at death or\nlifetime medication history was not feasible in our dataset, we assessed the potential overlap\nbetween medication signatures and the DEGs identified in our study. In the four neuronal\nclusters with the largest numbers of DEGs, there was generally no substantial overlap with\nmedication signatures for alcohol, antipsychotics, antidepressants, lithium, or valproic acid,\nexcept in a few instances (Fig. 8g, Supplementary Data 10). The upregulated genes in\nTha_ExN_CALB significantly overlapped with genes upregulated by olanzapine in medium\n10\nspiny neurons (MSNs)62 (Bonferroni-adjusted P = 0.00355, n of tests = 18 per cell type). A\nsimilar, though weaker, overlap was observed in Tha_ExN_PVALB_1 (adjusted P = 0.0772).\nAnother finding was the overlap between downregulated genes in FrC_ExN_L6ITCar3 and\ngenes downregulated by valproic acid in forebrain organoids63 (adjusted P = 0.0573). Although\nMSNs represent a distinct neuronal subtype from thalamic excitatory neurons, the consistent\noverlap between Tha_ExN_CALB/Tha_ExN_PVALB_1 and olanzapine-treated MSNs for\nupregulated genes suggests that these transcriptional changes could be influenced by lifetime\nmedication exposure. However, no notable overlaps were observed for other combinations.\nImportantly, the downregulated genes in Tha_ExN_CALB—which were a central focus of our\nresults—were unlikely to be attributable to medication or substance effects in this analysis.\nLine 646 in the revised manuscript\nDonors exhibited substantial heterogeneity in substances present at death, precluding\nstratification into subgroups of more than two individuals. Summary statistics of substances\npresent at death are presented in Supplementary Fig. 10a.\nLine 1155 in the original manuscript\nMedication signature\nWe evaluated the potential influence of medications and substances on DEGs by assessing\ntheir overlap with established medication signatures. A medication signature was defined as\nthe set of genes reported in previous studies to be differentially expressed in response to a\ngiven medication or substance. Target medications/substances included those detected in\nmore than 20% of BD donors and at least twice as frequently as in controls (alcohol98-100,\nantipsychotics62,101, and antidepressants102; Supplementary Fig. 10a), as well as medications\ncommonly prescribed for BD (lithium103 and valproic acid63). The reference data are\nsummarized in Supplementary Data 25. Because thalamic snRNA-seq data were unavailable\nfor these medications/substances, we instead used bulk RNA-seq data derived from\npostmortem brain tissue, human induced neurons generated from induced pluripotent stem\ncells, or animal brain samples, selecting the most appropriate publicly available dataset for\neach medication/substance. We used possible DEGs for this analysis to increase statistical\npower. This approach has inherent limitations due to brain region and cell-type differences,\nand the results should be interpreted as referential.\nFigure 8. Comparative analysis between the thalamus and cortex.\ng) Heatmap of P-values from hypergeometric tests for the overlap between DEGs in each cell\ntype and previously reported DEGs induced by specific medications/substances (medication\n11\nsignatures). Color indicates -log10(P-value) with Bonferroni-corrected significance threshold\nfor each panel (0.05/18). SSRI, selective serotonin reuptake inhibitor; MSN, medium spiny\nneuron; Bulk, bulk RNA-seq data; Cau, caudate; FrCExN, frontal cortex excitatory neuron;\nPFC,\nprefrontal cortex.\nSupplementary Table 1 (abridged)\nSupplementary Figure 10. Sample status and anatomical dissection.\na) Number of donors with detected substances at the time of death among 20 BD cases and\n19 controls, excluding those with missing data.\nSupplementary Table 1\nREMOVED: individual data including medication/substance information\nc. Some samples that were analyzed have high PMI, >24hours, and there are many studies\nshowing the impact of high PMI on the mRNA content, including benchmarking studies for\nRNA-seq experiments. The authors should perform a thorough analysis for potential\nconfounding impact of PMI and if necessary exclude samples that are too different from low\nPMI samples and interfere with data analysis.\nWe thank the reviewer for raising this important point. Prior to analysis, we also considered\nPMI to be a critical covariate, particularly as the majority of our samples had relatively high\nPMI values (>24 hours). Somewhat unexpectedly, however, incorporating PMI into the\n12\nstatistical models increased AIC values in both compositional and DEG analyses. Based on\nthese model fit metrics, we selected Dx + Age + Sex as our primary covariate combination,\nas previously described. Nonetheless, given the known importance of PMI in postmortem studies\nand the fact that Dx + Age + Sex + PMI consistently ranked as the second-best model, we\ndesignated it as our secondary covariate set. Although we cannot provide a definitive\nexplanation for this result, one possible interpretation is that the impact of PMI may be more\npronounced in samples with shorter postmortem intervals (e.g., PMI < 10 hours), where\nenzymatic degradation and other rapid postmortem processes are most active.\nIn response to the reviewer’s suggestion, we have added a more detailed assessment of the\nimpact of Dx, Age, Sex, and PMI in our statistical models (Supplementary Data 15/16, and\nExternal Data). The External Data is provided via figshare due to its large file sizes. The effects\nof Age, Sex, and PMI on cell composition are summarized in Supplementary Fig. 15b and\nSupplementary Data 18. The DEG results associated with each of these covariates are newly\nlisted in Supplementary Data 19 (in addition to the original Supplementary Fig. 16), enabling\nevaluation of covariate-specific effects.\nAs shown in these data, PMI had no significant effect on cell composition in the sccomp\nanalysis, while Dx, Age, and Sex were associated with changes in several cell types. The scCODA\nmodels incorporating PMI failed to converge reliably. In DEG analysis, Age and Sex exhibited\nconsiderable effects on gene expression, while PMI had less impact. Notably, PMI had no\nobservable effect on either cell composition or gene expression in the PVT neuron cluster\n(Tha_ExN_CALB), suggesting that our principal findings in PVT neurons are unlikely to be driven\nby PMI-related confounding.\nTo further confirm the robustness of our results, we performed an additional analysis\nrestricted to samples with relatively lower PMI values (<50 hours; n = 21). In the context of our\nsample set, a PMI of <50 hours represents a comparatively short interval, although we are\naware that this value is not considered as low in other studies. Nevertheless, both the\ncompositional shifts and DEG patterns remained consistent with our original results, further\nsupporting the validity of our conclusions.\nWhile the observed effects of Age and Sex on brain transcriptomes are intriguing, a more\ndetailed investigation is beyond the scope of the current study. We are currently pursuing this\ndirection in a follow-up project.\nLine 1088 in the revised manuscript\nExtended analysis incorporating PMI as an additional covariate for robustness check\nAs described above, we extendedly analyzed compositional and transcriptional changes by\nincorporating a covariate of PMI into the main covariates (Dx, age, and sex) to check the\nrobustness of our results. This extended analysis yielded results consistent with our initial\nfindings (Supplementary Fig. 15 and 16). Specifically, we observed a compositional decrease\nof thalamic excitatory neurons in BD, including PVT neurons (FDR = 0.0115, FC = 0.582),\nand\nFrC_ExN_L6ITCar3 (FDR = 0.0257, FC = 0.578) (Supplementary Fig. 15a, Supplementary Data\n18). Dx, sex, and age affected the cellular compositions, while PMI did not (Supplementary\nFig. 15b). The scCODA model did not perform reliably when PMI was included alongside Dx,\n13\nage, and sex, yielding inclusion probabilities of 1.0 for all variables across all clusters, even\nif using total cell count normalization. This likely reflects overparameterization due to the\nnumber of covariates relative to our sample size.\nIn this extended analysis, Dx, sex, and age prominently affected cell composition and\ngene expression, while PMI less affected them (Supplementary Fig. 15b and Supplementary\nData 19). The DEGs by age or PMI were defined by those with absolute value of FC > 1.1 for\nthe period of two standard deviations in our dataset (i.e. 26.7 years for age, and 38.7 h for\nPMI). Notably, PMI showed no detectable impact on either cell composition or gene\nexpression in the PVT neurons (Tha_ExN_CALB), suggesting our key findings in PVT neurons\nare unlikely to be confounded by PMI.\nWhen incorporating PMI into the statistical models, the Dx-based DEG-GOs were again most\nobserved in PVT neurons (Supplementary Fig. 16a). The DEGs in PVT neurons (Supplementary\nFig. 16b, Supplementary Data 20) continued to exhibit enrichment in synaptic and neuronal\nGO terms with the same top GOs—transporter complex (GO:1990351) and monoatomic ion\nchannel complex (GO:0034702)—including SHISA9, CACNA1C, and KCNQ3 (Supplementary\nFig. 16c, Supplementary Data 21). Again, the downregulated DEGs were linked to the co-\nexpression module Excit-M4, while upregulated DEGs were linked to Excit-M11 and Excit-M9\n(Supplementary Fig. 16d). Among the DEGs identified, 43 of 64 downregulated DEGs (67.2%),\nand 100 of 157 upregulated DEGs (63.7%) from the primary analysis were also classified as\nDEGs in this extended analysis (Supplementary Data 20). Notably, downregulated BD risk\ngenes—SHISA9, CACNA1C, and CTNN5—remained detectable as downregulated DEGs and\nconstituted core PPI networks in this extended analysis (Supplementary Fig. 16e). SHISA9 was\nalso detected as a downregulated DEGs in Tha_ExN_PVALB_1 in the extended analysis.\nThe\nSNP-based BD genomic risk was again most enriched in the downregulated genes in PVT\nneurons (P = 0.00419), whereas upregulated genes included no BD risk genes despite a greater\ntotal number of DEGs (69 downregulated vs. 110 upregulated). The DEG-GOs in the frontal\ncortex cell clusters were virtually not observed, consistent with the primary analysis\n(Supplementary Fig. 16f). Finally, GSEA performed on DEGs from this extended analysis\nyielded results similar to those of the primary analysis (Supplementary Data 22), with\ndownregulated genes in PVT neurons again enriched for synaptic and neuronal GO terms.\nLine 1125 in the revised manuscript\nExtended analysis to the samples with low PMI\nWhile PMI had no apparent effect on our findings of PVT neurons, we analyzed the samples\nwith relatively low PMI (< 50 h), consisting of 12 BD cases and 9 controls, to check the\nrobustness of our findings for PVT neurons and thalamic microglia. The covariates were set as\nthe same to the primary analysis: Dx, age, and sex. PVT neurons and thalamic excitatory\nneurons robustly exhibited compositional decreases (FC = 0.563, FDR = 0.0223 for PVT\nneurons; FC = 0.541, FDR = 0.0193 for thalamic excitatory neurons) by sccomp analysis. DEGs\nin PVT neurons were also consistent with the main analysis, showing comparable FCs\n(Supplementary Data 23). For instance, KCNQ3, SHISA9, and CACNA1C showed FC values of\n0.755, 0.840, and 0.843, respectively. GSEA also yielded similar results, with the top enriched\nterm being potassium ion transmembrane transporter activity (GO:0015079), including KCNQ3\n(Supplementary Data 24). In this analysis, we employed GSEA instead of DEG-based GO\nenrichment, as the smaller sample size (12 BD cases vs. 9 controls) reduced the statistical\npower and FDR control.\na) Compositional changes of cell clusters in the thalamus (left) and the cortex (right), with\nlog2(Fold-change) on the x-axis and -log10(FDR) on the y-axis. Circle size corresponds to the\nnuclei count within each cluster, and color represents major cell classes. b) Effect (c_effect) of\neach covariate on cell clusters, as estimated by sccomp. The endpoints of the lines represent\nthe lower and upper bounds of the c_effect estimate. Clusters with a false discovery rate (FDR)\n< 0.05 are indicated by red points and lines.\na) DEG and DEG-GO count across thalamic and cortical clusters. b) Volcano plot of the genes\nexpressed in Tha_ExN_CALB, displaying log2(Fold-change) on the x-axis and DEG robustness\nscore on the y-axis. c) Downregulated (left) and upregulated (right) robust DEG-GOs in the\nthalamus, plotted with fold-change on the x-axis and robustness score on the y-axis. Circle\nsize reflects the DEG count within each GO, with colors indicating respective cell clusters.\nAmong robust DEGs, only enriched GO terms among possible DEGs are displayed for fold-change\ncalculation. d) Proportions of co-expression modules in downregulated and upregulated DEGs\nin Tha_ExN_CALB. e) PPI network of downregulated possible DEGs in Tha_ExN_CALB. Colored\nnodes indicate overlap with DEGs from the primary analysis. f) Upregulated robust DEG-GOs\nin the cortex, plotted with fold-change on the x-axis and robustness score on the y-axis. Circle\nsize reflects the DEG count within each GO, with colors indicating respective cell clusters. No\ndownregulated robust DEG-GO was observed in this analysis.\nSupplementary Data\n16\n15.DataS15.CompositionalAnalysis.GLM.AIC.txt\nAkaike information criterion (AIC) assessment for the GLM in compositional analysis.\n16.DataS16.CompositionalAnalysis.GLM.Anova_ChiSq.txt\nAnalysis of variance (ANOVA) Chi-square assessment for the GLM in compositional analysis.\n18.DataS18.Sccomp_PMI.txt\nsccomp analysis including postmortem interval (PMI) as an additional covariate.\n19.DataS19.DEG_by_other_covariates.txt\nDEGs by sex, age, and PMI, using PMI as an additional covariate.\n20.DataS20.DEG_Dx_withPMI.txt\nDEGs by disease status, using PMI as an additional covariate.\n21.DataS21.DEG-GO_Dx_withPMI.txt\nGene ontology (GO) enrichment of DEGs by disease status, using PMI as an additional\ncovariate.\n22.DataS22.DEG_GSEA_withPMI.txt\nGSEA results for thalamic and cortical DEGs incorporating PMI as an additional covariate.\n23.DataS23.lowPMI.DEG.txt\nDEGs identified using only low-PMI samples.\n24.DataS24.lowPMI.GSEA.txt\nGSEA results using only low-PMI samples.\nExternal Data (figshare).\nBDsnRNAseq.thalamus.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each thalamic cluster.\nBDsnRNAseq.cortex.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each cortical cluster.\nd. Additionally, the authors should statistically analyze potential impact of all metadata\nvariables, including age, sex, drug exposure etc.\nWe thank the reviewer for raising this important issue. As described earlier, we have added\nanalyses to evaluate the potential impact of age, sex, and PMI. The effects of these covariates\non cell composition are summarized in Supplementary Table 3 and Supplementary Data 16,\nand their effects on gene expression are also detailed in Supplementary Data 19. As the files\nfor the relative contribution of each covariate for gene expression are too large, we have\nprepared them in figshare as External Data.\nvariety of substance use across samples and incomplete information made it infeasible to\nstatistically model medication effects within our dataset. Instead, as noted above, we\ninvestigated the potential influence of pharmacological exposure by comparing our\ntranscriptional findings to known expression signatures of major psychiatric medications.\nPlease, refer to the above responses.\nLine 935 in the revised manuscript\nThe relative contribution of each covariate was assessed by Type III likelihood ratio chi-square\ntests (Anova(type = 3) from the car package, Supplementary Data 16).\nLine 997 in the revised manuscript\nThe relative contribution of each covariate was further quantified using Type III chi-square\nstatistics (Anova(type = 3) from the car package), as reported in External Data (figshare).\nLine 1099 in the revised manuscript\nAs described above, we extendedly analyzed compositional and transcriptional changes by\nincorporating a covariate of PMI into the main covariates (Dx, age, and sex) to check the\nrobustness of our results. This extended analysis yielded results consistent with our initial\nfindings (Supplementary Fig. 15 and 16). Specifically, we observed a compositional decrease\nof thalamic excitatory neurons in BD, including PVT neurons (FDR = 0.0115, FC = 0.582),\nand\nFrC_ExN_L6ITCar3 (FDR = 0.0257, FC = 0.578) (Supplementary Fig. 15a, Supplementary Data\n17). Dx, sex, and age affected the cellular compositions, while PMI did not (Supplementary\nFig. 15b).\nLine 1101 in the revised manuscript\nIn this extended analysis, Dx, sex, and age prominently affected cell composition and\ngene expression, while PMI less affected them (Supplementary Fig. 15b and Supplementary\nData 19). The DEGs by age or PMI were defined by those with absolute value of FC > 1.1 for\nthe period of two standard deviations in our dataset (i.e. 26.7 years for age, and 38.7 h for\nPMI). Notably, PMI showed no detectable impact on either cell composition or gene\nexpression in the PVT neurons (Tha_ExN_CALB), suggesting our key findings in PVT neurons\nare unlikely to be confounded by PMI.\nSupplementary Table 3 (abridged)\n16.DataS16.CompositionalAnalysis.GLM.Anova_ChiSq.txt\nAnalysis of variance (ANOVA) Chi-square assessment for the GLM in compositional analysis.\n19.DataS19.DEG_by_other_covariates.txt\nDEGs by sex, age, and PMI, using PMI as an additional covariate.\nExternal Data (figshare).\nBDsnRNAseq.thalamus.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each thalamic cluster.\nBDsnRNAseq.cortex.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each cortical cluster.\n2. Although authors did include filter of ambient RNA, it should be improved.\na. There is some ambiguity with the analysis of detection and removal of ambient RNA. In\nMethods Line 565:\nDecontX was chosen over SoupX and no ambient RNA decontamination\nDo the authors refer to “no ambient RNA decontamination” as a tool or package? The authors\nshould be more specific.\nWe apologize for the ambiguity. In this context, “no ambient RNA decontamination” was\nintended to mean that no ambient RNA decontamination tool or package was applied. In the\nrevised Methods section, we have clarified this sentence.\nWe selected DecontX for ambient RNA removal in preference to SoupX85, CellBender86, or\nanalysis performed without any ambient RNA decontamination package (no_decont_package)\nfrom the following reasons. SoupX and no_decont_package could not effectively distinguish\nNeuCyto clusters (details are described later) from oligodendrocytes in the UMAP projection,\nprobably reflecting the effect of ambient RNA to clusters expressing a small number of genes\n(Supplementary Fig. 11d).\nSupplementary Figure 11. Quality control (QC) for cell clustering.\nd) UMAP projection displaying thalamic nuclei clusters processed with SoupX decontamination\nand without any ambient RNA decontamination package (no_decont_package).\nb. Overall, given potential contamination with ambient RNA, more robust analysis is necessary.\nWe propose to use CellBender and provide the following data in Supplementary figures:\n- Sample-wise expression of ambient RNAs\n- CellBender ELBO (test and train) plots\n- CellBender cell calling curves displaying cell probabilities of droplets\nWe appreciate this comment. CellBender has been reported as one of the most effective tools\nfor ambient RNA removal, as demonstrated in Janssen et al. (Genome Biology 2023) and\nCaglayan et al. (Neuron 2022).\nmasked cell types in brain single-nuclei datasets”, Neuron 110:4043-4056, 2022) indicated that\nneuronal cell clusters are difficult to decontaminate even with CellBender, as follows.\n“Strikingly, CellBender did not reduce ambient RNA contamination from the neurons (Figures\nS6C and S6E). As expected, non-nuclear ambient RNA markers NRGN and CHN1 levels were\nhigher in more contaminated neuronal nuclei, whereas nuclear-retained MALAT1 levels were\nlower (Figures S6D and S6F). Contamination patterns were similar among the previously\nannotated neuronal subtypes of NSD1, indicating that ambient RNA contamination in neurons\nis a cell-type agnostic problem and, unlike glia, is not accounted for by CellBender (Figures S7A\nand S7B). Other ambient RNA contamination removal tools were also more effective in glia\nthan neurons, indicating a general deficiency in the current methods to remove ambient RNA\ncontamination from the dominant cell type in the tissue (Figures S7C and S7D).“ (p.4049)\nThus, there is currently no definitive method for ambient RNA decontamination in snRNA-seq\ndatasets focusing on neurons. In our analysis, we heuristically selected the most appropriate\ndecontamination approach based on clustering results, acknowledging that there is no\nabsolute ground truth for ambient RNA removal in our study.\nWe initially applied CellBender to our dataset, recognizing its potential advantages for ambient\nRNA removal. We have deposited “Sample-wise expression of ambient RNAs”, “CellBender\nELBO (test and train) plots”, and “CellBender cell calling curves displaying cell probabilities of\ndroplets” in External Data. Unexpectedly, it did not perform well with our thalamic samples.\nIn several cases, CellBender appeared to fail in clearly distinguishing nuclei-containing droplets\nfrom empty droplets in several samples, as illustrated in Supplementary Fig. 12a and files in\nExternal Data.\nThe downstream integration of CellBender-processed data resulted in atypical clusters\n(Tha_ExN_PVALB_X and Tha_InN_X), which were annotated as putative excitatory or\ninhibitory neurons based on marker gene expression but exhibited abnormally low total RNA\ncounts (Supplementary Fig. 12b-g). These atypical clusters are not interpreted as cytoplasmic\ncomponents. In contrast, no such atypical clusters were observed in other pipelines using\nDecontX, SoupX, or no RNA decontamination packages, even if increasing cluster resolution.\nWe interpreted these atypical clusters in the CellBender-based pipeline as likely artifacts\nintroduced by miscalculation during ambient RNA removal.\nThis suboptimal performance of CellBender might stem from the presence of cytoplasmic\nfragment-containing droplets in our dataset. Based on the presence of NeuCyto clusters and\nprobable cytoplasmic components in microscopic observation, we classify droplets into three\ncategories:\n(A) droplets containing intact nuclei,\n(B) droplets containing cytoplasmic components without nuclei\n(C) droplets containing only liquid ambient RNA without nuclei or cytoplasmic components.\nAlthough type B droplets are not ideal for nuclear transcriptomics, they seemed to exist in our\nsamples and HBCA dataset, likely due to the fragility of postmortem human brain tissue—\n21\nparticularly in thalamic samples. In fact, TapeStation QC during library preparation also\nindicated greater RNA degradation in thalamic than cortical tissue. To preserve nuclear\nintegrity, we employed a gentler homogenization protocol using 0.05% Tween-20 instead of\nthe more commonly used 0.1% Triton X-100. This adjustment improved RNA integrity, yielding\nhigher RIN scores during pilot testing. Additionally, to avoid laser-induced RNA damage, we\nrelied on centrifugation rather than fluorescence-activated cell sorting (FACS) for nuclei\npurification. While these protocols helped preserve RNA quality, they may also contribute to\nthe retention of cytoplasmic fragments (type B droplets), which correspond to NeuCyto\nclusters in our analysis. While microscopic observations were limited to our samples, we\nnoticed NeuCyto clusters in HBCA samples. This implies that cytoplasmic components may also\nexist in other snRNA-seq datasets derived from the human brain.\nCellBender’s underlying model assumes a binary classification between nuclei-containing (A)\nand empty (C) droplets. Consequently, type B droplets—particularly those enriched in\nneuronal cytoplasmic RNA—may be difficult to resolved to A or C. This difficulty could bias\nambient RNA estimation and result in the removal of biologically meaningful transcripts.\nAlthough this provides a plausible explanation for CellBender’s suboptimal performance, we\nacknowledge that the exact cause of the aberrant clusters remains uncertain.\nGiven this uncertainty at the time of initial submission, we did not include CellBender results\nin the original manuscript. Additionally, we had limited opportunities to optimize CellBender\nparameters due to the availability of minimal GPU resources. Processing a single snRNA-seq\nlibrary sometimes required many hours on conventional CPUs, making parameter tuning\nimpractical without adequate GPU resources.\nIn response to the reviewer’s comments, we have now included these results in the revised\nversion. While cytoplasmic fragments likely contributed to CellBender’s limitations in our\ndataset, we acknowledge that this interpretation remains speculative and have therefore only\nminimally mentioned it in the main text.\nLine 764 in the revised manuscript\nCellBender was applied using the expected cell count estimated by Cell Ranger along with\nother default parameters. However, in several cases, CellBender appeared to fail in clearly\ndistinguishing nuclei-containing droplets from empty droplets (Supplementary Fig. 12a,\nExternal Data [figshare]). CellBender-based clustering resulted in atypical clusters that\nexpressed neuronal markers but exhibited extremely low total RNA counts (Supplementary\nFigs. 12b-g). These clusters, labeled as Tha_ExN_PVALB_X and Tha_InN_X did not correspond\nto any known thalamic cell types and were not consistent with cytoplasmic droplets, as\ndetermined by nuclear fraction and mitochondrial RNA content. Although a definitive ground\ntruth for cell-type identity was not available, we interpreted these clusters as likely artifacts\nresulting from inappropriate ambient RNA removal, leading to biologically implausible\ntranscriptomic profiles. In contrast, cell-type annotations were consistent across DecontX,\nSoupX, and the no_decontamination_package, and no unexplained anomalous clusters\nemerged under these approaches, even if increasing clustering resolution. While the exact\ncause of this phenomenon remains uncertain, we suspect that the presence of cytoplasmic\ncomponents—discussed later and supported by NeuCyto clustering—may interfere with\ndroplets.\nSupplementary Figure 12. CellBender-based ambient RNA decontamination.\na) Representative plots illustrating unsuccessful cell probability inference by CellBender, with\nbarcodes ranked by total UMI counts. b) UMAP projection of thalamic cell nuclei clustered\nfollowing CellBender-based ambient RNA removal. c) Zoomed-in view of excitatory neuron\nclusters from panel b, highlighting Tha_ExN_PVALB_X using distinct coloring. d) Zoomed-in\nview of inhibitory neuron clusters from panel b, highlighting Tha_InN_X. e) Dot plot showing\ncell-type marker expression across thalamic clusters identified using CellBender. Dot size\nindicates the proportion of expressing nuclei, and color intensity represents average\nexpression levels. f) Nuclear fractions and mitochondrial RNA fractions across CellBender-\ndefined clusters. g) Number of detected genes per cluster across four RNA decontamination\n23\nmethods: CellBender, DecontX, SoupX, and no decontamination. Red and blue triangles\nindicate Tha_ExN_PVALB_X and Tha_InN_X, respectively.\nExternal Data\nCellBender.output_pdf.zip\nCellBender output PDF files for all samples.\nCellBender.removedRNA.zip\nCellBender removed RNA data for all samples.\nIn addition to the arguments for CellBender, we have revised the supplementary figures\nillustrating the possible presence of cytoplasmic fragments (presumably neuronal in origin).\nLine 789 in the revised manuscript\nClusters designated as NeuCyto were considered neuronal cytoplasmic components from the\nfollowing properties: low nuclear fraction, high content of mitochondrial RNA (Supplementary\nFig. 13a), high expression of neuronal markers, virtually no expression of other glial markers\n(Supplementary Fig. 1a, 3a, and 6a), and distinctive clustering from probable empty droplets\n(Supplementary Fig. 13b-c). The probable empty droplets, designated as “Empty”, were\ndefined as droplets filtered out by Cell Ranger with detected genes of 100 – 300. This range\nwas set by QC procedure and barcode rank plot for UMI counts by Cell Ranger. NeuCyto\ndroplets were robustly detected from 40/41 of our thalamic samples (occupying 1.46% of total\ncell nuclei) and 40/41 of our cortex samples (occupying 1.46% of total cell nuclei) as well as\n20/20 of HBCA thalamic samples (occupying 1.66% of total cell nuclei). The clear separation of\nVGLUT+ and GAD+ NeuCyto clusters in the thalamus indicated that they were distinct\ncytoplasmic components (Fig. 3a). We removed NeuCyto clusters for the subsequent analysis\nas non-nuclei droplets but retained for UMAP visualization and dot plots as a technical\nreference. Indeed, in the final nuclei suspension, we observed objects stained with Nonyl\nAcridine Orange (NAO) but negative for Propidium Iodide (PI) and Hoechst 33342\n(Supplementary Fig. 13d), likely representing mitochondria-rich, nucleus-free cytoplasmic\nfragments corresponding to the NeuCyto clusters.\nLine 861 in the revised manuscript\nFor reference, we also mapped the location of thalamic NeuCyto clusters. The NeuCyto\nclusters within the thalamus could be further categorized into three sub-clusters based on\nmarker expression (Supplementary Fig. 13g). Spatial deconvolution revealed that\nNeuCyto_VGLUT2_CALB was relatively localized in the PVT region, while\nNeuCyto_VGLUT2_PVALB was sparsely distributed across other areas (Supplementary Fig.\n13h). NeuCyto_GAD clusters seemed relatively located in proximity to inhibitory neurons.\nThe\nthree NeuCyto sub-clusters likely represented neuronal cytoplasm from excitatory and\ninhibitory neurons.\na) Nuclear fractions and mitochondrial RNA fractions across samples in the thalamic\nintegration dataset (see Fig. 3a). b) UMAP projection illustrating “Empty” clusters (ambient\nRNA, 1,000 droplets per sample) and canonical cell nuclei clusters shown in Fig. 3a. c) UMAP\nprojection of the same clusters as in panel b, with distinct coloring to highlight\nTha_NeuCyto_VGLUT, Tha_NeuCyto_GAD, and Empty clusters. d) Microscopic images of\ncytoplasmic fragments (white arrowheads) in two representative spots. Cytoplasmic\nfragments are indicated by NAO staining (green) in the absence of PI or Hoechst 33342 signals.\nCell nuclei appear magenta due to combined PI and Hoechst 33342 staining. e) CytoSPACE\ndeconvolution based solely on our snRNA-seq dataset (excluding the HBCA dataset). f)\nCytoSPACE deconvolution using only the CT06 snRNA-seq dataset. g) UMAP projection\n25\ndisplaying three NeuCyto sub-clusters within the thalamus. h) Predicted spatial distribution of\nthe three NeuCyto sub-clusters. Color intensity in the side bars reflects the estimated number\nof cells per spatial spot. For panels g and h, only control samples were used to prevent\npotential bias from BD-associated alterations in spatial deconvolution.\n3. Data integration\nThe authors used iterative integration with Seurat and Harmony. There is no justification\nprovided for why Seurat-only or Harmony-only integration are not enough. Simple comparison\nplots will help the reader to understand the rationale for choosing iterative integration.\nWe apologize for the lack of clarity in our original description. Our approach, referred to as\n“iterative integration” in the original manuscript, consists of the following two-step process:\n1. We first merged the target data and performed Seurat normalization, scaling, and PCA\nusing all nuclei, including doublets and low-represented clusters (i.e., clusters detected in\n<75% of samples). Clustering and doublet identification are then conducted on this dataset\nafter Harmony batch-correction. At this stage, the scaled expression values still reflect\ncontributions from doublets and potentially artifactual clusters.\n2. After removing doublets and low-represented clusters, we re-normalized, re-scaled, and\nre-clustered the dataset using Seurat and performed Harmony batch-correction again. This\nsecond round ensures that scaling and downstream analyses are no longer biased by\ndoublet-contaminated or artifactual clusters.\nWe adopted this two-step procedure because omitting the second round would allow doublets\nand low-represented clusters to influence scaling and clustering, potentially distorting batch\ncorrection and downstream analyses such as differential expression.\nIn retrospect, we acknowledge that the terms “iterative” and “integration” were ambiguous.\nOur intention was not to imply simultaneous use of both Seurat and Harmony batch correction\nmethods. Rather, we used only Harmony for batch correction in both rounds, and the word\n“integration” referred to the consecutive processing of merging, normalization, scaling, PCA,\nand batch correction. To improve clarity, we have revised the terminology in the manuscript:\n“iterative integration” has been removed. A word “batch-correction” is now used to describe\nmore specifically to describe how Harmony was used.\nLine 741 in the revised manuscript\nIntegration of snRNA-seq data\nCount data from Cell Ranger were decontaminated using the DecontX tool in the celda\npackage80 (v1.16.1) and subjected to multiplet identification via scDblFinder81 (v1.14.0). These\nsteps were performed for each count data in our and the HBCA6 dataset prior to integration.\nAfter decontamination and multiplet labeling, only nuclei with > 300 and < 15,000 detected\ngenes and < 5% reads that mapped to mitochondrial genome were retained. The resulting\ncount data were processed by Seurat (v5.0/v5.1)82 with Harmony83 (v1.2.0) batch-correction.\nMultiplets were defined as cells labeled by scDblFinder, as well as those within clusters\nexhibiting a high fraction of scDblFinder-identified multiplets and expressing markers of\n26\nmultiple major cell types (e.g., oligodendrocytes and astrocytes) (Supplementary Fig. 11a-b).\nAfter removing multiplets and low-represented clusters, each dataset was re-processed using\nSeurat with Harmony batch-correction. This two-step procedure was consistently applied\nacross the study to minimize the influences of doublets and low-represented clusters. PCA\ndimensions and clustering resolution parameters were tuned to obtain a reasonable number\nof clusters for each dataset. Clusters from < 75% of samples were removed as low-represented\nclusters after both the first and second rounds of clustering.\nThe details of the scripts are available at Github repository.\nProvided that the word “integration” in the reviewer’s comment means “batch-correction”,\nthe integration (batch-correction) method in our main pipeline is Harmony-only. Seurat-only\nintegration (batch-correction) would correspond to RPCA-based batch-correction using Seurat.\nThe results of this procedure were described in the original manuscript. RPCA (reciprocal\nprincipal component analysis) batch-correction is a native algorithm in Seurat package without\nother external package. A simple comparison was described in the main text and\nSupplementary Fig. 11.\nLine 756 in the original manuscript\nFor batch-correction, Harmony was preferred over Seurat RPCA, because the latter failed to\nclearly separate thalamic and cortical OPCs/astrocytes (Supplementary Fig. 11c), a result\nmisaligned with findings by Siletti et al.6\nSupplementary Figure (not changed from the original manuscript set)\nSupplementary Figure 11. Quality control (QC) for cell clustering.\nc) UMAP projection of major cell types (left) and their respective brain region origins (right)\nfor the thalamus-cortex integration dataset, using Seurat RPCA integration.\nThis is contrasted to Fig. 1b (below)\n4. Differential gene expression\na. We could not find information whether p values for differentially expressed genes were\nadjusted. The authors should outline how they did the adjustment.\nWe apologize for the omission and confusion caused by our original description. In the initial\nmanuscript, we reported only false discovery rate (FDR) values for differentially expressed\ngenes (DEGs) in the main text. All p-values from the differential expression analysis using MAST\nwere adjusted for multiple comparisons using the Benjamini-Hochberg method to control the\nFDR. These FDR-adjusted values and the original non-adjusted P-values are provided in the\noriginal version of Supplementary Table 4 (column M and J) . To improve clarity, we have now\nadded a description of the FDR calculation method in the revised manuscript.\nLine 1002 in the revised manuscript\nDEG analysis was limited to protein-coding genes with Entrez IDs to streamline downstream\nanalyses. Clusters with fewer than 1,000 nuclei were excluded to ensure sufficient counts per\nindividual. All the P-values from the MAST analysis were adjusted for multiple comparisons\nusing the Benjamini-Hochberg method to control the false discovery rate (FDR). The raw P-\nvalues and FDRs are provided in Supplementary Table 4.\nSupplementary Table 4 (abridged, not changed from the original manuscript set)\n28\nb. While GO analysis of DE genes is often used in snRNA-seq studies, it works better for larger\ngene lists than the authors obtained in this study. GSEA is better applicable here, since GSEA\nimplements all genes, not only short list of significantly changed DE genes.\nWe agree that GSEA is generally a powerful and widely used method for DEGs. We initially\nconsidered using GSEA, but opted for a GO enrichment approach based on the following\nrationales.\n1. Primary objective and DEG quantification\nOne of the primary aims of this study was to identify cell types exhibiting substantial\ntranscriptional changes, as represented by the number of DEGs per cluster. This required us to\ndefine and count DEGs explicitly. As DEG calling inherently depends on thresholding, we sought\nto minimize arbitrariness by adopting the approach used in our original manuscript: applying\na range of thresholds and evaluating the robustness of DEG classification based on the\nfrequency with which genes are identified as DEGs across thresholds. Given that DEGs were\ndefined as discrete entities through the above procedure, a conventional over-representation\nanalysis using GO terms was a natural and direct downstream analysis. This allowed us to\ninterrogate the biological functions of explicitly defined gene sets.\n2. Limitations of single-metric ranking in GSEA:\nGene Set Enrichment Analysis (GSEA) requires ranking genes based on a single metric, typically\nfold change (FC). However, in our dataset, FC alone may not reliably capture gene relevance\ndue to substantial inter-individual variability and our moderate sample size (N = 41),\nwhich\nmakes the analysis susceptible to outlier effects. In addition, genes with low expression rates\n(e.g., ER = 0.1) can exhibit unstable or spurious FC values, further complicating interpretation.\nWhile we considered composite ranking approaches that integrate FC, FDR, and ER (e.g., via\nmultiplication), we ultimately opted against them due to the inherent subjectivity and\narbitrariness involved in selecting weighting schemes.\nDespite the above considerations, we acknowledge that GSEA remains a robust and\ninformative approach and can serve as a valuable complement to traditional GO enrichment\nanalysis. In response to the reviewer’s suggestion, we have incorporated GSEA based on FC as\na complementary analysis. Notably, the GSEA results for the downregulated genes in PVT\nneurons are relatively consistent with the DEG-GO enrichment findings, further supporting the\nrobustness of our results. In contrast, the GSEA results for the upregulated genes in PVT\nneurons are less consistent with the DEG-GO enrichment findings, indicating a difference in\nbiological interpretability. We again thank the reviewer for this insightful suggestion.\nThese GSEA results are provided in Supplementary Data 7 and summarized in the revised\nResults section. To illustrate the notable genes identified by DEG-GO and GSEA analyses,\nwe\nhave revised Supplementary Table 5 and created a new Supplementary Table 6. In response\nto the Reviewer 4’s suggestion, we have also incorporated SynGO enrichment analysis. We\nhave also performed GSEA for secondary analysis.\nLine 193 in the revised manuscript\nThis analysis revealed 56 downregulated and 69 upregulated probable DEGs in the medial\nthalamus (Fig. 3e, Supplementary Table 4). Tha_ExN_CALB (PVT neurons) displayed the\nhighest DEG count, with 22 downregulated and 43 upregulated probable DEGs (64\ndownregulated and 157 upregulated possible DEGs). Tha_ExN_PVALB_1 exhibited the second-\nlargest number of upregulated DEGs, with 81.4% of these genes overlapping with those\nupregulated in Tha_ExN_CALB. There were 29 downregulated and 67 upregulated robust DEG-\nGOs in the thalamus, detected only in Tha_ExN_CALB, Tha_ExN_PVALB_1, and Tha_Micro_1\n(Fig. 3f, Supplementary Table 5). Tha_ExN_CALB also demonstrated the highest number of\nDEG-GOs, 24 downregulated and 58 upregulated. The robust DEG-GOs were predominantly\nenriched in neuronal and synaptic processes, and several DEG-GOs demonstrated overlap\nbetween downregulated and upregulated categories in Tha_ExN_CALB, such as the regulation\nof trans-synaptic signaling (GO:0099177) and modulation of chemical synaptic transmission\n(GO:0050804). Given the enrichment of synaptic terms in the DEG-GOs, we additionally\nperformed SynGO enrichment analysis45. The downregulated DEGs in the Tha_ExN_CALB\ncluster again showed the strongest SynGO enrichment, including three terms related to\nsynaptic structure: integral component of postsynaptic membrane (GO:0099055), integral\ncomponent of presynaptic membrane (GO:0099056), and integral component of postsynaptic\ndensity membrane (GO:0099061) (Supplementary Table 5).\nComplementary GSEA based on fold changes supported the enrichment of neuronal\nand synaptic processes among downregulated genes in Tha_ExN_CALB (Supplementary Data\n7). GSEA identified potassium channel activity, including KCNQ3, KCND2, and KCNC2 as the top\nGO term for the downregulated genes in Tha_ExN_CALB. The overlap between DEG-GO and\nGSEA results (GSEA-GO) was more pronounced for downregulated than upregulated genes\n(9/24 vs. 6/58; two-sided Fisher’s exact test, P = 0.00931, odds ratio = 5.07) (Supplementary\nTable 5). Notably, the top three DEG-GO terms among downregulated genes—transporter\ncomplex (GO:1990351), monatomic ion channel complex (GO:0034702), and main axon\n(GO:0044304)—were also identified by GSEA. The downregulated genes in Tha_ExN_CALB\nexhibit greater biological coherence and interpretability across enrichment methods.\nSupplementary Table 6 summarizes the downregulated probable DEGs indicated by both DEG-\nGO and GSEA as biologically interpretable genes for BD. Among them, KCNQ3, CACNA1C,\nSHISA9, KCNC2, LRRC7, ROBO2, and CNTNAP4 were consistently enriched across all three\nanalytical approaches (DEG-GO, GSEA, and SynGO).\nLine 269 in the revised manuscript\nIn PVT neurons, the downregulated genes showed consistent GO term enrichment in both GO\nanalysis and GSEA, and the MAGMA-based BD genomic risk was relatively prominent for\ndownregulated genes. Thus, we focused on the downregulated genes. When extending the\nPPI analysis to possible DEGs, the downregulated possible DEGs in PVT neurons also\ndemonstrated robust PPI networks (enrichment P < 1.0 x 10⁻¹⁶), encompassing genes within\nthe top three DEG-GOs as biological processes (BP) or cellular components (CC) (Fig. 4d). These\nkey genes, including CACNA1C and SHISA9 as putative BD risk genes, were identified as core\nDEGs based on their PPI and GO enrichment. Particularly, KCNQ3, GRIN2A, and GRM5\nintersected with the top three BP GOs—modulation of chemical synaptic transmission\n(GO:0050804), regulation of trans-synaptic signaling (GO:0099177), and regulation of\nmonoatomic ion transmembrane transport (GO:0034765). KCNQ3 and KCNC2 intersected\nwith the top three CC GOs—transporter complexes (GO:1990351), monoatomic ion channel\ncomplexes (GO:0034702), and main axon (GO:0044304). These synaptic and ion-channel-\nrelated genes as PPI network hubs aligned with the biological pathways indicated by BD\ngenomic studies1,2,4 and our DEG-GO/GSEA/SynGO results. Particularly, the overlap of KCNQ3\nin both BP and CC top GOs underscores its significant relevance to BD, besides putative BD risk\ngenes.\nLine 311 in the revised manuscript\nSupplementary Table 6 summarizes the downregulated genes in Tha_ExN_CALB and\nTha_Micro_1 that were enriched in one or more GO terms in both DEG-GO and GSEA analyses.\nNotably, these genes belonged to the same co-expression modules in each cluster.\nLine 1030 in the revised manuscript\nWe also performed FC-based GSEA separately for the downregulated and upregulated gene\nsets to further evaluate the robustness of the identified biological signals, using gseGO\nfunction in clusterProfiler with the enrichmentScore threshold of 0.7.\nLine 1121 in the revised manuscript\nFinally, GSEA performed on DEGs from this extended analysis yielded results similar to those\nof the primary analysis (Supplementary Data 22), with downregulated genes in PVT neurons\nagain enriched for synaptic and neuronal GO terms.\nLine 1133 in the revised manuscript\nGSEA also yielded similar results, with the top enriched term being potassium ion\ntransmembrane transporter activity (GO:0015079), including KCNQ3 (Supplementary Data\n24). In this analysis, we employed GSEA instead of DEG-based GO enrichment, as the smaller\nsample size (12 BD cases vs. 9 controls) reduced the statistical power and FDR control.\nSupplementary Table 5 (abridged)\nSupplementary Table 5. Robust DEG-GOs in BD (FDR < 0.05). Add GSEA overlap information\n(column G).\nSupplementary Table 6. Downregulated genes in Tha_ExN_CALB and Tha_Micro_1 with\nbiological interpretations\nSupplementary Data\n07.DataS7.DEG_GSEA.txt\nGene set enrichment analysis (GSEA) results for thalamic and cortical differentially expressed\n32\ngenes (DEGs).\n22.DataS22.DEG_GSEA_withPMI.txt\nGSEA results for thalamic and cortical DEGs incorporating PMI as an additional covariate.\n24.DataS24.lowPMI.GSEA.txt\nGSEA results using only low-PMI samples.\nIn contrast, cortical clusters showed less enrichment by GSEA.\nLine 349 in the revised manuscript\nOur DEG analysis revealed 55 downregulated and 73 upregulated probable DEGs in the frontal\ncortex (Fig. 6c, Supplementary Table 4). FrC_ExN_L5IT_4 and FrC_ExN_L6Car3 exhibited the\nhighest numbers of probable DEGs (Supplementary Table 4) and 11 robust DEG-GOs\n(Supplementary Table 5). While the DEG-GOs primarily pertained to synaptic and neuronal\nfunctions, their number in the frontal cortex was lower than that observed in the thalamus\n(Fig. 6d vs. Fig. 3d). Among cortical DEG-GO, only GO:0045211 (postsynaptic membrane), an\nupregulated GO term in FrC_ExN_L6ITCar3, was identified by FC-based GSEA, and no SynGO\nenrichment was observed. These findings suggest that cortical DEGs exhibit lower biological\ninterpretability compared to thalamic DEGs.\n5. Integration with genetics is rather preliminary and more sophisticated tools can be\nimplemented for integration with the latest GWAS BD study. This will help validating the\nfindings in the study for their relevance to BD. The thalamus and cortex should be analyzed\nseparately by integration with GWAS data to show potential regional specialization with\nrespect to BD-linked brain dysfunction.\nFinally, at the moment, it seems that most association of DE genes in this study with BD risk\nfactors is based only on 1 gene from 36 credible genes from the last GWAS study - SHISA9.\nWe would caution the authors for not to overinterpret such association.\nWe appreciate the reviewer’s suggestion to strengthen the genetic integration. As correctly\nnoted, interpreting GWAS results based solely on MAGMA is preliminary. The recent PGC4\nstudy (O’Connell Nature 2025) identified 36 “credible genes” through integrative approaches,\nof which only SHISA9 overlapped with the genes identified in our analysis. We acknowledge\nthat only one gene among the 36 credible genes was associated with our results. We should\nbe cautious in interpreting this result.\nHowever, given the polygenic nature of BD, it is unlikely that only these 36 genes are involved\nin its pathophysiology. Although no definitive method currently exists for exhaustive gene\nprioritization, it remains important to consider a broader set of BD-associated gene candidates\nand assess their overlap with DEGs. Our original intent in using MAGMA was to provide a\npreliminary estimate of SNP-based genomic risk, and we now clarify this motivation in the\nrevised manuscript.\nIn line with the reviewer’s comment that MAGMA alone is insufficient for robust gene\nprioritization, we have now incorporated additional integrative approaches to complement the\nMAGMA findings:\n- Fine-mapping using PolyFun+SuSiE/FINEMAP\n- Transcriptome-wide association study (TWAS) using FUSION\nWe used fine-mapping to narrow down the putative BD risk genes from MAGMA analysis.\nTWAS data was not used to narrow down the putative BD risk genes, but only to support the\nfindings of CACNA1C, because TWAS data is not comprehensive across all brain regions,\nlacking\nthalamic data. These methods were employed to refine candidate BD-associated genes.\nNevertheless, we emphasize that the results from these approaches remain putative, and the\nidentified genes should be interpreted with appropriate caution, given the inherent challenges\nof mapping disease-associated genes from GWAS.\nAccordingly, we have removed the previously described list of “164 BD risk genes” derived\nfrom MAGMA. While we acknowledge the preliminary nature of MAGMA, we have retained\nits results to provide a rough estimate of SNP-based genomic risk and to address Reviewer #4’s\nsuggestion regarding the extension of MAGMA analyses to other psychiatric disorder GWAS\ndatasets. Finally, as further recommended, we have conducted and visualized separate\nanalyses for the thalamus and cortex to improve anatomical specificity.\nWe have revised the following points. Throughout the manuscript, we replaced “BD risk genes”\nwith “putative BD risk genes.”\nLine 223 in the revised manuscript\nUpon closer inspection of probable DEGs, MAGMA-based BD genomic risk1,46 was most\nenriched in the downregulated DEGs in Tha_ExN_CALB (P = 0.00123). Additional fine-mapping\nby PolyFun+SuSiE/FINEMAP47-49 further supported CACNA1C and SHISA9 as putative BD risk\ngenes within the downregulated set. Especially, SHISA9 is one of the 36 credible genes for BD\nin the largest GWAS to date1. While CACNA1C is not listed in these 36 credible genes, CACNA1C\nhas consistently ranked among the top genes by p-value in prior PGC1-4 studies1,2,50,51, and\nwas also implicated by our FUSION-based TWAS integration. In contrast, the upregulated\nprobable DEGs in Tha_ExN_CALB exhibited weaker enrichment for MAGMA-based BD\ngenomic risk (P = 0.0916), including no putative BD risk genes despite the greater number of\nupregulated genes (n = 22 for downregulated vs. 43 for upregulated, Supplementary Table 4).\nLine 387 in the revised manuscript\nWe examined which cell types most expressed putative BD risk genes. Thalamic excitatory\nneurons, including PVT neurons, showed the highest expression of the putative BD risk genes\nconcerning the number of expressed genes and their average expression rate (Fig. 8b). For\nexample, CACNA1C and SHISA9 exhibited relatively higher expression in thalamic excitatory\nneurons than in cortical excitatory neurons (Fig. 8c). However, the overall differences were\nmodest, and both genes showed the highest expression in Tha_ExN_PVALB_1, suggesting that\nadditional factors likely contribute to the transcriptional alterations observed in PVT neurons,\nas discussed later.\nLine 1034 in the revised manuscript\nAssessment of BD risk genes\nTo evaluate the enrichment of SNP-based genomic risk for BD, we first applied MAGMA46 to\nroughly assess the overlap between BD GWAS signals and probable DEGs, using multi-ancestry\nPGC4 BD summary statistics1. Enrichment analyses were conducted for cell clusters containing\nfive or more probable DEGs (n = 14 clusters). To assess specificity, we extended the analysis to\nGWAS summary statistics from other recent Psychiatric Genomics Consortium (PGC) datasets,\nincluding schizophrenia (SCZ)33, major depressive disorder (MDD)59, attention-\ndeficit/hyperactivity disorder (ADHD)88, post-traumatic stress disorder (PTSD)89, obsessive-\ncompulsive disorder (OCD)90, substance use disorder (SUD)91, eating disorder (ED)92, panic\ndisorder (PD)93, and autism spectrum disorder (ASD)94 (https://pgc.unc.edu/).\nBeyond assessing SNP-based genomic enrichment, we further prioritized specific BD\nrisk genes. As MAGMA alone provides only preliminary gene-level association signals, we\nadopted a multi-step approach to define “putative BD risk genes”. First, we designated the “36\ncredible genes” reported in the PGC4 BD study1 as the most reliable BD-associated genes\ncurrently available. Given the highly polygenic architecture of BD, it is unlikely that only 36\ngenes account for all the BD genetic basis. Additional BD risk gene candidates were defined as\nthose overlapping between MAGMA-significant genes and those harboring fine-mapped SNPs\nidentified by PolyFun+SuSiE/FINEMAP47-49, using the multi-ancestry PGC4 BD summary\nstatistics1. MAGMA identified 167 genes associated with BD at Bonferroni-corrected\nsignificance (adjusted P < 0.05, based on 17,487 genes with Entrez IDs). Fine-mapping was\nconducted following the developer’s pipeline, using UKBB_LD as the reference panel. Fine-\nmapped SNPs were annotated using promoter regions (PsychENCODE TSS±2kb; file:\ntss.sites.codingOnly.gencode.v19.annotation.bed), enhancer regions (PsychENCODE GRN2;\nfile: INT-14_ElasticNet_Filtered_Cutoff_0.1_cross_assembly_GRN_2.csv), and SnpEff-4.3 for\ncoding consequences. We then defined “putative BD risk genes” as the union of the 36 credible\ngenes and 75 MAGMA-based genes containing fine-mapped SNPs in promoters, enhancers, or\nas missense variants, resulting in a total of 107 genes (four overlapped, Supplementary Data\n17). To explore regional and cell-type specificity, we evaluated the expression of these 107\ngenes across thalamic and cortical cell clusters using two metrics: (1) the number of putative\nBD risk genes expressed per cluster (expression rate > 0.1), and (2) the mean expression level\nof those genes.\nTo further support gene prioritization, we conducted transcriptome-wide association\nstudies (TWAS) using FUSION95 with precomputed weights from the GTEx v8 brain datasets\n(http://gusevlab.org/projects/fusion/). Due to the absence of thalamic samples in GTEx, the\nTWAS results were used only for complementary reference and were not included in the\ndefinition of putative BD risk genes. TWAS-significant genes were identified using Bonferroni\ncorrection for each precomputed weight set (adjusted P < 0.05, across 2,557-7,732 tested\ngenes per dataset). Given that GTEx represents bulk RNA-seq data from multiple brain regions\n(not including thalamus), we interpreted these results broadly, without strict anatomical\nspecificity.\nFigure 4. Transcriptional disturbances in PVT neurons.\nc) PPI networks highlighting downregulated and upregulated probable DEGs in PVT neurons.\nNode color indicates the direction of gene regulation (Red, down; blue, up). Rectangles\nindicate BD risk genes. PPI enrichment P-value < 1.0 × 10-16. d) PPI networks highlighting\ndownregulated possible DEGs in PVT neurons. Rectangles represent BD risk genes, and colors\ndenote the genes in top DEG-GOs. Font size indicates DEG-GO enrichment. Note that the two\nPPI networks are identical but have different focuses on the genes. PPI enrichment P-value <\n1.0 × 10-16.\nFigure 8. Comparative analysis between the thalamus and cortex.\nb) Expression of putative BD risk genes across cell nuclei clusters in the thalamus and cortex\nset (see Fig. 3a and 6a). The x-axis represents the count of expressed putative BD risk genes,\nand the y-axis denotes their average expression rates (percent expressed) among the\nrespective clusters. The genes are limited to protein-coding genes with expression rates > 0.1\nwithin each cluster.\nSupplementary Data\n17.DataS17.PutativeBDriskgenes.txt\nPutative bipolar disorder risk genes identified in this study and their sources.\nWe have also revised the relevant phrasing (from “and” to “and/or”) to reflect the putative\nnature of BD risk genes, acknowledging that KCNQ3 is not included among the putative BD risk\ngenes.\nAbstract in the revised manuscript\nNotably, PVT neurons showed significant downregulation of genes associated with synaptic\nfunction and ion channel regulation, including SHISA9, CACNA1C, and KCNQ3—genes linked\nto putative BD risk and/or implicated as critical hubs in the downregulated molecular\nnetworks.\nLine 40 in the revised manuscript\nIn PVT neurons, SHISA9, CACNA1C, and KCNQ3 exhibited notable changes as putative BD risk\ngenes and/or critical hubs in the downregulated synaptic/neuronal gene networks.\n6. Compositional analysis\nWhile it is justified to use sccomp package for compositional analysis, the package does not\ncontrol for inter-dependency of proportional changes, and the authors admit it when\ncomparing to immunohistochemical assessment, Lines 169-172:\nNote that immunohistochemistry directly assessed cell density, while compositional analysis\nusing snRNA-seq data assessed relative proportional changes. The immunohistochemical\nassessment indicated that the compositional increase in oligodendrocytes and astrocytes in BD\nshould be mainly derived from prominent reductions of excitatory neurons, not reflecting\nabsolute cell density increases.\nTo avoid such issue, there are packages that allow for independent compositional analysis of\ncell proportions, such as scCODA. The authors should use one of such packages to cross-\ncompare snRNA-seq with immunohistochemistry.\nWe agree with the reviewer’s point and have now conducted a parallel compositional analysis\nusing scCODA. As correctly noted, compositional methods rely on specific assumptions, and\napplying multiple complementary methods is important to ensure robustness. We have now\nincluded the scCODA results in the revised manuscript. To define an appropriate reference,\nwe\nfirst selected oligodendrocyte precursor cells (OPCs), as they exhibited no significant case-\ncontrol differences in both the thalamus and cortex and appeared relatively stable across\nsamples. We also tested oligodendrocytes as the reference to compare snRNAseq to\nimmunohistochemistry (IHC) data.\nNotably, the scCODA results were concordant with those obtained using sccomp in the\nthalamus, while not concordant in the cortex. In cortex analysis, total nuclei counts were\nnormalized to 10,000 per individual for scCODA, because scCODA did not work well with raw\ncounts. This supports our argument that the thalamus exhibited more notable BD-associated\nchanges than the cortex. We include these scCODA results in a new Supplementary Data file\nand have updated the Results and Discussion sections accordingly. Unfortunately, scCODA did\nnot perform reliably when PMI was included alongside diagnosis, age, and sex, yielding\ninclusion probabilities of 1.0 for all variables across all clusters, even if using normalized counts.\nThis unlikely reflects the true phenomena, which might be due to too many covariates for\nscCODA analysis with our sample size.\nLine 160 in the revised manuscript\nComplementary analysis using scCODA43 also supported a consistent compositional reduction\nin thalamic excitatory neurons, including PVT neurons (Supplementary Data 4). In the scCODA\nanalysis, Tha_ExN_CALB and Tha_ExN_PVALB_1 exhibited comparable compositional\ndecreases (FC = 0.686 and 0.663).\nLine 182 in the revised manuscript\nIf we set Tha_Oligo_1 as a reference cluster in scCODA analysis, assuming that the absolute\ncounts of oligodendrocytes were comparable between BD and control from IHC results, we\nconsistently observed compositional decreases in thalamic excitatory neurons\n(Supplementary Data 4).\nLine 339 in the revised manuscript\nThe scCODA analysis was not consistent to the sccomp analysis, exhibiting the compositional\ndecreases of FrC_ExN_L2/3IT_1 (FC = 0.797), FrC_ExN_L5IT_4 (FC = 0.866), and\nFrC_InN_CGE_LAMP5 (FC = 0.848). In contrast to findings in the thalamic excitatory neurons,\nthe cortical neurons showed inconsistent results and smaller reductions (Fig. 6b vs. Fig. 3b).\nLine 885 in the revised manuscript\nAs compositional analysis methods rely on specific assumptions, we applied two methods of\na compositional analysis for single-cell data, sccomp42 (v2.1.6) and scCODA43 (v0.1.9), to\ncompare nuclei counts across cell clusters between BD and control groups. The sccomp\nemploys a sum-constrained beta-binomial model, while the scCODA employs a Bayesian\nmodel using a reference cluster. The threshold of inclusion probability in scCODA was set as\n0.95. We considered the analysis to have failed when the observed inclusion probability was\n1.0 for all clusters within a given covariate, indicating a lack of model convergence or\ninformativeness. In cases where raw counts could not be used successfully (cortex analysis\nand PMI analysis), total nuclei counts were normalized to 10,000 per individual, and cluster-\nspecific counts were rounded to integers to reduce inter-individual variability in total nuclei\ncounts. The compositional analysis was limited to the clusters in six major classes:\nexcitatory/inhibitory neurons, oligodendrocytes, OPCs, astrocytes, and microglia.\nLine 1095 in the revised manuscript\nThe scCODA model did not perform reliably when PMI was included alongside Dx, age, and\nsex, yielding inclusion probabilities of 1.0 for all variables across all clusters, even if using total\ncell count normalization. This likely reflects overparameterization due to the number of\ncovariates relative to our sample size.\nSupplementary Data\n04.DataS4.scCODA.txt\nscCODA results for thalamic and cortical clusters, including estimated fold-changes and\nposterior probabilities for each covariate.\n7. Cell-cell interaction\n38\na. Based on the text lines 270-283, it sounds as the authors propose that SHISA9 protein on\nneurons interacts with SYNDIG1 on microglia? If yes, then some literature search should help\nconfirming this potential interaction. Otherwise, we suggest concentrating interpretation on\ngeneral discussion on how microglia modulate synaptic processes in excitatory neurons via\nAMPA receptor complex. Additionally, LIANA analysis should help with such discussion,\nhowever the data is just shown in Figure 5e, without any interpretation.\nAs correctly pointed out, there is currently no direct evidence in the literature that SHISA9 and\nSYNDIG1 proteins physically interact, and the original manuscript contained an\noverinterpretation in retrospect. Studies such as Kalashnikova et al. (Neuron 2010) and\nKhodosevich et al. (Neruon 2014) describe SHISA9 and SYNDIG1 proteins as auxiliary\nmembrane proteins associated with AMPA receptor (AMPAR) complexes. While both proteins\nare involved in regulating synaptic function, there is no published evidence indicating a direct\nmolecular binding between them. We have now revised this section to accurately reflect\ncurrent knowledge. The revised statement emphasizes that microglia might modulate synaptic\nprocesses in excitatory neurons through pathways involving the AMPA receptor complex,\nwithout suggesting a direct interaction between SHISA9 and SYNDIG1.\nAbstract in the original manuscript\nREMOVED: This interaction was further exemplified by the concurrent downregulation of\nSYNDIG1 in microglia and SHISA9 in PVT neurons, both of which are involved in AMPA receptor\nsignaling pathways.\nLine 279 in the original manuscript\nREMOVED: The coordinated downregulation of SYNDIG1 in microglia and SHISA9 in excitatory\nneurons exemplifies a probable disturbance in microglia-neuron interactions for synapses in\nthe BD thalamus, indicating a plausible pathological mechanism of SHISA9, one of the 36\nputative BD-causal genes1.\nLine 415 in the original manuscript\nREMOVED: The downregulation of SYNDIG1 in microglia and SHISA9 in PVT neurons might be\nan example of dysregulated molecular interactions within the Excit-Micro meta-gene programs\nassociated with BD.\nLine 314 in the revised manuscript\nAt the gene level, SYNDIG1 and SHISA9 encode auxiliary subunits for AMPA receptor (AMPAR)\ncomplex56,57, with both genes’ robustness scores exceeding 0.7 (Supplementary Table 4).\nNotably, SYNDIG1 was prominently expressed in microglia, while SHISA9 was primarily\nexpressed in excitatory neurons (Fig. 5d). SHISA9 consisted of core PPI networks of BD-\nassociated downregulated genes in PVT neurons (Fig. 4d) and overlapped in the\ndownregulated genes in Tha_ExN_PVALB_1, as one of the 36 putative BD-causal genes1.\nBeyond microglia and excitatory neurons, thalamic T cells were included in this meta-gene\nprogram. Fig. 5e illustrates potential ligand-receptor interactions between neurons and\nmicroglia/T cells, as analyzed through LIANA58. Although we did not detect concurrently\ndownregulated DEGs as ligand–receptor pairs, KCND2 in Tha_ExN_CALB, LAMB1 in\nTha_ExN_PVALB_2, and KCNQ1 in Tha_Micro_1 were identified as downregulated DEGs\namong the LIANA-based ligand-receptor pairs. These genes encode ion channels or\nextracellular matrix proteins that may potentially contribute to Excit–Micro interactions.\nHowever, as this analysis is an inference, the precise ligand–receptor mechanisms remain to\nbe elucidated in future studies.\nLine 556 in the revised manuscript\nBeyond cell-type-based analysis, we identified trans-cellular meta-gene programs between\nexcitatory neurons and microglia (Excit-Micro meta-gene programs), likely reflecting\nphysiological synaptic interactions between these cell types. Microglia are known to modulate\nsynaptic plasticity and neuronal circuitry through dynamic interactions with neurons76.\nNotably, such meta-gene programs in the thalamus are disrupted in BD, further characterized\nby a significant decrease in synaptic gene expression in both PVT neurons and microglia. Ling\net al. reported downregulation of the synaptic neuron and astrocyte program in\nschizophrenia10; our findings highlight downregulation of the synaptic genes in the thalamic\nneuron and microglia in BD. The concurrent downregulation of SYNDIG1 in microglia and\nSHISA9 in excitatory neurons suggests disruptions of AMPAR signaling, as both genes encode\nauxiliary subunits of the AMPAR complex56,57. TIAM1, another downregulated gene in\nmicroglia, has also been implicated in regulating synaptic AMPAR function77. The preferential\nlocalization of SYNDIG1-positive microglia in cortical gray matter is consistent with this\ninterpretation. Nonetheless, these meta-gene programs remain data-driven constructs, and\nthe roles of SYNDIG1 in microglia—as well as their mechanistic links to SHISA9 in PVT\nneurons—are currently unknown. Future studies employing animal and in vitro cellular\nmodels will be crucial for elucidating these mechanisms.\nAdditionally, LIANA does not support a SHISA9-SYNDIG1 interaction. In the manuscript, LIANA\nwas employed to investigate potential microglia-excitatory neuron interactions suggested by\nNMF-based module analysis. However, no ligand-receptor pair was identified in which both\nthe ligand and receptor genes were downregulated.\nNevertheless, several DEGs were detected that may participate in intercellular signaling as\neither the ligand or receptor: KCND2 in Tha_ExN_CALB, LAMB1 in Tha_ExN_PVALB_2, and\nKCNQ1 in Tha_Micro_1. These genes encode ion channels or extracellular matrix proteins that\nmay indirectly contribute to cell-cell interactions. However, the precise ligand-receptor\nmechanisms remain unresolved, and we now explicitly state that cell-cell communication\nanalysis remains a subject for future investigation.\nLine 320 in the revised manuscript\nFig. 5e illustrates potential ligand-receptor interactions between neurons and microglia/T\ncells, as analyzed through LIANA58. Although we did not detect concurrently downregulated\nDEGs as ligand–receptor pairs, KCND2 in Tha_ExN_CALB, LAMB1 in Tha_ExN_PVALB_2, and\nKCNQ1 in Tha_Micro_1 were identified as downregulated DEGs among the LIANA-based\nligand-receptor pairs. These genes encode ion channels or extracellular matrix proteins that\nmay potentially contribute to Excit–Micro interactions. However, as this analysis is an\ninference, the precise ligand–receptor mechanisms remain to be elucidated in future studies.\ne) Predicted ligand-receptor interactions between excitatory neurons and microglia/T cells.\nDot plots represent expression magnitudes and interaction specificities, with interactions\nlabeled as [ligand gene] -> [receptor gene]. Blue or red stars indicate ligand–receptor pairs in\nwhich the ligand or receptor is a DEG. Left: ligands expressed in microglia or T cells targeting\nreceptors in thalamic excitatory neurons. Right: ligands expressed in thalamic excitatory\nneurons targeting receptors in microglia or T cells.\nb. One missing opportunity in the study is analysis of changes in cortico-thalamic connectivity\nin BD. This can be done using similar cell-cell interaction tools.\nThank you for your excellent suggestion. In the original manuscript, we refrained from\npresenting cell-cell communication between cortical and thalamic clusters because we did not\nidentify meta-gene programs supporting such cortico-thalamic interactions. Nevertheless,\ninteractions between cortical and thalamic cells could be of interest to readers.\nAlthough our snRNA-seq approach does not permit direct inference of long-range connectivity,\nLIANA can provide indirect evidence by evaluating potential ligand-receptor interactions\nbetween cortical and thalamic cells. Given the large number of possible cortical-thalamic cell-\ntype combinations, we restricted our analysis to deep-layer cortical excitatory neurons and\nthalamic excitatory neurons, the most plausible candidate cell types for cortex-thalamus\nconnectivity. These findings are now presented in Supplementary Fig. 9.\nLine 431 in the revised manuscript\nAlthough we did not identify trans-regional cortico-thalamic meta-gene programs, we\nextended the LIANA analysis to explore potential interactions between the cortex and\nthalamus, focusing specifically on connections between cortical deep-layer excitatory neurons\n(ExN_FrC_DL) and thalamic excitatory neurons (ExN_Tha) (Fig. 1b). Among the top 20 ligand-\nreceptor pairs in the cortex-to-thalamus direction, eight involved DEGs in either ExN_Tha or\nExN_FrC_DL (Supplementary Fig. 9), and four DEG-containing pairs were detected in the\nthalamus-to-cortex direction. These results suggest potential disruptions in cortico-thalamic\nconnectivity mediated by altered ligand-receptor interactions. However, this analysis is\nexploratory and hypothesis-generating; further studies are required to more definitively\ncharacterize these connections.\nSupplementary Figure 9. Thalamus-cortex communication estimated by LIANA.\nPredicted ligand-receptor interactions between thalamic excitatory neurons (ExN_Tha) and\ncortical deep-layer excitatory neurons (ExN_FrC_DL). Dot plots represent expression\nmagnitudes and interaction specificities, with interactions labeled as [ligand gene] -> [receptor\ngene]. Blue or red stars indicate ligand-receptor pairs in which the ligand or receptor is a DEG.\nTop: ligands expressed in ExN_FrC_DL targeting receptors in ExN_Tha. Bottom: ligands\nexpressed in ExN_Tha targeting receptors in ExN_FrC_DL.\n8. Such snRNA-seq studies usually require validation of some findings by single-molecular FISH.\nWe propose selecting the most interesting finding in gene expression and validate it by\nRNAscope or similar technologies.\nfluorescence in situ hybridization (smFISH). In response, we performed RNAscope Multiplex\nFluorescent in situ hybridization to validate one of our most notable gene expression\nobservations—the unexpected expression of SYNDIG1 in microglia. Using probes for SYNDIG1\nand CX3CR1 (a microglial marker), we confirmed SYNDIG1 expression in CX3CR1-positive cells\nwithin both the thalamus and cortex in representative specimens (new Fig. 9). We further\nmapped the spatial distribution of SYNDIG1-expressing microglia in these regions and\nquantified SYNDIG1 puncta within CX3CR1-positive cells in representative cases and controls.\nCX3CR1/SYNDIG1 double-positive cells were primarily localized to the gray matter of the\nfrontal cortex, consistent with our hypothesis—derived from geneNMF analysis (Fig. 5b and\nFig. 8e)—that SYNDIG1-positive microglia might interact with excitatory neurons. In the\nthalamus, SYNDIG1-positive microglia appeared enriched in ventricular regions, although this\ntrend was not definitive. Quantification of SYNDIG1 expression in ten specimens revealed a\ntendency toward reduced mean expression in BD compared with controls. Although sample\nsize was limited due to the scarcity of available specimens (many of which had been used),\nthese results support and extend our snRNA-seq findings, providing orthogonal single-cell\nresolution validation.\nWhile SNCA expression in PVT neurons represents another noteworthy observation relevant\nto the basic understanding of the human brain, this has already been validated by Schulmann\net al. (bioRxiv, doi: https://doi.org/10.1101/2024.09.03.611112). We therefore did not\nreplicate this experiment in our study.\nLine 461 in the revised manuscript\nSingle-molecule fluorescence in situ hybridization for SYNDIG1 in microglia\nThe unexpected detection of SYNDIG1 in thalamic and cortical microglia prompted further\nvalidation of its expression and spatial distribution using RNAscope. SYNDIG1 transcripts were\ndetected in CX3CR1-positive cells, consistent with a microglial identity (Fig. 9a).\nCX3CR1/SYNDIG1 double-positive cells—representing SYNDIG1-expressing microglia—were\npreferentially localized to the gray matter of the frontal cortex (Fig. 9b). In the thalamus, these\ncells appeared more abundant in medial regions, although this pattern was not definitive.\nComparison of representative thalamic male samples from five BD cases and five controls\nshowed relatively lower SYNDIG1 expression in BD microglia (raw FC = 0.801, Fig. 9c). While\ngeneralized linear mixed model (GLMM) analysis with a negative binomial distribution did not\nreach statistical significance, it indicated a trend toward reduction (GLMM P = 0.153, estimated\nFC = 0.783).\nLine 1170 in the revised manuscript\nSingle-molecule fluorescence in situ hybridization (smFISH)\nWe performed smFISH using the RNAscope Multiplex Fluorescent Reagent Kit v2 (Advanced\nCell Diagnostics [ACD], Newark, USA) to detect SYNDIG1 and CX3CR1 transcripts in fresh-\nfrozen human frontal cortex (FrC) and thalamus (Tha) sections from CT11 as well as additional\nrepresentative thalamic sections from five BD male cases and five male controls. SYNDIG1\nserved as the primary target, while CX3CR1 was used to identify microglia. Fresh-frozen 12-\n43\nμm-thick sections were cut onto SuperFrost Plus Gold slides and post-fixed overnight in 4%\nformaldehyde prepared in 0.1 M PBS at 4 °C. The sections were washed with 1× PBS, and\ndehydrated through graded ethanol series (30% EtOH in 1× PBS, 50% EtOH in 0.5× PBS, 70%\nEtOH in water, and 100% EtOH twice), baked at 60 °C for 30 min, and circumscribed with a\nhydrophobic barrier. Sections were treated with hydrogen peroxide and Protease IV, then\nhybridized with prewarmed RNAscope probes complementary to SYNDIG1 (C3; #804391-C3 /\nACD) and CX3CR1 (C1; #411251 / ACD) transcripts for 2 h at 40 °C. After washing, sections were\nphotobleached in 5× SSC solution for 60 min by using the autofluorescence Quenching\nilluminator TiYO (Nepa Gene, Chiba, Japan)104. Signal was amplified via AMP1–AMP3\nhybridization, and sequentially developed per probe by HRP-C#/fluorophore-tyramide (TSA\nVivid, 1:750) reactions and HRP blocking. Sections were counterstained with DAPI, mounted in\nProLong Gold Antifade, and cured at RT for 48 h. Imaging was performed on a NanoZoomer\n2.0-RS (Hamamatsu Photonics) after brightfield and fluorescence calibration each day. While\ncorrecting a focal plane (z = 0 μm) every 12 min by autofocus on brightfield, wide-field (~ 8 mm\nx 12 mm) images were captured at a high-resolution (227 nm/pixel) as multi-channel (4\nchannels; brightfield, Cy5 for TSA Vivid 650, TRITC for TSA Vivid 570, and DAPI) z-stacks (3\nslices; 0 ± 1.5 μm) in NDPI/NDPIS format. The signal quality of smFISH on RNAscope was\nconfirmed per specimen with positive (#320861 / ACD) and negative (#320871 / ACD) control\nprobes. CX3CR1 and SYNDIG1 probes were visualized with TSA Vivid 570 and TSA Vivid 650,\nrespectively. For case-control sets, one case and one control were processed in parallel to\nminimize batch effects.\nImage preprocessing was performed on multi-channel z-stacks using Fiji/ImageJ105\n(v2.16.0/1.54p). To prevent out-of-memory errors, each multi-channel z-stack was split into\nfour z-stacks per channel without importing the entire image, by editing NDPIS files. NDPIS\nfiles are text-formatted lists that specify the order in which individual channel images (NDPI\nfiles) are read and reconstructed into a multi-channel z-stack. The brightfield z-stack was\nexcluded from subsequent analyses. The remaining three fluorescent z-stacks (monochrome,\n8-bit grayscale) were sequentially imported via the Bio-Formats Importer, underwent\nbackground subtraction by the sliding paraboloid algorithm (radius: 300 pixels, disable\nsmoothing: true), converted into maximum intensity z-projection images (MIzPs), and saved\nto disk to free memory. The fluorescent MIzPs were then assembled into a multi-channel stack\nwith pseudo-color.\nObject detection and quantification were performed by using QuPath106 (v0.5.1). A\nconsistent fluorescent signal threshold of 45 in 8-bit (0–255) grayscale was applied to each\nchannel across all specimens. Cells were detected using the WatershedCellDetection plugin,\ndefined as areas containing a nucleus (DAPI) and a peripheral region within 3-μm of the\nnucleus. By the SubcellularDetection plugin, the smFISH signals of SYNDIG1 (TSA Vivid 650, Cy5\nchannel) and CX3CR1 (TSA Vivid 570, TRITC channel) were detected as puncta or clusters in\ncells, and counted as the number of spots estimated, assuming that single spots are\napproximately 1 μm2 (range: 0.5 μm2 to 3.0 μm2). Additional details are available at\nhttps://github.com/msk240/BDsnRNAseq\nData aggregation was conducted in R (v4.1.0). Positive cells were defined as those containing\nthree or more spots corresponding to the target gene, based on the number of spots estimated\nfrom smFISH signals. SYNDIG1 expression in CX3CR1-positive cells, representing putative\nmicroglia, was assessed using a generalized linear mixed model (GLMM) implemented in the\nglmmTNB package (v1.1.12), with diagnosis (Dx) and age as fixed effects and sample ID as a\nrandom effect, employing the nbinom2 parameterization of the negative binomial distribution.\nSex was not included as a covariate, as all samples were from male donors.\nFigure 9. Single-molecule fluorescence in situ hybridization of SYNDIG1 transcript\na) RNAscope validation of SYNDIG1 expression in CX3CR1-positive cells within the thalamus\nand cortex. Left panels: Low-magnification view of the thalamus and cortex with DAPI staining.\nWhite rectangles denote the regions magnified in the middle panels. Middle panels: Magnified\nviews highlighting CX3CR1-positive cells (likely microglia), indicated by white arrowheads.\nMagenta signals represent SYNDIG1 expression, which varies among CX3CR1-positive cells.\nThe white square indicates the region further magnified in the right panels. Right panels:\nHigher magnification of CX3CR1/SYNDIG1 double-positive cells, corresponding to\nCX3CR1(+)/SYNDIG1(+) cells in panel b. b) Spatial distribution of cell types classified by\nCX3CR1/SYNDIG1 expression. CX3CR1(+)/SYNDIG1(+) cells are considered SYNDIG1-expressing\nmicroglia. CX3CR1(+)/SYNDIG1(-) cells represent microglia with little or no SYNDIG1\nexpression. CX3CR1(-)/SYNDIG1(+) cells likely correspond to non-microglial SYNDIG1-\nexpressing cells, such as excitatory neurons. CX3CR1(-)/SYNDIG1(-) cells are classified as other\ncell types. Cell counts per hexagonal bin are assigned a color code at log2 scale. White contours\nare the median of cell density for each cell classification. c) Quantification of SYNDIG1 puncta\nin CX3CR1-positive cells (likely microglia) across 10 representative male samples (five BD cases\n45\nand five controls). Each dot represents an individual cell.\n9. The cortex does not show large changes in gene expression that is why the authors argue to\nfocus on the thalamus. Nevertheless, this is surprising given previous role of the cortex that has\nbeen proposed in brain dysfunction in BD. There is a lot of focus in the Discussion section on\nthe thalamus, but such surprisingly minor changes in the cortex are not discussed. We think\nthese findings are very important, since the study can pave the way for future analyses of\npostmortem brain tissue in psychiatric disorders. Currently, almost all studies focus on the\nprefrontal cortex and maybe in the future the thalamus should be the major region for such\nsnRNA-seq analyses.\nWe thank the reviewer for highlighting this important point. In the revised Discussion, we have\nadded a dedicated paragraph addressing the unexpectedly subtle changes observed in the\ncortex, in contrast to the well-established cortical emphasis in prior BD studies. While we\nacknowledge the significance of cortical pathology, we also note that some degree of\n“streetlight effect” (i.e., the tendency to investigate areas that are more accessible or better\ncharacterized) may have contributed to the predominant focus on the cortex in previous\nresearch.\nOur findings suggest that the thalamus may play a more central role in the pathophysiology\nof BD than previously recognized. This is especially notable given the relative inaccessibility of\nsubcortical regions in human postmortem research. We believe that this shift in focus could\nhelp guide future transcriptomic and neuroimaging investigations in psychiatric disorders,\npotentially revealing underexplored yet biologically critical brain regions.\nLine 572 in the revised manuscript\nAlthough BD is also associated with reduced cortical volume, the transcriptional and\ncompositional changes in the cortex were less pronounced than those in the thalamus.\nThe\nfrontal cortex has been frequently analyzed in psychiatric research, not only due to its\nsuspected involvement in disease pathology but also because of its relative accessibility. This\nmay reflect a form of the “streetlight effect,” wherein research tends to focus on more\naccessible regions. However, higher accessibility does not necessarily equate to higher\nbiological significance. Our findings underscore the importance of extending investigations to\nsubcortical regions, which may harbor critical yet underexplored components of BD\npathophysiology. It is noteworthy that the cortical region analyzed in this study was BA10.\nResults may differ in other cortical regions, such as the anterior cingulate cortex, emphasizing\nthe need for broader sampling across diverse cortical areas. Integration of GWAS data with\nbrain-wide snRNA-seq indicated heritability enrichment in cortical neuronal populations,\nincluding interneurons and deep-layer intra-telencephalic neurons1. While such genetic\ninsights are informative, caution is warranted when inferring cell-type-specific dysfunction\nbased solely on heritability enrichment. Genetic enrichment does not necessarily guarantee\ntranscriptional dysregulation. Moreover, technical variability in defining cell-type-specific gene\nexpression profiles can influence enrichment outcomes and lead to divergent\ninterpretations78. Therefore, direct observations of transcriptional and compositional changes\nremain essential to elucidate disease-associated cellular alterations.\n1. Discussion lacks similarity of BD and schizophrenia, although CACNA1C is one of the major\ngenes in the current study and it has been identified as genetic risk factor by schizophrenia\nGWAS and not BD.\nWe thank the reviewer for this important point. In the revised Discussion, we now highlight\nthe shared genetic architecture between bipolar disorder (BD) and schizophrenia, with\nparticular emphasis on CACNA1C, which has been consistently implicated in both disorders.\nAlthough CACNA1C is not included among the 36 causal genes for BD listed in the PGC4 study\n(O’Connell et al., Nature 2025), it has remained one of the top loci associated with BD across\nmultiple GWAS phases based on p-value rankings:\nPGC1: ranked 1st in Sklar et al., Nat Genet 2011 (Table 3)\nPGC2: ranked 4th in Stahl et al., Nat Genet 2019 (Table S4)\nPGC3: ranked 3rd in Mullins et al., Nat Genet 2021 (Table S4)\nPGC4: ranked 6th in O’Connell et al., Nature 2025 (based on our in-house MAGMA analysis)\nThe exclusion of CACNA1C from the 36 causal gene list in PGC4 may reflect methodological\nfactors. Causal gene identification remains largely inferential and is influenced by analytical\nparameters and reference datasets, rather than definitive experimental validation. Notably,\nour own analyses using fine-mapping and TWAS integration based on previous BD GWAS—an\napproach comparable to that used in the PGC4 study—also identified CACNA1C as one of\npromising BD-associated genes. This convergence might support the hypothesis that BD and\nschizophrenia share molecular mechanisms, including dysregulation of calcium channel\nsignaling, particularly within thalamic excitatory neurons. We have expanded the revised\nDiscussion to address this potential overlap in pathophysiology.\nWe have also added the relative expression profiles of CACNA1C alongside SHISA9 to highlight\nthese genes as probable DEGs in PVT neurons (Fig. 8). Both genes show the highest expression\nin thalamic excitatory neurons within the thalamus-cortex integration set.\nLine 227 in the revised manuscript\nWhile CACNA1C is not listed in these 36 credible genes, CACNA1C has consistently ranked\namong the top genes by p-value in prior PGC1-4 studies1,2,50,51, and was also implicated by our\nFUSION-based TWAS integration.\nLine 489 in the revised manuscript\nAmong all the thalamic and cortical cell types in our study, PVT neurons exhibited the highest\nnumber of DEGs in BD. Notably, downregulated DEGs in PVT neurons were more biologically\ninterpretable than upregulated ones and included putative BD risk genes. Among the PVT\nneuron DEGs, SHISA9, CACNA1C, and KCNQ3 stood out as particularly relevant to BD\npathophysiology, as they were centrally positioned in the PPI network and consistently\nhighlighted in GO, GSEA, and SynGO analyses. SHISA9 is one of the 36 credible genes identified\nin the largest BD GWAS. Although CACNA1C is not among these 36, it has repeatedly ranked\namong the top BD-associated loci by p-value across PGC1-4 studies1,2,50,51 and was also\n47\nimplicated in our FUSION-based TWAS analysis. CACNA1C has additionally been associated\nwith schizophrenia in prior GWAS33, suggesting a shared genetic architecture between BD and\nschizophrenia. Consistently, downregulated DEGs in PVT neurons were enriched for\nschizophrenia polygenic risk as assessed by MAGMA. While KCNQ3 is not classified as a\nputative BD risk gene, its downregulation was reported in bulk brain tissue from BD patients61.\nKCNQ3 may represent a BD-associated gene identified through transcriptomic rather than\ngenomic approaches.\nFigure 8. Comparative analysis between the thalamus and cortex.\nc) Relative expression of CACNA1C and SHISA9 across all clusters in the thalamus-cortex\nintegration dataset.\n2. In general, although Discussion is quite long, some points are extensively discussed, whereas\nothers are missing. For instance, close relation of thalamic excitatory neurons with GABAergic\nneurons is interesting. Is it related to cell fate or neuronal specialization? Similarly, cardiac\ntissue similarity can be moved from Results to Discussion section.\nWe appreciate the reviewer’s insightful comment. In response, we have expanded the\nDiscussion to address the transcriptional similarity between thalamic excitatory neurons and\ncortical inhibitory neurons. We now suggest that this resemblance may reflect a unique\nsubtype specialization potentially related to electrophysiological function, rather than\nindicating a shared developmental origin. In line with the reviewer’s suggestion, we have also\nrelocated the discussion of cardiac gene expression profiles in thalamic neurons from the\nResults section to the Discussion section.\nLine 72 in the original manuscript\nREMOVED: When compared to ExN_FrC, ExN_Tha showed a preferential expression of genes\nassociated with cardiac muscle development, such as SHOX2, PROX1, HCN4, and CACNA1G, as\ndetermined by gene set enrichment analysis (GSEA) (Fig. 1e, Supplementary Table 2,\nSupplementary Data 1). The unexpected enrichment of genes associated with cardiac muscle,\n48\nparticularly atrial muscle cells, may reflect a shared characteristic of frequent and rhythmic\nCa²⁺ regulation between thalamic neurons and cardiac muscle. CACNA1G, which encodes a\nsubunit of T-type calcium channels involved in electrical pace-making activity, is specifically\nexpressed in the thalamus and the heart30, and has been identified as a core gene associated\nwith schizophrenia31. Among the ExN_Tha-specific genes, CASQ2 was highly characteristic to\nExN_Tha, with virtually no expression in other cell clusters (Supplementary Fig. 2a).\nCASQ2\nencodes calsequestrin-2, a high-capacity calcium-binding protein particularly expressed in\ncardiac muscle cells32, further underscoring the commonality between thalamic excitatory\nneurons and cardiac muscle cells, probably reflecting the frequent and rhythmic activities with\nintra-cellular Ca²⁺ changes.\nLine 102 in the original manuscript\nREMOVED: These findings may suggest shared characteristics of frequent and rhythmic\nactivities with intra-cellular Ca²⁺ changes between thalamic neurons and cardiac muscle cells.\nLine 543 in the revised manuscript\nIn addition to PVT neurons, other thalamic excitatory neurons express CALB1, CALB2,\nand/or\nPVALB, all of which encode calcium-binding proteins. The widespread expression of these\nproteins may underlie the unexpected transcriptomic similarity observed between thalamic\nexcitatory neurons and cortical inhibitory neurons. Calcium-binding proteins are critical for\nbuffering rapid fluctuations in intracellular calcium concentrations associated with high-\nfrequency neuronal firing. Kim et al. reported that thalamic excitatory neurons originate from\nradial glia derived from prosomere 2 of the embryonic diencephalon35, whereas cortical\ninhibitory neurons arise from the ganglionic eminence, highlighting their distinct\ndevelopmental origins. The observed transcriptomic similarity is unlikely to reflect shared\nembryonic origins but may instead suggest functional convergence, particularly in intracellular\ncalcium regulation. Interestingly, both excitatory and inhibitory thalamic neurons showed\npreferential expression of genes typically associated with cardiac muscle. This observation\nmay reflect a shared physiological trait—frequent and rhythmic intracellular calcium\nregulation—between thalamic neurons and cardiac myocytes.\n3. Lines 137-141 compare data from the current study to previously published organoid data.\nThis comparison is completely separated from all other analyses in the study and in the current\nstate should be removed from the manuscript.\nWe agree that this comparison was not sufficiently integrated with the main narrative.\nAccordingly, we have removed this section, including supplementary figures, from the\nmanuscript to maintain focus and cohesion.\nLine 137 in the original manuscript\nREMOVED: To explore opportunities for simulating human PVT neurons in vitro, we\ninvestigated whether human thalamic organoids contain PVT-like neuronal populations. PVT-\nlike neurons were identified in the thalamic organoids created by Kiral et al.41 (Supplementary\nFig. 2d-e). These neurons expressed PVT marker genes, including SLC17A6, CALB2, and SNCA.\nTheir expression profile suggests the feasibility of modeling human PVT neurons in vitro.\nREMOVED: Despite these limitations, this study indicates the disturbances of PVT neurons as\na core pathology in BD, providing a promising foundation for developing new diagnostic and\ntherapeutic strategies. Neuroimaging techniques targeting the PVT, such as MRI and PET,\ncould enhance clinical assessment for BD. PVT-like neurons in thalamic organoids underscore\nthe feasibility of modeling PVT neurons from patient-derived peripheral samples, which can\naccelerate the discovery of novel pharmacological interventions to the molecular mechanisms\nunderlying BD. Integrating multi-modal methodologies, encompassing cellular, molecular, and\nneuroimaging techniques, will provide a more comprehensive understanding of the\npathological contributions of PVT in BD. PVT neurons represent a promising target for the\ndevelopment of advanced diagnostic and therapeutic strategies aimed at improving patient\noutcomes.\nLine 833 in the original manuscript\nREMOVED: Kiral et al. data can be obtained from GEO with accession number GSE210720.\nSupplementary Figure\nREMOVED: Supplementary Figure S2d-e in the original manuscript\n4. In Figure 3b, it will help to add direction of \"increased/decreased\" for BD.\nWe thank the reviewer for this helpful suggestion. We have now added directional annotations\nin Figs. 3b/5b and Supplementary Fig.15a to clarify compositional changes for readers.\nFigure 3b\nSupplementary Figure 15a\n51\n5. We are not sure why the authors decided to include reference to NeuCyto in Lines 158-159.\nIf this “cell type” is ambient RNA contamination, then it is difficult to make conclusions based\non such “cell type”.\nWe appreciate this observation. The NeuCyto cluster represents cytoplasmic components\nderived from neuronal populations and was used solely as a technical reference. Due to their\nindirectness, it was excluded from quantitative analyses. We have revised the text to clarify\nthis point and prevent any misinterpretation.\nLine 157 in the original manuscript\nREMOVED: For reference, a decrease in Tha_NeuCyto_VGLUT2 (putative excitatory neuron\ncytoplasm, see Methods) indirectly supported these reductions (Supplementary Fig. 3c).\nNonetheless, the presence of such cluster would be interesting for snRNA-seq communities.\nThus, we have revised the method sections and Supplementary Fig. 13 describing NeuCyte for\nclarity.\nLine 789 in the revised manuscript\nClusters designated as NeuCyto were considered neuronal cytoplasmic components from the\nfollowing properties: low nuclear fraction, high content of mitochondrial RNA (Supplementary\nFig. 13a), high expression of neuronal markers, virtually no expression of other glial markers\n(Supplementary Fig. 1a, 3a, and 6a), and distinctive clustering from probable empty droplets\n(Supplementary Fig. 13b-c). The probable empty droplets, designated as “Empty”, were\ndefined as droplets filtered out by Cell Ranger with detected genes of 100 – 300. This range\nwas set by QC procedure and barcode rank plot for UMI counts by Cell Ranger. NeuCyto\ndroplets were robustly detected from 40/41 of our thalamic samples (occupying 1.46% of total\ncell nuclei) and 40/41 of our cortex samples (occupying 1.46% of total cell nuclei) as well as\n20/20 of HBCA thalamic samples (occupying 1.66% of total cell nuclei). The clear separation of\nVGLUT+ and GAD+ NeuCyto clusters in the thalamus indicated that they were distinct\n52\ncytoplasmic components (Fig. 3a). We removed NeuCyto clusters for the subsequent analysis\nas non-nuclei droplets but retained for UMAP visualization and dot plots as a technical\nreference. Indeed, in the final nuclei suspension, we observed objects stained with Nonyl\nAcridine Orange (NAO) but negative for Propidium Iodide (PI) and Hoechst 33342\n(Supplementary Fig. 13d), likely representing mitochondria-rich, nucleus-free cytoplasmic\nfragments corresponding to the NeuCyto clusters.\nSupplementary Figure 13. NeuCyto clusters and additional spatial deconvolution.\na) Nuclear fractions and mitochondrial RNA fractions across samples in the thalamic\nintegration dataset (see Fig. 3a). b) UMAP projection illustrating “Empty” clusters (ambient\nRNA, 1,000 droplets per sample) and canonical cell nuclei clusters shown in Fig. 3a. c) UMAP\n53\nprojection of the same clusters as in panel b, with distinct coloring to highlight\nTha_NeuCyto_VGLUT, Tha_NeuCyto_GAD, and Empty clusters. d) Microscopic images of\ncytoplasmic fragments (white arrowheads) in two representative spots. Cytoplasmic\nfragments are indicated by NAO staining (green) in the absence of PI or Hoechst 33342 signals.\nCell nuclei appear magenta due to combined PI and Hoechst 33342 staining. e) CytoSPACE\ndeconvolution based solely on our snRNA-seq dataset (excluding the HBCA dataset). f)\nCytoSPACE deconvolution using only the CT06 snRNA-seq dataset. g) UMAP projection\ndisplaying three NeuCyto sub-clusters within the thalamus. h) Predicted spatial distribution of\nthe three NeuCyto sub-clusters. Color intensity in the side bars reflects the estimated number\nof cells per spatial spot. For panels g and h, only control samples were used to prevent\npotential bias from BD-associated alterations in spatial deconvolution.\n6. Please correct Lines 255-256:\nrelated to synaptic function, although microglia are primarily regarded as an immune-related\ncell type\nMicroglia has well known role in synaptic plasticity and modification of neuronal connectivity,\nwhich is also one of their primary roles in the brain.\nWe thank the reviewer for this helpful suggestion. In response, we have revised the relevant\nsentence (moved to the Discussion) to better reflect the established role of microglia in\nregulating synaptic plasticity and neuronal circuitry. The revised sentence strengthens our\ninterpretation and supports the overall argument.\nLine 254 in the original manuscript\nREMOVED: Furthermore, all five DEG-GOs identified within Tha_Micro_1 were related to\nsynaptic function, although microglia are primarily regarded as an immune-related cell type\n(Supplementary Table 5).\nLine 556 in the revised manuscript\nBeyond cell-type-based analysis, we identified trans-cellular meta-gene programs between\nexcitatory neurons and microglia (Excit-Micro meta-gene programs), likely reflecting\nphysiological synaptic interactions between these cell types. Microglia are known to modulate\nsynaptic plasticity and neuronal circuitry through dynamic interactions with neurons76.\n7. Figure 8a: show p-adjusted line on the figure\nWe have updated Fig. 8a to include a line indicating the Bonferroni-corrected significance\nthreshold. This improves visual interpretation of the statistical results.\nFigure 8a\n54\n8. Lines 452-460 repeat lines 463-472.\nThank you for pointing this out. We have removed the redundant text to improve clarity and\nflow.\nLine 452 in the original manuscript\nREMOVED: Despite the limitations of this study, our findings suggest that disturbances in PVT\nneurons constitute a central pathology in BD. Neuroimaging techniques targeting the PVT can\nenhance clinical assessments of BD. The presence of PVT-like neurons in thalamic organoids\nhighlights the feasibility of modeling PVT neurons from patient-derived peripheral samples.\nThis approach could facilitate the discovery of pharmacological interventions that target the\nmolecular mechanisms underlying BD. The integration of multi-modal methodologies,\nincluding cellular, molecular, and neuroimaging techniques, will offer a more comprehensive\nunderstanding of the pathological role of PVT neurons in BD. PVT neurons emerge as a\npromising target for the development of novel diagnostic and therapeutic strategies for\npatients with BD.\n9. Correct misprint:\nFig1a Neruon-microglia synaptic interaction in BD (Figure 5)\nrecognized that the original legend for Fig. 1a did not adequately explain the figure. We have\ntherefore revised the legend of Fig. 1a to improve clarity.\nFigure 1a\na) Schematic representation of the study design. Thalamic (magenta) and cortical (green)\nspecimens from 21 BD cases and 20 controls underwent snRNA-seq to identify BD-associated\ncellular changes. A dashed line indicates the medial thalamic region analyzed in this study.\nNishioka et al. present a snRNA-seq study of PVT in 21 bipolar disorder (BD) patients and 20\ncontrols. BD is a severe mental health condition marked by recurring manic and depressive\nepisodes, requiring new treatments due to resistance, adverse effects, and high suicide rates.\nRecent research has identified the thalamus, particularly the paraventricular thalamic nucleus\n(PVT), as a crucial brain region for mood regulation, with significant differences in BD patients\ncompared to controls. Authors report that thalamus samples from the BD patients they studied\nexhibited a 50% reduction in PVT neurons, which also showed notable downregulation in genes\nrelated to synaptic function and ion channel regulation. Additionally, disruptions in interactions\nbetween PVT neurons and microglia, particularly in AMPA receptor signaling pathways, were\nobserved. These findings suggest that targeting AMPA signaling in PVT neurons could pave the\nway for innovative treatments and diagnostics aimed at improving outcomes for BD patients.\nThis paper has several strengths: Well written, addresses important questions, focuses on an\nunderstudied brain region, uses mostly state-of-the-art methods.\nThe biggest concern has to do with overinterpretation of results, with too little consideration\nfor confounding and other limitations. There is also insufficient methodological detail in several\nareas.\nWe sincerely thank the reviewer for their thoughtful and constructive comments. In response,\nwe have carefully revised several sections to avoid overinterpretation and to better reflect the\nlimitations of our data. Specifically, we have removed or revised speculative statements\nconcerning SYNDIG1 and SHISA9 proteins.\nThe Methods section has also been expanded to explicitly address potential confounders,\nincluding disease status (Dx), age, sex, postmortem interval (PMI), and pH. Additionally, we\nnow provide a clear rationale for excluding certain other potential confounders, such as\nexperimental/sequencing batch, substances at death, race, and cause of death, based on their\nminimal impact or the risk of introducing bias due to incomplete metadata. We have further\nrevised the Methods to enhance clarity and reproducibility, particularly in describing our\nthalamic sampling strategy.\nThe reviewer’s comment led us to identify and correct an error in our integration of the\nHuman\nBrain Cell Atlas (HBCA) dataset (Siletti et al., Science 2023). We had inadvertently included\nhypothalamic samples instead of the intended epithalamic ones. Specifically:\n- Incorrect samples used:\n10X190_3, 10X190_4, 10X193_7, 10X193_8, 10X348_3, 10X348_4, 10X354_5, 10X354_6,\n10X361_3, 10X361_4, 10X362_1, 10X362_2, 10X375_3, 10X375_4, 10X377_1, 10X377_2,\n10X389-1, 10X389-2, 10X393_1, 10X393_2\n- Correct samples:\n10X190_3, 10X190_4, 10X193_7, 10X193_8, 10X348_3, 10X348_4, 10X354_5, 10X354_6,\n10X361_3, 10X361_4, 10X362_1, 10X362_2, 10X375_3, 10X375_4, 10X377_1, 10X377_2,\n10X389-5, 10X389-6, 10X393_1, 10X393_2\nThe samples 10X389-1 and 10X389-2 are derived from hypothalamus tissue, whereas 10X389-\n5 and 10X389-6 are derived from the epithalamus, as originally intended. We apologize for this\n57\noversight.\nTo correct this, we re-performed the integration using the correct epithalamus-derived\nsamples, followed by the full downstream analyses. The original integration, which included\nhypothalamic samples, was enriched in Splatter neuron clusters due to the nature of\nhypothalamic tissues, as illustrated in Siletti et al. 2023. As expected, the revised integration\nshowed a reduction in Splatter neuron clusters and yielded clearer cluster distinction.\nMoreover, we adjusted the clustering resolution accordingly, since the previous resolution was\noverly high. Given that our aim is to characterize PVT neurons, we optimized the resolution to\nbetter distinguish PVT neurons from other thalamic excitatory populations.\nSubsequently, we revised the downstream analyses based on the corrected HBCA-integrated\ndataset. As a result, the characteristics of CytoSPACE-based deconvolution, the\ncorrespondence to mouse PVT neurons, and the gene set enrichment analysis (GSEA) of PVT-\nneuron-specific genes became clearer and more interpretable, as detailed in the revised\nsections. We thank the reviewer for helping us identify and address this issue, which improved\nthe robustness and clarity of our dataset and analyses. Notably, the excluded datasets\n(10X389-1 and 10X389-2) were not critical to our main analyses and appeared to reflect off-\ntarget contamination. Their exclusion does not affect the validity of the original results.\nLine 733 in the revised manuscript\nThe original sample IDs were 10X190_3, 10X190_4, 10X193_7, 10X193_8, 10X348_3,\n10X348_4, 10X354_5, 10X354_6, 10X361_3, 10X361_4, 10X362_1, 10X362_2, 10X375_3,\n10X375_4, 10X377_1, 10X377_2, 10X389-5, 10X389-6, 10X393_1, and 10X393_2, referred to\nas HBCA01-20 in our study (Supplementary Data 11).\na) UMAP projection displaying cell clusters identified in the HBCA-integrated thalamic dataset.\n88.6% of ExN_Tha_CALB cells in Fig. 1b are included within the ExN_CALB cluster here. b) Top:\nRepresentative paraventricular thalamus with a dashed black rectangle indicating the target\nregion for Visum analysis. Bottom: Illustrated version with topographical orientation. D, dorsal.\nV, ventral. M, medial. L, lateral. c) Left: Spatial transcriptomic data delineating 15 distinct\nclusters, each represented by a different color, illustrating the spatial organization of gene\nexpression across the paraventricular thalamus. Right: Calretinin immunohistochemistry\nstaining of the adjacent section of Visium analysis. d) Spatial localization of cell clusters as\npredicted by CytoSPACE deconvolution, illustrating their estimated anatomical distribution\nwithin the thalamus. Color intensity in the sidebars indicates the estimated number of cells\nper spot, providing a quantitative measure of cell density across regions. e) Merged view of\nthe spatial clusters from panel d for comparison with panel c. f) Matching scores of each cell\ncluster to mouse PVT neurons. g) Proportions of mouse PVT neuron subtypes from Gao et al.\nand Shima et al. corresponding to human ExN_CALB clusters. All PVT subtypes refer to mouse\nannotated PVT subtypes in pPVT, posterior PVT. AL, anterior lateral. AM, anterior medial.PL,\nposterior lateral. PM, posterior medial. PVT# in Gao et al. denotes the original subtype\nclassification. h) GO network visualization for PVT-enriched genes, illustrating major functional\npathways related to dopamine metabolism and calcium/ion channels. GO terms with\nenrichment scores > 0.5 and network connectivity > 2 are shown; core genes are highlighted.\nCore genes involved in these pathways are highlighted. i) Spatial expression profiles of key PVT\nmarker genes, with darker red indicating higher expression. PVALB serves as a negative control\nto demonstrate spatial specificity.\nConcerns in order of importance:\n• Anatomical definition of PVT\n• CALB1 and CALB2 were used to define PVT, but they are not exclusive to PVT; CALB1 is broadly\nexpressed in all ‘matrix’ nuclei; CALB2 is expressed in all midline and all intralaminar nuclei\n• Needs more detailed description (ideally with illustration) how the “medial thalamus” was\ndissected; e.g. more anterior part bordering anterior thalamic nuclei or more posterior part\nbordering mediodorsal thalamus\nWe agree with the reviewer that CALB1 and CALB2 are not exclusive markers of the\nparaventricular thalamic nucleus (PVT). In the revised manuscript, we have clarified that these\ngenes were initially used as putative markers for PVT neurons. This preliminary identification\nwas subsequently refined through anatomical localization using spatial transcriptomic data\nand anchor-based transfer from mouse PVT neuron data. Our anatomical definition of the PVT\nfollows Uroz et al. 2004.\nOur workflow for identifying human PVT neurons proceeded as follows. First, excitatory\nneuron clusters co-expressing CALB1 and CALB2 were identified as putative PVT neuron\ncandidates. These were then validated by (i) spatial localization using CytoSPACE\ndeconvolution and (ii) transcriptomic correspondence to mouse PVT neurons. While the\nmouse PVT is anatomically well-defined, the human PVT lacks a sharp boundary with\nsurrounding regions. This anatomical ambiguity is also reflected in transcriptomic space,\nwhere PVT clusters form a continuum with adjacent excitatory populations, consistent with\nobservations in Schulmann et al. (bioRxiv 2024). We have revised the text to more clearly\nreflect this stepwise identification process.\nTo further clarify our dissection protocol, we have added a new Supplementary Figure\n(Supplementary Fig. 10b-d) illustrating how blocks were cut from coronal thalamic slabs\ncontaining the mammillary body, targeting the mid-PVT. A more detailed description of the\ndissection strategy has also been incorporated into the Methods section. In addition, we have\nattached a demonstration video illustrating the dissection of slabs containing the PVT from the\nhuman postmortem brain, narrated by Tadafumi Kato. This video was made for the brain\ncutting in this project in 2013, as an internal communication between Japan and Canada. As\nwe have no copyright for the photo in this video, this video is not to be public. We attach this\nvideo only for review purposes.\nThis description was derived from the CT06 section for Visium analysis. We aimed at the mid-\nPVT and inherent anatomical variation resulted in slight differences in sampling positions along\nthe anterior-posterior axis. As anatomical studies of the human PVT remain limited, we hope\nthat our findings, including candidate molecular markers, may contribute to future efforts to\nbetter characterize this region.\nLine 483 in the original manuscript\nREMOVED: We selected slabs including a relatively anterior side of PVT and cut the specimen\nwith enough margins from PVT in the coronal section (Fig. 1a).\nLine 102 in the revised manuscript\nCytoarchitecture of PVT neurons\nAmong the thalamic excitatory neurons, ExN_Tha_CALB cells in Fig. 1b were putatively\nsupposed to be PVT neurons, as characterized by the expression of CALB1 and CALB2—\nmarkers for PVT neurons in humans36 and rodents37,38—while these two markers were also\nexpressed in other clusters. These putative PVT neurons were also robustly identified in the\nHuman Brain Cell Atlas (HBCA) dataset6, as demonstrated by the integration of our and HBCA\ndatasets (Fig. 2a). ExN_CALB in Fig. 2a corresponded to ExN_Tha_CALB in Fig. 1b. The putative\nPVT neurons (ExN_CALB cells) were detected from all 20 our samples and 20 HBCA samples\n(see Methods). This integrated dataset was used for subsequent analyses to enhance\ngeneralizability.\nAs CALB1 and CALB2 expression alone was insufficient to confirm that the ExN_CALB\nclusters represented PVT neurons, we validated this classification using two approaches: (1)\nspatial transcriptomic deconvolution and (2) anchor-based label transfer using mouse PVT\nneuron reference data. We obtained spatial transcriptomic data of the medial thalamus from\none control sample as a representative specimen (Fig. 2b). The PVT was characterized through\nunbiased clustering of spatial transcriptomic data and partially corresponded to the areas\nimmunostained by calretinin, a protein coded by CALB236 (Fig. 2c). Deconvolution of thalamic\nspatial transcriptomic data using the integrated cell clusters by CytoSPACE39 indicated that\nExN_CALB cells were predominantly localized in the paraventricular regions of the thalamus\nwith a clear demarcation from other excitatory neuron populations (Fig. 2d-e), aligning with\nthe anatomical definition of the PVT by Uroz et al.36. Oligodendrocytes were concentrated in\nthe lateral side of the thalamus, while astrocytes were primarily found in the medial side.\nOligodendrocyte precursor cells (OPCs) and microglia displayed a diffuse distribution\nthroughout the thalamus. Expectedly, ExN_CALB corresponded most closely with mouse PVT\nneurons when employing anchor transfer with published single-cell datasets from mice23,40\n(Fig. 2f). These data confirmed the validity of ExN_CALB as PVT neurons. Upon closer\nexamination of PVT neuron subtypes, the majority of human PVT neurons, particularly those\nfrom our dataset, corresponded to mouse posterior PVT neurons (Fig. 2g).\nLine 649 in the revised manuscript\nBulk specimens of the medial thalamus and frontal cortex (BA10) were dissected from archived\nfresh-frozen slabs and then sectioned at a thickness of 200 μm using a CM3050S cryostat-\nmicrotome (Leica Microsystems, Wetzlar, Germany). First, the whole brain was bisected into\nleft and right hemispheres. The left hemisphere was snap-frozen in 2-methylbutane precooled\n61\nto -40 °C and stored at -80 °C, while the right hemisphere was fixed in 10% neutral buffered\nformalin (phosphate buffer, pH 7.0; Fisher Scientific, Cat. #SF-100-20) for paraffin embedding\nand sectioning. Within a half hemisphere, the PVT was macroscopically identified, with its\nanterior boundary defined by the coronal plane at the anterior edge of the mammillary body\nand its posterior boundary by the coronal plane at the anterior edge of the posterior\ncommissure, spanning approximately 20 mm along the anterior-posterior axis. A 10-mm-thick\ncoronal slab encompassing the PVT was then selected and excised. Each slab extended\nanteriorly to the center of the mammillary body and posteriorly to midbrain structures\nincluding the red nucleus and substantia nigra. From each slab, tissue blocks were prepared to\ninclude the PVT together with major adjacent landmarks (putamen, caudate, and insular\ncortex), ensuring accurate anatomical orientation and precise localization of the PVT (Fig. 1a,\nSupplementary Fig. 10b-d). Dissections were typically performed from the posterior aspect of\nthe slab, with each block carefully inspected to confirm that the slab contained the appropriate\nrostrocaudal level of the PVT. Although the procedure was designed to target the mid-PVT,\nunavoidable interindividual variation in gross neuroanatomy occasionally resulted in slight\nshifts in sampling position along the anterior-posterior axis. The medial region of the thalamic\ncryosections and the gray matter part of the cortical cryosections were excised for snRNA-seq.\nApproximately 50 mg of cryosections, corresponding to two to four cryosections, underwent\nthe subsequent cell nuclei isolation.\nb) Schematic illustration of the targeted thalamic region in the left hemisphere. The 10-mm-\nthick slab selected for PVT dissection is indicated by violet lines. The center of the mammillary\nbody and the red nucleus/substantia nigra, used as anatomical landmarks, are indicated by red\nand green circles, respectively. c) Representative coronal brain slab encompassing the medial\nthalamus (black dashed rectangle), the region of interest in this study. d) Representative\nthalamic tissue block for cryosectioning.\nReview material\nDissection_PVT.wmv (video):\nDemonstration video illustrating the dissection of slabs containing the PVT from the human\npostmortem brain.\n63\n• Gao et al. and Shima et al. have identified major differences within PVT, especially along\nanterior-posterior axis; the author only included a general scoring of cells relative to these\nmouse datasets, not specific to PVT subdivisions\nWe thank the reviewer for this insightful suggestion. Our anatomical dissection did not\nexplicitly distinguish between anterior and posterior subregions in the human PVT due to\nlimitations in tissue availability. In response to this comment, we performed a new analysis by\nmapping our identified ExN_CALB cluster (human PVT neurons) onto the Gao et al. and Shima\net al. mouse datasets using an anchor-based integration approach. We used the corrected\nHBCA-integrated dataset for this new analysis.\nThis analysis revealed that the ExN_CALB population in our dataset aligns more closely with\nposterior PVT (pPVT) neurons in both reference datasets. Notably, despite methodological\ndifferences (Gao et al. used single-nucleus RNA-seq with the 10x Chromium platform, whereas\nShima et al. employed single-cell RNA-seq with Quartz-Seq2), the observed correspondence\nwas consistent across datasets. These findings have been added to the Results and Methods\nsections, and the mapping results are presented in Fig. 2g.\nLine 123 in the revised manuscript\nExpectedly, ExN_CALB corresponded most closely with mouse PVT neurons when employing\nanchor transfer with published single-cell datasets from mice23,40 (Fig. 2f). These data\nconfirmed the validity of ExN_CALB as PVT neurons. Upon closer examination of PVT neuron\nsubtypes, the majority of human PVT neurons, particularly those from our dataset,\ncorresponded to mouse posterior PVT neurons (Fig. 2g).\nLine 513 in the revised manuscript\nIn mice, recent work has revealed functional and transcriptomic differences between the\nanterior and posterior PVT23,38,40. In our study, most human PVT neurons corresponded to the\nposterior PVT neurons in mice. Future investigations examining the entire PVT will be essential\nto clarify the evolutionary correspondence of PVT neurons between humans and mice.\nLine 869 in the revised manuscript\nMatching mouse PVT neurons to the human PVT neurons\nTo identify cell types in our human snRNA-seq data that correspond most closely to mouse\nPVT neurons, we performed anchor-based label transfer using the Seurat functions\nFindTransferAnchors and TransferData. As reference datasets of mouse PVT neurons, we used\nthe snRNA-seq data from Gao et al. (2023)40 and the single-cell RNA-seq data from Shima et\nal. (2023)23. The mouse gene name was converted to the corresponding human gene name in\nthe Ensemble database. Prior to anchor transfer, each dataset was independently processed\nusing Seurat, including normalization, scaling, and principal component analysis (PCA). For the\nGao et al. dataset, we followed the authors’ published analytical pipeline to identify PVT\nsubtypes, with the exception of PCA dimensions to improve clustering resolution, and\nannotated PVT subtypes based on the marker genes presented in Figure 3 of Gao et al. (2023).\nFor the Shima et al. dataset, we used their original cell type annotations.\nf) Matching scores of each cell cluster to mouse PVT neurons. g) Proportions of mouse PVT\nneuron subtypes from Gao et al. and Shima et al. corresponding to human ExN_CALB clusters.\nAll PVT subtypes refer to mouse PVT neurons; human PVT neurons were not subclassified in\nthis study. aPVT, anterior PVT. The annotated PVT subtypes in pPVT, posterior PVT. AL, anterior\nlateral. AM, anterior medial. PL, posterior lateral. PM, posterior medial. PVT# in Gao et al.\ndenotes the original subtype classification.\n• Consideration of alternative explanations for differences in cell type composition\n• Dissection likely influences which excitatory populations were captured. How are cell\npopulations represented across the 40 individuals?\nWe fully agree with the reviewer that dissection variability may influence the observed\ncompositional differences. To address this concern, we systematically examined the\ndistribution of major cell types and excitatory neuron subpopulations across all 41 individual\nsamples, as presented in the newly added Supplementary Fig. 3b-c. The Tha_ExN_CALB\ncluster was consistently detected across samples, with two outliers: BD13 and CT10.\nImportantly, even when these two samples were excluded, the compositional reduction of\nTha_ExN_CALB remained, supporting the robustness of our findings.\nTo further address the potential influence of dissection variability or anatomical outliers,\nwe\nperformed an additional analysis using repeated random subsampling. Specifically, we\nrandomly selected 10 BD and 10 control samples from the full dataset (21 BD vs. 20 controls)\nand assessed disease-associated fold-changes (FC) using sccomp, while adjusting for Age and\nSex. This approach allowed us to test whether the originally observed Dx-based FCs might have\nbeen driven by a small number of anatomical outliers.\nWe repeated this procedure 1,000 times, and the resulting FC distributions are presented in\nthe newly added Supplementary Fig. 3e. As shown, over 97.5% of the random subsamples\nyielded Dx-based FCs < 0.79 for both Tha_ExN_CALB and thalamic excitatory neurons in total,\nsuggesting that the observed compositional reductions are unlikely to be artifacts driven by a\nsmall number of outlier samples with atypical anatomical positions.\nIn addition, we have added scCODA analysis as suggested by Reviewer 1.\nLine 141 in the revised manuscript\nCellular compositional changes in the medial thalamus of BD\nWe subsequently examined the cellular compositional and transcriptional alterations in the\nmedial thalamus and frontal cortex between 21 BD patients and 20 controls. The thalamus\nand frontal cortex were analyzed separately to enhance precision. Cellular compositions refer\nto the relative proportions of cell clusters within the total population of cell nuclei.\nWe obtained 146,643 cell nuclei in 18 unique cell clusters and 2,139 droplets in two\nNeuCyto clusters after QC (Fig. 3a and Supplementary Fig. 3a). Two NeuCyto clusters were\nremoved as non-nuclei from the subsequent analysis. Tha_ExN_CALB in Fig. 3a corresponded\nto PVT neurons. All major cell classes were consistently detected across the 41 samples\n(Supplementary Fig. 3b). The Tha_ExN_CALB cluster was also robustly identified in all samples\namong thalamic excitatory neurons with two outliers, BD13 and CT10 (Supplementary Fig.\n3c).\nSignificant compositional reductions were observed exclusively in thalamic excitatory\nneurons (FDR < 0.05, FC = 0.54 – 0.59), including PVT neurons (Tha_ExN_CALB, FDR = 0.00488,\nFC = 0.55), as determined by sccomp42 (Fig. 3b, Supplementary Table 3). In contrast, five\noligodendrocytes and astrocytes clusters showed compositional increases (FDR < 0.05, FC =\n1.22 - 1.55). The reduction of Tha_ExN_CALB remained significant even after excluding the\ntwo outlier samples, BD13 and CT10 (FDR = 0.011, FC = 0.592). In examining six major cell\nclasses, we observed a significant compositional reduction in excitatory neurons (FDR =\n0.00525, FC = 0.58) and compositional increases in oligodendrocytes and astrocytes in BD (FDR\n= 0.0168 and 0.004, FC = 1.29 and 1.44, respectively) (Supplementary Fig. 3d). Complementary\nanalysis using scCODA43 also supported a consistent compositional reduction in thalamic\nexcitatory neurons, including PVT neurons (Supplementary Data 4). In the scCODA analysis,\nTha_ExN_CALB and Tha_ExN_PVALB_1 exhibited comparable compositional decreases (FC =\n0.686 and 0.663).\nVariability in tissue dissection may also influence cell-type composition; however,\nsuch variability is difficult to control during sampling and to explicitly model statistically.\nTo\naddress this potential confound, we randomly selected 10 BD and 10 control samples from\nthe full dataset (21 BD vs. 20 controls) and used sccomp to evaluate BD-associated FCs across\n1,000 random subsamples. Notably, over 97.5% of these subsamples yielded diagnosis-\nassociated FCs < 0.79 for both Tha_ExN_CALB and the broader thalamic excitatory neuron\npopulation (Supplementary Fig. 3e), suggesting that the observed compositional reductions\nare unlikely to be driven by a small number of anatomically atypical outlier samples.\nTo validate this finding, we compared the cell density of thalamic excitatory neurons\nand oligodendrocytes by immunohistochemical staining of VGLUT2 and SOX10, calculating cell\ndensity from 548,423 excitatory neurons and 1,819,912 oligodendrocytes in 36 thalamic\nsamples (Supplementary Data 5). This immunohistochemical assessment revealed significant\nreductions of the thalamic excitatory neurons in BD (P = 0.00781, FC = 0.53) and no significant\n66\ndifferences in the oligodendrocytes (P = 0.252, FC = 0.91) (Fig. 3c-d). If incorporating PMI\nand/or pH to the GLM models, the results remained consistent (Supplementary Data 6).\nNote\nthat immunohistochemistry directly assessed cell density, while compositional analysis using\nsnRNA-seq data assessed relative proportional changes. The immunohistochemical\nassessment indicated that the compositional increase in oligodendrocytes and astrocytes in\nBD should be mainly derived from prominent reductions of excitatory neurons, not reflecting\nabsolute cell density increases. If we set Tha_Oligo_1 as a reference cluster in scCODA analysis,\nassuming that the absolute counts of oligodendrocytes were comparable between BD and\ncontrol from IHC results, we consistently observed compositional decreases in thalamic\nexcitatory neurons (Supplementary Data 4).\nFigure 3. Compositional and transcriptional alterations in the thalamus.\nb) Left: Proportional distribution of thalamic cell clusters, comparing BD and control. Total\nnuclei counts are normalized to 10,000, and the total nuclei count in BD and control is adjusted\nto an equivalent sum. Right: Compositional changes estimated by sccomp with log2(Fold-\nchange) on the x-axis and -log10(FDR) on the y-axis. Circle size reflects the number of nuclei\nper cluster, and color denotes major cell classes.\nSupplementary Data\n04.DataS4.scCODA.txt\nscCODA results for thalamic and cortical clusters, including estimated fold-changes and\nposterior probabilities for each covariate.\n7\n05.DataS05.IHC_count_data.txt\nImmunohistochemical cell count data for thalamic excitatory neurons and oligodendrocytes.\n06.DataS06.IHC_GLM.txt\nGeneralized linear model (GLM) results for immunohistochemistry of thalamic excitatory\nneurons and oligodendrocytes.\n• Spatial transciptomics shows Exc_PVALB3, but no Exc_PVALB1/2. Why?\nWe appreciate this important question. As described above, we corrected the HBCA\nintegration and used the revised dataset for the CytoSPACE analysis. In the original analysis,\nthe CytoSPACE deconvolution was performed at an excessively high resolution, partly due to\nthe erroneous inclusion of hypothalamic samples containing numerous Splatter clusters. In the\nupdated integration, only ExN_CALB, ExN_PVALB_1, and ExN_PVALB_2 were observed as\nthalamic excitatory neurons. The original ExN_PVALB_1 and ExN_PVALB_3 clusters were highly\nsimilar and are now indistinguishable in the revised integration, consistent with the case-\ncontrol thalamus analysis (Fig. 3). In retrospect, distinguishing ExN_PVALB_1 from\nExN_PVALB_3 in the original Fig. 2 was inadequate, given our focus on PVT neurons.\nIn the revised Fig. 2, ExN_PVALB_1 was primarily detected outside the PVT area. The\nExN_PVALB_2 subtype was not detected, consistent with our original submission. This finding\naligns with the absence of ExN_PVALB_2 in the CT06 snRNA-seq dataset, thereby supporting\nthe validity of the CytoSPACE-based spatial deconvolution. In the HBCA integration dataset,\nthe CT06 dataset comprised 4,510 nuclei in total, including 127 ExN_CALB, 364 ExN_PVALB_1,\nand 0 ExN_PVALB_2, as shown in Supplementary Data 11. The absence of ExN_PVALB_2 in\nCT06 is likely attributable to the anatomical position of the cryosection.\nTo corroborate this finding, we performed additional CytoSPACE deconvolution using two\nsnRNA-seq datasets: an integrated set comprising our 20 controls (excluding HBCA samples)\nand the single CT06 snRNA-seq dataset. Both analyses yielded cellular distributions similar to\nthose observed in the original deconvolution (Supplementary Fig. 13).\nLine 851 in the revised manuscript\nDeconvolution of spatial transcriptomic data with snRNA-seq cell clusters\nWe deconvoluted the Visium spatial transcriptomic data using cell-type expression profiles\nobtained from snRNA-seq via CytoSPACE39 (v1.1.0) with default parameters to estimate the\nspatial distribution of each cell type. The reference cell types included excitatory/inhibitory\nneurons, astrocytes, OPCs, oligodendrocytes, and microglia, as defined in the integrated set\ncomprising the HBCA and our dataset (see Fig. 2a). The low representation of the\nExN_PVALB_2 cluster in Fig. 2d was consistent with its low/no proportion in the CT06 snRNA-\nseq dataset (Supplementary Data 11), supporting the validity of our deconvolution approach.\nAs a robustness check, we repeated the deconvolution using either our control thalamic\ndataset or CT06-derived data alone, both of which yielded consistent results (Supplementary\nFig. 13e-f).\nFigure 2. Characterization of human PVT neurons.\nd) Spatial localization of cell clusters as predicted by CytoSPACE deconvolution, illustrating\ntheir estimated anatomical distribution within the thalamus. Color intensity in the sidebars\nindicates the estimated number of cells per spot, providing a quantitative measure of cell\ndensity across regions. e) Merged view of the spatial clusters from panel d for comparison with\npanel c.\ne) CytoSPACE deconvolution based solely on our snRNA-seq dataset (excluding the HBCA\ndataset). h) CytoSPACE deconvolution using only the CT06 snRNA-seq dataset.\nSupplementary Data\n11.DataS11.HBCA.profiles.txt\nSample profiles of HBCA data used in this study, including cell cluster counts.\n• Is neuronal loss in BD limited to PVT (as defined by Exc_CALB cluster) or only detectable there\ndue to more consistent dissection of the midline?\nWhile the prominent reduction was observed in the ExN_CALB population (PVT neurons in our\nstudy), other thalamic excitatory neuron subtypes (Tha_ExN_PVALB_1 and\nTha_ExN_PVALB_2) also exhibited comparable compositional reductions, as illustrated in Fig.\n3b. In response to Reviewer 1’s suggestion, we additionally applied scCODA, which supported\nthese findings and revealed a slightly more pronounced decrease in Tha_ExN_PVALB_1.\nImmunohistochemistry (IHC) further confirmed a reduction in excitatory neurons across the\nmedial thalamus. This analysis was not restricted to the PVT neurons but encompassed all\nthalamic excitatory neurons. Given that PVT neurons represent only a small fraction of total\nexcitatory neurons (16.2%), the observed IHC reduction cannot be attributed solely to changes\nwithin the PVT. These findings suggest that the loss of PVALB-positive neurons also contributes\nto the broader excitatory neuron decline, indicating a more generalized neuronal vulnerability\nwithin the medial thalamus.\nAs the central focus of our study is PVT, tissue sampling was anatomically centered on this\nregion. While neighboring nuclei, such as the mediodorsal (MD) nucleus, are of significant\ninterest, they were not explicitly targeted in this analysis. Nonetheless, the observed\ncompositional changes may implicate reductions in excitatory neuron populations extending\ninto adjacent nuclei, potentially including the MD. A more definitive characterization of these\nchanges will require future studies with targeted sampling of these adjacent regions.\nWe have revised the compositional analysis section to present the results more clearly.\nLine 141 in the revised manuscript\nCellular compositional changes in the medial thalamus of BD\nWe subsequently examined the cellular compositional and transcriptional alterations in the\nmedial thalamus and frontal cortex between 21 BD patients and 20 controls. The thalamus\nand frontal cortex were analyzed separately to enhance precision. Cellular compositions refer\nto the relative proportions of cell clusters within the total population of cell nuclei.\nWe obtained 146,643 cell nuclei in 18 unique cell clusters and 2,139 droplets in two\nNeuCyto clusters after QC (Fig. 3a and Supplementary Fig. 3a). Two NeuCyto clusters were\nremoved as non-nuclei from the subsequent analysis. Tha_ExN_CALB in Fig. 3a corresponded\nto PVT neurons. All major cell classes were consistently detected across the 41 samples\n(Supplementary Fig. 3b). The Tha_ExN_CALB cluster was also robustly identified in all samples\namong thalamic excitatory neurons with two outliers, BD13 and CT10 (Supplementary Fig.\n3c).\nSignificant compositional reductions were observed exclusively in thalamic excitatory\nneurons (FDR < 0.05, FC = 0.54 – 0.59), including PVT neurons (Tha_ExN_CALB, FDR = 0.00488,\nFC = 0.55), as determined by sccomp42 (Fig. 3b, Supplementary Table 3). In contrast, five\noligodendrocytes and astrocytes clusters showed compositional increases (FDR < 0.05, FC =\n1.22 - 1.55). The reduction of Tha_ExN_CALB remained significant even after excluding the\ntwo outlier samples, BD13 and CT10 (FDR = 0.011, FC = 0.592). In examining six major cell\nclasses, we observed a significant compositional reduction in excitatory neurons (FDR =\n0.00525, FC = 0.58) and compositional increases in oligodendrocytes and astrocytes in BD (FDR\n= 0.0168 and 0.004, FC = 1.29 and 1.44, respectively) (Supplementary Fig. 3d). Complementary\nanalysis using scCODA43 also supported a consistent compositional reduction in thalamic\nexcitatory neurons, including PVT neurons (Supplementary Data 4). In the scCODA analysis,\nTha_ExN_CALB and Tha_ExN_PVALB_1 exhibited comparable compositional decreases (FC =\n0.686 and 0.663).\nVariability in tissue dissection may also influence cell-type composition; however,\nsuch variability is difficult to control during sampling and to explicitly model statistically.\nTo\naddress this potential confound, we randomly selected 10 BD and 10 control samples from\nthe full dataset (21 BD vs. 20 controls) and used sccomp to evaluate BD-associated FCs across\n1,000 random subsamples. Notably, over 97.5% of these subsamples yielded diagnosis-\nassociated FCs < 0.79 for both Tha_ExN_CALB and the broader thalamic excitatory neuron\npopulation (Supplementary Fig. 3e), suggesting that the observed compositional reductions\nare unlikely to be driven by a small number of anatomically atypical outlier samples.\nTo validate this finding, we compared the cell density of thalamic excitatory neurons\nand oligodendrocytes by immunohistochemical staining of VGLUT2 and SOX10, calculating cell\ndensity from 548,423 excitatory neurons and 1,819,912 oligodendrocytes in 36 thalamic\nsamples (Supplementary Data 5). This immunohistochemical assessment revealed significant\nreductions of the thalamic excitatory neurons in BD (P = 0.00781, FC = 0.53) and no significant\n70\ndifferences in the oligodendrocytes (P = 0.252, FC = 0.91) (Fig. 3c-d). If incorporating PMI\nand/or pH to the GLM models, the results remained consistent (Supplementary Data 6).\nNote\nthat immunohistochemistry directly assessed cell density, while compositional analysis using\nsnRNA-seq data assessed relative proportional changes. The immunohistochemical\nassessment indicated that the compositional increase in oligodendrocytes and astrocytes in\nBD should be mainly derived from prominent reductions of excitatory neurons, not reflecting\nabsolute cell density increases. If we set Tha_Oligo_1 as a reference cluster in scCODA analysis,\nassuming that the absolute counts of oligodendrocytes were comparable between BD and\ncontrol from IHC results, we consistently observed compositional decreases in thalamic\nexcitatory neurons (Supplementary Data 4).\nFigure 3. Compositional and transcriptional alterations in the thalamus.\nb) Left: Proportional distribution of thalamic cell clusters, comparing BD and control. Total\nnuclei counts are normalized to 10,000, and the total nuclei count in BD and control is adjusted\nto an equivalent sum. Right: Compositional changes estimated by sccomp with log2(Fold-\nchange) on the x-axis and -log10(FDR) on the y-axis. Circle size reflects the number of nuclei\nper cluster, and color denotes major cell classes.\nSupplementary Data\n04.DataS4.scCODA.txt\nscCODA results for thalamic and cortical clusters, including estimated fold-changes and\nposterior probabilities for each covariate.\n71\n05.DataS05.IHC_count_data.txt\nImmunohistochemical cell count data for thalamic excitatory neurons and oligodendrocytes.\n06.DataS06.IHC_GLM.txt\nGeneralized linear model (GLM) results for immunohistochemistry of thalamic excitatory\nneurons and oligodendrocytes.\n• Overinterpretation of mechanistic link between excitatory neurons and microglia\n• SHISA9 and SYNDIG1 are both AMPAR auxiliary subunits, but an interaction between them is\nnot known; SYNDIG1 encodes a transmembrane protein that is known to interact with AMPAR in\nneurons. Its function in microglia is not known. More evidence is needed of any claimed\ninteractions in microglia.\n• Discussion would benefit from a more nuanced handling of the data, given its limitations,\nand more attention to areas where evidence is less clear and multiple conclusions are possible,\ne.g.,\nWe fully agree with the reviewer’s comment and acknowledge that our previous interpretation\noverstated the relationship between SYNDIG1 and SHISA9 proteins. As correctly noted, both\nproteins are AMPA receptor (AMPAR) auxiliary subunits, but there is no evidence supporting a\ndirect protein-protein interaction between them. To our knowledge, the mechanistic link\nbetween excitatory neurons and microglia regarding SYNDIG1 and SHISA9 proteins remains\nunknown. We have thoroughly revised the description of this point in the manuscript and the\nrelevant section of the Discussion accordingly.\nIn response to Reviewer #1’s comment, we have incorporated RNAscope validation of SYNDIG1\nexpression in microglia. The preferential localization of SYNDIG1-positive microglia in cortical\ngray matter aligns with our meta-gene programs, which implicate interactions between\nmicroglia and excitatory neurons. Nonetheless, we recognize that this data-driven construct\nwarrants further experimental validation.\nAbstract in the original manuscript\nREMOVED: This interaction was further exemplified by the concurrent downregulation of\nSYNDIG1 in microglia and SHISA9 in PVT neurons, both of which are involved in AMPA receptor\nsignaling pathways.\nLine 279 in the original manuscript\nREMOVED: The coordinated downregulation of SYNDIG1 in microglia and SHISA9 in excitatory\nneurons exemplifies a probable disturbance in microglia-neuron interactions for synapses in\nthe BD thalamus, indicating a plausible pathological mechanism of SHISA9, one of the 36\nputative BD-causal genes1.\nLine 415 in the original manuscript\nREMOVED: The downregulation of SYNDIG1 in microglia and SHISA9 in PVT neurons might be\nan example of dysregulated molecular interactions within the Excit-Micro meta-gene programs\nassociated with BD.\nLine 314 in the revised manuscript\nAt the gene level, SYNDIG1 and SHISA9 encode auxiliary subunits for AMPA receptor (AMPAR)\ncomplex56,57, with both genes’ robustness scores exceeding 0.7 (Supplementary Table 4).\nNotably, SYNDIG1 was prominently expressed in microglia, while SHISA9 was primarily\nexpressed in excitatory neurons (Fig. 5d). SHISA9 consisted of core PPI networks of BD-\nassociated downregulated genes in PVT neurons (Fig. 4d) and overlapped in the\ndownregulated genes in Tha_ExN_PVALB_1, as one of the 36 putative BD-causal genes1.\nBeyond microglia and excitatory neurons, thalamic T cells were included in this meta-gene\nprogram. Fig. 5e illustrates potential ligand-receptor interactions between neurons and\nmicroglia/T cells, as analyzed through LIANA58. Although we did not detect concurrently\ndownregulated DEGs as ligand–receptor pairs, KCND2 in Tha_ExN_CALB, LAMB1 in\nTha_ExN_PVALB_2, and KCNQ1 in Tha_Micro_1 were identified as downregulated DEGs\namong the LIANA-based ligand-receptor pairs. These genes encode ion channels or\nextracellular matrix proteins that may potentially contribute to Excit–Micro interactions.\nHowever, as this analysis is an inference, the precise ligand–receptor mechanisms remain to\nbe elucidated in future studies.\nLine 556 in the revised manuscript\nBeyond cell-type-based analysis, we identified trans-cellular meta-gene programs between\nexcitatory neurons and microglia (Excit-Micro meta-gene programs), likely reflecting\nphysiological synaptic interactions between these cell types. Microglia are known to modulate\nsynaptic plasticity and neuronal circuitry through dynamic interactions with neurons76.\nNotably, such meta-gene programs in the thalamus are disrupted in BD, further characterized\nby a significant decrease in synaptic gene expression in both PVT neurons and microglia. Ling\net al. reported downregulation of the synaptic neuron and astrocyte program in\nschizophrenia10; our findings highlight downregulation of the synaptic genes in the thalamic\nneuron and microglia in BD. The concurrent downregulation of SYNDIG1 in microglia and\nSHISA9 in excitatory neurons suggests disruptions of AMPAR signaling, as both genes encode\nauxiliary subunits of the AMPAR complex56,57. TIAM1, another downregulated gene in\nmicroglia, has also been implicated in regulating synaptic AMPAR function77. The preferential\nlocalization of SYNDIG1-positive microglia in cortical gray matter is consistent with this\ninterpretation. Nonetheless, these meta-gene programs remain data-driven constructs, and\nthe roles of SYNDIG1 in microglia—as well as their mechanistic links to SHISA9 in PVT\nneurons—are currently unknown. Future studies employing animal and in vitro cellular\nmodels will be crucial for elucidating these mechanisms.\nAdditionally, LIANA does not support a SHISA9-SYNDIG1 interaction. In the manuscript, LIANA\nwas employed to investigate potential microglia-excitatory neuron interactions suggested by\nNMF-based module analysis. However, no ligand-receptor pair was identified in which both\nthe ligand and receptor genes were downregulated.\nNevertheless, several DEGs were detected that may participate in intercellular signaling as\neither the ligand or receptor: KCND2 in Tha_ExN_CALB, LAMB1 in Tha_ExN_PVALB_2, and\nKCNQ1 in Tha_Micro_1. These genes encode ion channels or extracellular matrix proteins that\nmay indirectly contribute to cell-cell interactions. However, the precise ligand-receptor\nmechanisms remain unresolved, and we now explicitly state that cell-cell communication\nanalysis remains a subject for future investigation.\nLine 320 in the revised manuscript\nFig. 5e illustrates potential ligand-receptor interactions between neurons and microglia/T\ncells, as analyzed through LIANA58. Although we did not detect concurrently downregulated\nDEGs as ligand–receptor pairs, KCND2 in Tha_ExN_CALB, LAMB1 in Tha_ExN_PVALB_2, and\nKCNQ1 in Tha_Micro_1 were identified as downregulated DEGs among the LIANA-based\nligand-receptor pairs. These genes encode ion channels or extracellular matrix proteins that\nmay potentially contribute to Excit–Micro interactions. However, as this analysis is an\ninference, the precise ligand–receptor mechanisms remain to be elucidated in future studies.\nd) UMAP projection illustrating the expression patterns of SYNDIG1 and SHISA9.\nFigure 8. Comparative analysis between the thalamus and cortex.\nc) Relative expression of CACNA1C and SHISA9 across all clusters in the thalamus-cortex\nintegration dataset.\n• Findings in PVT may not apply to the rest of the thalamus\nAs pointed out, the findings from the PVT may not be directly applicable to the entire thalamus.\nOther areas outside of the medial thalamus need further investigations. We have added this\nargument in the discussion.\nLine 483 in the revised manuscript\nThe observed reduction in thalamic excitatory neurons, including PVT neurons, may be\nassociated with the reduced thalamus volume reported in BD19, although our analysis was\nlimited to the medial thalamus. Given that variability in tissue dissection represents a potential\nconfounder, future studies encompassing the entire thalamus will be required to determine\nwhether this reduction is specific to the medial thalamus or extends to the whole structure. It\nremains unclear whether the present findings alone can account for the volumetric reductions\nobserved in BD.\n• Avoid anthropomorphic statements such as “Mice with recurrent depressive episodes”\n• Claims about high-frequency firing and role of KCNQ3 in potential response to high firing are\ntoo speculative; citations would be needed to support such a claim for this gene\nWe have revised the Discussion section to adopt a more cautious tone. Specifically, we\nremoved anthropomorphic phrases such as “mice with recurrent depressive episodes” and\nadded a detailed explanation for KCNQ3 and potential hyper-excitability. We appreciate the\nreviewer’s emphasis on scientific rigor and nuance.\nKCNQ3 encodes the Kv7.3 subunit, which forms heteromers with Kv7.2 to generate the non-\ninactivating, slowly activating M-current near resting membrane potential (Jentsch. Nature\nReviews Neuroscience 2000, 1:21-30; Wang et al. Science 1998, 282:1890-1893). These\nchannels modulate neuronal excitability by affecting input resistance, the action potential\nthreshold, and the medium afterhyperpolarization (mAHP), thereby controlling firing\nfrequency (Yue and Yaari. Journal of Neurophysiology 2006, 95:3480-95). Suppression of the\nM-current lowers the action potential threshold, increases afterdepolarization, and\ndepolarizes the axonal resting potential, leading to enhanced excitability and a predisposition\nto high-frequency firing (Brown and Adams. Nature 1998, 283:673-676; Delmas and Brown.\nNature Reviews Neuroscience 2005, 6: 850-862; Selyanko and Brown. Neuron 1996, 16:151-\n162).\nKv7.2/7.3 channels are enriched at the axon initial segment and nodes of Ranvier, where they\nstabilize the resting potential. This stabilization, which is distinct from the inhibitory effect of\nM current at the soma and dendrites, inhibits Na+ channel inactivation by preventing excessive\ndepolarization, thereby maintaining their availability for reliable action potential initiation and\npropagation and preserving normal, rather than excessive, excitability (Shah et al. PNAS 2008,\n105:7869-7874; Battefeld et al. Journal of Neuroscience 2004, 34:3719-3732). Thus,\ndownregulation of Kv7.3 could be associated with excessive excitability.\nAs the above explanation is too lengthy for this manuscript, we have added the shorter version\nin the main manuscript.\nLine 502 in the revised manuscript\nOur prior findings showed that regulation of neurotransmission in the PVT via tetanus toxin or\nDesigner Receptors Exclusively Activated by Designer Drugs (DREADDs) induced recurrent\ndepression-like behavioral changes even in wildtype mice24,25.\nLine 526 in the revised manuscript\nPotassium channels, such as KCNQ3—a gene identified as a key DEG in BD—are pivotal\nregulators of high-frequency firing. KCNQ3 encodes the Kv7.3 subunit, which assembles with\nKv7.2 to form heteromeric channels mediating the non-inactivating, slowly activating M-\ncurrent near the resting membrane potential73. The M-current regulates neuronal excitability\nby modulating input resistance, action potential threshold, and the medium\nafterhyperpolarization, thereby influencing firing frequency74. Inhibition of the M-current\nlowers the action potential threshold, enhances afterdepolarization, and depolarizes the\naxonal resting potential, resulting in increased excitability and a propensity for high-frequency\nfiring75. Accordingly, reduced Kv7.3 expression may contribute to neuronal hyperexcitability.\n• Lack of relevance of thalamic organoid integration analysis\n• Authors did not sample the entire thalamus to claim that organoids recapitulate properties\n74\nof PVT\n• This analysis is not really related to the rest of the paper\n• Insufficient methodological detail on Quality Control\n• Authors should describe QC and filtering process in more detail\nWe agree with the reviewer that the organoid analysis, while interesting, was peripheral to the\ncentral findings. Another reviewer also pointed out that this section is irrelevant and should\nbe removed. Initially, we had incorporated this part for the final discussion. In retrospect, this\npart is not essential for our manuscript. We have removed the organoid comparison from the\nmain text and supplementary figures.\nLine 137 in the original manuscript\nREMOVED: To explore opportunities for simulating human PVT neurons in vitro, we\ninvestigated whether human thalamic organoids contain PVT-like neuronal populations. PVT-\nlike neurons were identified in the thalamic organoids created by Kiral et al.41 (Supplementary\nFig. 2d-e). These neurons expressed PVT marker genes, including SLC17A6, CALB2, and SNCA.\nTheir expression profile suggests the feasibility of modeling human PVT neurons in vitro.\nLine 463 in the original manuscript\nREMOVED: Despite these limitations, this study indicates the disturbances of PVT neurons as\na core pathology in BD, providing a promising foundation for developing new diagnostic and\ntherapeutic strategies. Neuroimaging techniques targeting the PVT, such as MRI and PET,\ncould enhance clinical assessment for BD. PVT-like neurons in thalamic organoids underscore\nthe feasibility of modeling PVT neurons from patient-derived peripheral samples, which can\naccelerate the discovery of novel pharmacological interventions to the molecular mechanisms\nunderlying BD. Integrating multi-modal methodologies, encompassing cellular, molecular, and\nneuroimaging techniques, will provide a more comprehensive understanding of the\npathological contributions of PVT in BD. PVT neurons represent a promising target for the\ndevelopment of advanced diagnostic and therapeutic strategies aimed at improving patient\noutcomes.\nLine 833 in the original manuscript\nREMOVED: Kiral et al. data can be obtained from GEO with accession number GSE210720.\nSupplementary Figure\nREMOVED: Supplementary Figure S2d-e in the original manuscript\n• No vascular or ependymal cells seen in the final annotated cell types: were they filtered due\nto low quality or not included because too few cells were represented?\nVascular cells were annotated in the integrated datasets shown in Fig. 2a, Fig. 3a, and Fig. 6a,\nbut with low proportions. Specifically, Endo (endothelial cells) and Peri (pericytes) in Fig. 2a\ncomprised 0.70% of total nuclei; the Tha_Vascular cluster in Fig. 3a accounted for 0.216%; and\nthe FrC_Vascular cluster in Fig. 6a represented 0.196% of total nuclei in each respective\nintegration. Although further subclassification into endothelial cells and pericytes would be\n75\npossible at higher resolution, we did not pursue this, as our primary focus was the\ncharacterization of PVT neurons.\nThe integration of our thalamic samples with HBCA thalamic samples (Fig. 2a) contained a\nhigher proportion of vascular cells than the integration of our thalamic samples alone (Fig. 3a).\nThis likely reflects a technical limitation, as low-abundance cell types such as vascular cells are\nmore susceptible to loss during clustering and integration procedures.\nVascular cells were not retained in the final annotated dataset (Fig. 1b), likely due to the\nintegration of a large number of clusters from both thalamus and cortex, which comprise\ndistinct cell types. Such integration may have lacked the resolution needed to preserve low-\nabundance populations, including vascular cells.\nEpendymal cells were not annotated in our dataset, likely because they are difficult to isolate\nusing our nuclear extraction protocol. These cells are embedded within a dense extracellular\nmatrix, and our protocol employed a relatively mild detergent concentration (0.05% Tween\n20) to preserve nuclear integrity and RNA quality in long-stored postmortem brain samples,\nincluding those with a postmortem interval (PMI) exceeding 50 hours. This relatively mild\ndetergent condition may also have reduced the yield of vascular cells, as vascular tissue is\nenriched in collagen fibers, which are more resistant to dissociation into single cells. Resistance\nto detergent-based cell nuclei isolation procedures would be another reason for the low\nrepresentation of vascular cells in our study.\nFigure 2a (Integration of HBCA and our thalamic dataset)\nFigure 3a (Integration of our thalamic dataset)\n• Insufficient methodological detail on anchor-based integration\n• Was CCA or RPCA approach used, and why?\nFor anchor-based processing, we employed Seurat’s “anchor-based transfer” approach using\nthe FindTransferAnchors and TransferData functions, as described in the official tutorial\n“Mapping and annotating query datasets”\nhttps://satijalab.org/seurat/articles/integration_mapping\nWe used the web-based Azimuth platform (https://azimuth.hubmapconsortium.org/), which\ninternally uses the same functions. We have revised the manuscript to include this information.\nWhile CCA- or RPCA-based integration could be used to align mouse and human cell types,\nwe\nopted for anchor-based transfer using FindTransferAnchors and TransferData functions for the\nfollowing reasons:\n1. It allows for explicit, cell-to-cell mapping scores, providing a quantitative measure of how\n77\nwell a human cell type corresponds to a specific mouse cell type.\n2. Mouse and human expression profiles differ substantially, and anchor-based transfer is\nbetter suited for cross-species annotation without enforcing artificial alignment.\n3. The mouse data from Shima et al. were generated using Quartz-seq2, a platform with\ndifferent technical characteristics from 10x Genomics Chromium system, further\nsupporting the use of transfer-based methods over direct integration.\nAccordingly, we have revised the Methods section and uploaded a script for processing Gao et\nal. data in our GitHub repository.\nLine 782 in the revised manuscript\nNeuronal cell clusters in the frontal cortex and glial cell clusters were annotated using cell-type\nlabel transfer with Azimuth, using a cortex reference29. Azimuth internally employed anchor-\nbased transfer using FindTransferAnchors and TransferData functions in Seurat.\nLine 869 in the revised manuscript\nMatching mouse PVT neurons to the human PVT neurons\nTo identify cell types in our human snRNA-seq data that correspond most closely to mouse\nPVT neurons, we performed anchor-based label transfer using the Seurat functions\nFindTransferAnchors and TransferData. As reference datasets of mouse PVT neurons, we used\nthe snRNA-seq data from Gao et al. (2023)40 and the single-cell RNA-seq data from Shima et\nal. (2023)23. The mouse gene name was converted to the corresponding human gene name in\nthe Ensemble database. Prior to anchor transfer, each dataset was independently processed\nusing Seurat, including normalization, scaling, and principal component analysis (PCA). For the\nGao et al. dataset, we followed the authors’ published analytical pipeline to identify PVT\nsubtypes, with the exception of PCA dimensions to improve clustering resolution, and\nannotated PVT subtypes based on the marker genes presented in Figure 3 of Gao et al. (2023).\nFor the Shima et al. dataset, we used their original cell type annotations.\n• Insufficient exploration of confounding variables for DGE (and cellular composition)\n• Why were the covariates selected as described? The authors did not appear to have\ninvestigated other quality-related variables and explore variance they explained.\nWe appreciate the reviewer’s insightful comment and apologize for omitting the detailed\nrationale for covariate selection. In response, we have now added a comprehensive\ndescription of the covariate selection process and thoroughly revised the relevant section in\nthe Methods.\nCovariate Selection Strategy\nInitially, we considered a wide range of potential covariates: diagnosis (Dx), Sex, Age,\nPostmortem interval (PMI), brain pH (pH), race, cause of death, refrigerator delay, substances\nat death, experimental/sequencing batch. Given our modest sample size (N = 41), covariates\nshould be selected to minimize the risk of overfitting and multicollinearity.\n- Race: The cohort was heavily biased toward Caucasian individuals, with only one Asian\nsample. Thus, race was excluded from statistical modeling due to lack of variability.\n- Cause of death: A significant imbalance was observed, with suicide overrepresented\namong cases (Fisher’s exact test, 2.33 × 10-8). To reduce multicollinearity, this covariate was\nexcluded.\n- Refrigerator delay (RD): RD contained three missing values, which is substantial for a\nsample of 41. Furthermore, RD was significantly correlated with PMI (Pearson’s correlation,\nP = 0.0181). We prioritized PMI as a relatively more relevant potential covariate.\n- Medication: Due to the wide heterogeneity of substances at death and missing information\nin several donors, subgrouping was not feasible with the available sample size.\n- Experimental and sequencing batches: We had 21 batches for nuclei isolation and 16\nsequencing batches in thalamic analysis. Modeling batch effects failed due to the large\nnumber of batches relative to sample size (N = 41).\nThe remaining candidate covariates were Dx, sex, age, PMI, and pH. Among them, pH was\nparticularly correlated with Dx (GLM P = 0.0074; see revised Supplementary Fig. 14b),\nsuggesting potential collinearity. Nonetheless, we tested models with all five covariates (Dx,\nSex, Age, PMI, and pH), as these are commonly used in postmortem brain studies.\nCompositional Analysis\nFor the compositional analysis, model fit was evaluated using Akaike Information Criterion\n(AIC) from generalized linear models (GLM) employing negative binomial distribution (glm.nb,\nMASS package), given that the sccomp package does not provide scores indicating model fit\nlike AIC. We compared multiple candidate models and ranked them by summed AIC rankings\n(rankSum) across both thalamic and cortical datasets.\nThe model including Dx, age, and sex (Dx + Age + Sex) consistently exhibited the lowest AIC by\nrankSum, indicating the best overall fit. The contributions of PMI and pH differed between the\nthalamic and cortical analyses. We selected PMI as the secondary covariate because pH\nshowed a stronger correlation with diagnosis and represents a less direct measure of tissue\nintegrity for both the thalamus and cortex (as described in the original manuscript, Line 672).\nBased on these findings, we selected Dx + Age + Sex as the primary covariate combination for\ncompositional modeling, and Dx + Age + Sex + PMI as the secondary covariate combination.\nIn line with the reviewer’s suggestion, we also assessed the relative contribution of each\ncovariate using Type III likelihood ratio chi-square tests (Anova(type = 3) from the car package).\nThe results are provided in Supplementary Data 16. Dx, Sex, and Age consistently\ndemonstrated significant contributions to compositional variation.\nDifferential Expression Analysis\nIn the differential gene expression (DEG) analysis, we adopted the same covariate combination\nused in the compositional analysis. In response to the reviewer’s comments, we calculated the\nAIC of generalized linear models (GLMs) for each gene’s log-transformed expression [log(p +\n1)] to specifically evaluate model fit when incorporating PMI and/or pH, as MAST does not\nprovide scores indicating model fit like AIC. Gene expression was aggregated at the specimen\nlevel and limited to genes expressed in at least 10% of cells within each cluster, consistent with\nthe DEG criteria. This GLM-based modeling was performed solely to assess the impact of\ncovariates; the DEG analysis itself was conducted using zlm in the MAST package.\nAmong the tested models, Dx + Age + Sex yielded the lowest total AIC across genes by clusters,\nfollowed by Dx + Age + Sex + PMI. This supports the use of Dx + Age + Sex as the primary\ncovariate combination and Dx + Age + Sex + PMI as a secondary model. In addition, we\nassessed the relative contribution of each covariate using Type III chi-square statistics\n(Anova(type = 3) from the car package) for both models. As these data are too huge to submit\nfor the journal web system, we deposited these data in figshare.\nWe have revised the Methods section to incorporate these details. The updated description\nenhances the clarity and transparency of our statistical approach.\nLine 897 in the revised manuscript\nCovariate selection in compositional analysis\nInitially, we considered a wide range of potential covariates: diagnosis (Dx), sex, age,\npostmortem interval (PMI), brain pH (pH), race, cause of death, refrigerator delay, substances\nat death, experimental batch, and sequencing batch. Given our modest sample size (N = 41),\ncovariates should be selected to minimize the risk of overfitting and multicollinearity. Race\nwas excluded from statistical modeling because the cohort was heavily biased toward\nCaucasian individuals, with only one Asian sample. Cause of death was excluded to reduce\nmulticollinearity, because it exhibited a significant imbalance with suicide overrepresented\namong BD cases (two-sided Fisher’s exact test, P = 2.33 × 10-8). Refrigerator delay contained\nthree missing values, and RD was correlated with PMI (Pearson’s correlation, P = 0.0181). We\nprioritized PMI as a relatively more relevant potential covariate. Due to the wide\nheterogeneity of substances at death and missing information in several donors, subgrouping\nbased on substances at death was not feasible in our dataset. Regarding experimental and\nsequencing batches, we had 21 batches for nuclei isolation and 16 sequencing batches in\nthalamic analysis. Modeling batch effects was unsuccessful due to the large number of batches\nrelative to the sample size (N = 41). Although experimental and sequencing batches could not\nbe incorporated into the subsequent statistical modeling due to their high dimensionality\nrelative to sample size, we confirmed that the identified cell clusters were broadly and\nconsistently represented across batches (Supplementary Fig. 14a, Supplementary Data 13\nand 14).\nThe remaining candidate covariates were Dx, sex, age, PMI, and pH. Among them, pH\nvalues obtained from cerebellum samples, may not directly represent the specimen\nconditions of the thalamus and frontal cortex; given probable regional variability, we treated\npH as an indirect indicator for this study. Notably, pH and PMI were associated with Dx, with\npH showing a stronger signal (Supplementary Fig. 14b-c). These associations likely reflect\ndifferences in cause of death, as BD samples predominantly originated from individuals who\ndied by suicide, unlike control samples. Nonetheless, we tested statistical models with all five\ncovariates (Dx, Sex, Age, PMI, and pH), as commonly used covariates in postmortem brain\nstudies.\nFor the compositional analysis, model fit was evaluated using Akaike Information\nCriterion (AIC) from generalized linear models (GLM) employing negative binomial distribution\n(glm.nb in MASS [v7.3-60] package), because the sccomp and scCODA packages do not provide\nscores indicating model fit like AIC. In the GLM, total nuclei counts were normalized to 10,000,\nwith each cluster’s count rounded to integers to reduce variability in total nuclei counts across\nindividuals. We compared multiple candidate models and ranked them by summed AIC\nrankings across both thalamic and cortical datasets (see rankSum in Supplementary Data 15).\nThe model including Dx, age, and sex (Dx + Age + Sex) showed the lowest overall AIC, indicating\nthe best model fit. While the contributions of PMI and pH differed between the thalamic and\ncortical analyses, we selected PMI as the secondary covariate because pH showed a stronger\ncorrelation with diagnosis and represents a less direct measure of thalamus and cortex tissues.\nBased on these findings, we selected “Dx + Age + Sex” as the primary covariate combination\nfor compositional modeling, and “Dx + Age + Sex + PMI” as the secondary covariate\ncombination. The relative contribution of each covariate was assessed by Type III likelihood\nratio chi-square tests (Anova(type = 3) from the car package, Supplementary Data 16).\nAs support for the above covariate selection, adding PMI did not increase the model\nfitness in the sccomp analysis, without changes of elpd_diff/se_diff > 5. More than five of\nelpd_diff/se_diff change is suggested for significant model improvement according to the\ndevelopers’ guideline42.\nLine 985 in the revised manuscript\nDifferentially expressed gene analysis\nDifferentially expressed genes (DEGs) were analyzed within each cell cluster in BD compared\nto controls using zlm function in MAST44 (v1.28.0), a statistical package for single-cell RNA-seq\nDEG analysis. The model included Dx, age, sex, and cngeneson (a covariate for cell quality) as\nfixed effects, and sample ID as a random effect. Sample-level covariates were selected based\non the prior compositional analysis and further validated by evaluating AIC values from GLMs\nfitted to log-transformed gene expression values [log(p + 1)]. In this GLM-based model to\nselect covariates, gene expression was aggregated at the specimen level and restricted to\ngenes expressed in at least 10% of cells within each cluster, consistent with the downstream\nDEG filtering criteria. These GLMs were used solely to assess covariate impact. Among the\ntested covariate combinations, “Dx + Age + Sex” yielded the lowest overall AIC across gene-\ncluster pairs, followed by “Dx + Age + Sex + PMI”, as determined by cumulative AIC rank sums\n(External Data, see Data Availability). These results support the use of “Dx + Age + Sex” as the\nprimary model and “Dx + Age + Sex + PMI” as a secondary model. The relative contribution of\neach covariate was further quantified using Type III chi-square statistics (Anova(type = 3) from\nthe car package), as reported in External Data. Following MAST guidelines, we incorporated\ncngeneson, a variable of scaled gene expression counts, to adjust for cell nuclei quality, where\nhigher gene detection rates suggest better nuclei quality within the same cluster.\nSupplementary Data\n15.DataS15.CompositionalAnalysis.GLM.AIC.txt\nAkaike information criterion (AIC) assessment for the GLM in compositional analysis.\n16.DataS16.CompositionalAnalysis.GLM.Anova_ChiSq.txt\nAnalysis of variance (ANOVA) Chi-square assessment for the GLM in compositional analysis.\nExternal Data (available in figshare)\nBDsnRNAseq.thalamus.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each thalamic cluster.\nBDsnRNAseq.cortex.allclusters.logLM_Anova_CR01.model_comparison.txt.gz\nAIC and ANOVA Chi-square assessments for covariate combinations in the GLM of log-\nnormalized summed gene expression in each cortical cluster.\no Defining BD risk genes\no Reliance on MAGMA-based gene level scores from GWAS has significant limitations\no What about genes implicated by rare lof alleles in exome studies?\nWe appreciate the reviewer’s suggestion to strengthen the genetic integration. As correctly\nnoted, interpreting GWAS results based solely on MAGMA is preliminary. The recent PGC4\nstudy (O’Connell Nature 2025) identified 36 “credible genes” through integrative approaches,\nof which only SHISA9 overlapped with the genes identified in our analysis. We acknowledge\nthat only one gene among the 36 credible genes was associated with our results.\nHowever, given the polygenic nature of BD, it is unlikely that only these 36 genes are involved\nin its pathophysiology. Although no definitive method currently exists for exhaustive gene\nprioritization, it remains important to consider a broader set of BD-associated gene candidates\nand assess their overlap with differentially expressed genes (DEGs). Our original intent in using\nMAGMA was to provide a preliminary estimate of SNP-based genomic risk, and we now clarify\nthis motivation in the revised manuscript.\nIn line with the reviewer’s comment that MAGMA alone is insufficient for robust gene\nprioritization, we have now incorporated additional integrative approaches to complement the\nMAGMA findings:\n- Fine-mapping using PolyFun+SuSiE/FINEMAP\n- Transcriptome-wide association study (TWAS) using FUSION\nWe used fine-mapping to narrow down the putative BD risk genes from MAGMA analysis.\nTWAS data was not used to narrow down the putative BD risk genes, but only to support the\nfindings of CACNA1C, because TWAS data is not comprehensive across all brain regions,\nlacking\nthalamic data. These methods were employed to refine candidate BD-associated genes.\nNevertheless, we emphasize that the results from these approaches remain putative, and the\nidentified genes should be interpreted with appropriate caution, given the inherent challenges\nof mapping disease-associated genes from GWAS.\nAccordingly, we have removed the previously described list of “164 BD risk genes” derived\nfrom MAGMA. While we acknowledge the preliminary nature of MAGMA, we have retained\nits results to provide a rough estimate of SNP-based genomic risk and to address Reviewer #4’s\nsuggestion regarding the extension of MAGMA analyses to other psychiatric disorder GWAS\ndatasets. Finally, as further recommended, we have conducted and visualized separate\nanalyses for the thalamus and cortex to improve anatomical specificity.\nWe have revised the following points. In addition, throughout the manuscript, we replaced\n“BD risk genes” with “putative BD risk genes.”\nLine 223 in the revised manuscript\nUpon closer inspection of probable DEGs, MAGMA-based BD genomic risk1,46 was most\nenriched in the downregulated DEGs in Tha_ExN_CALB (P = 0.00123). Additional fine-mapping\nby PolyFun+SuSiE/FINEMAP47-49 further supported CACNA1C and SHISA9 as putative BD risk\ngenes within the downregulated set. Especially, SHISA9 is one of the 36 credible genes for BD\nin the largest GWAS to date1. While CACNA1C is not listed in these 36 credible genes, CACNA1C\nhas consistently ranked among the top genes by p-value in prior PGC1-4 studies1,2,50,51, and\nwas also implicated by our FUSION-based TWAS integration. In contrast, the upregulated\nprobable DEGs in Tha_ExN_CALB exhibited weaker enrichment for MAGMA-based BD\ngenomic risk (P = 0.0916), including no putative BD risk genes despite the greater number of\nupregulated genes (n = 22 for downregulated vs. 43 for upregulated, Supplementary Table 4).\nLine 387 in the revised manuscript\nWe examined which cell types most expressed putative BD risk genes. Thalamic excitatory\nneurons, including PVT neurons, showed the highest expression of the putative BD risk genes\nconcerning the number of expressed genes and their average expression rate (Fig. 8b). For\nexample, CACNA1C and SHISA9 exhibited relatively higher expression in thalamic excitatory\nneurons than in cortical excitatory neurons (Fig. 8c). However, the overall differences were\nmodest, and both genes showed the highest expression in Tha_ExN_PVALB_1, suggesting that\nadditional factors likely contribute to the transcriptional alterations observed in PVT neurons,\nas discussed later.\nLine 1034 in the revised manuscript\nAssessment of BD risk genes\nTo evaluate the enrichment of SNP-based genomic risk for BD, we first applied MAGMA46 to\nroughly assess the overlap between BD GWAS signals and probable DEGs, using multi-ancestry\nPGC4 BD summary statistics1. Enrichment analyses were conducted for cell clusters containing\nfive or more probable DEGs (n = 14 clusters). To assess specificity, we extended the analysis to\nGWAS summary statistics from other recent Psychiatric Genomics Consortium (PGC) datasets,\nincluding schizophrenia (SCZ)33, major depressive disorder (MDD)59, attention-\ndeficit/hyperactivity disorder (ADHD)88, post-traumatic stress disorder (PTSD)89, obsessive-\ncompulsive disorder (OCD)90, substance use disorder (SUD)91, eating disorder (ED)92, panic\ndisorder (PD)93, and autism spectrum disorder (ASD)94 (https://pgc.unc.edu/).\nBeyond assessing SNP-based genomic enrichment, we further prioritized specific BD\nrisk genes. As MAGMA alone provides only preliminary gene-level association signals, we\nadopted a multi-step approach to define “putative BD risk genes”. First, we designated the “36\ncredible genes” reported in the PGC4 BD study1 as the most reliable BD-associated genes\ncurrently available. Given the highly polygenic architecture of BD, it is unlikely that only 36\ngenes account for all the BD genetic basis. Additional BD risk gene candidates were defined as\nthose overlapping between MAGMA-significant genes and those harboring fine-mapped SNPs\nidentified by PolyFun+SuSiE/FINEMAP47-49, using the multi-ancestry PGC4 BD summary\nstatistics1. MAGMA identified 167 genes associated with BD at Bonferroni-corrected\nsignificance (adjusted P < 0.05, based on 17,487 genes with Entrez IDs). Fine-mapping was\nconducted following the developer’s pipeline, using UKBB_LD as the reference panel. Fine-\nmapped SNPs were annotated using promoter regions (PsychENCODE TSS±2kb; file:\ntss.sites.codingOnly.gencode.v19.annotation.bed), enhancer regions (PsychENCODE GRN2;\nfile: INT-14_ElasticNet_Filtered_Cutoff_0.1_cross_assembly_GRN_2.csv), and SnpEff-4.3 for\ncoding consequences. We then defined “putative BD risk genes” as the union of the 36 credible\ngenes and 75 MAGMA-based genes containing fine-mapped SNPs in promoters, enhancers, or\nas missense variants, resulting in a total of 107 genes (four overlapped, Supplementary Data\n17). To explore regional and cell-type specificity, we evaluated the expression of these 107\ngenes across thalamic and cortical cell clusters using two metrics: (1) the number of putative\nBD risk genes expressed per cluster (expression rate > 0.1), and (2) the mean expression level\nof those genes.\nTo further support gene prioritization, we conducted transcriptome-wide association\nstudies (TWAS) using FUSION95 with precomputed weights from the GTEx v8 brain datasets\n(http://gusevlab.org/projects/fusion/). Due to the absence of thalamic samples in GTEx, the\nTWAS results were used only for complementary reference and were not included in the\ndefinition of putative BD risk genes. TWAS-significant genes were identified using Bonferroni\ncorrection for each precomputed weight set (adjusted P < 0.05, across 2,557-7,732 tested\ngenes per dataset). Given that GTEx represents bulk RNA-seq data from multiple brain regions\n(not including thalamus), we interpreted these results broadly, without strict anatomical\nspecificity.\nFigure 4. Transcriptional disturbances in PVT neurons.\nc) PPI networks highlighting downregulated and upregulated probable DEGs in PVT neurons.\nNode color indicates the direction of gene regulation (Red, down; blue, up). Rectangles\nindicate BD risk genes. PPI enrichment P-value < 1.0 × 10-16. d) PPI networks highlighting\ndownregulated possible DEGs in PVT neurons. Rectangles represent BD risk genes, and colors\ndenote the genes in top DEG-GOs. Font size indicates DEG-GO enrichment. Note that the two\nPPI networks are identical but have different focuses on the genes. PPI enrichment P-value <\n1.0 × 10-16.\nSupplementary Data\n17.DataS17.PutativeBDriskgenes.txt\nPutative bipolar disorder risk genes identified in this study and their sources.\nWe have also revised the relevant phrasing (from “and” to “and/or”) to reflect the putative\nnature of BD risk genes, acknowledging that KCNQ3 is not included among the putative BD risk\ngenes.\nAbstract in the revised manuscript\nNotably, PVT neurons showed significant downregulation of genes associated with synaptic\nfunction and ion channel regulation, including SHISA9, CACNA1C, and KCNQ3—genes linked\nto putative BD risk and/or implicated as critical hubs in the downregulated molecular\nnetworks.\nLine 40 in the revised manuscript\nIn PVT neurons, SHISA9, CACNA1C, and KCNQ3 exhibited notable changes as putative BD risk\ngenes and/or critical hubs in the downregulated synaptic/neuronal gene networks.\nIn response to this comment, we examined the expression patterns of HECTD2 and AKAP11\nacross cell types in both the thalamus and cortex. Both genes exhibited relatively high\nexpression in glial populations in these regions. We did not observe consistent differential\nexpression between BD and control samples for either gene in any cell clusters in this study.\nThis lack of differential expression may reflect distinct mechanisms underlying BD-associated\ngenes identified through rare variants, as opposed to those implicated by common variant\nanalyses. We have added a brief description of these findings to the Results section and\nincluded new panels in Supplementary Fig. 8.\nLine 394 in the revised manuscript\nWe examined the expression patterns of BD-associated genes identified through rare variant\nanalyses, specifically HECTD2 and AKAP113,60. Neither gene was differentially expressed in any\ncell type clusters in the thalamus or cortex at an FDR threshold of < 0.1. The expression levels\nof HECTD2 and AKAP11 were not notably elevated in thalamic excitatory neurons, and their\nspatial expression patterns showed no specificity to the PVT region (Supplementary Fig. 8a-\nb). These findings may reflect distinct mechanisms underlying BD-associated genes identified\nthrough rare variants compared to those implicated by common variant analyses.\nSupplementary Figure 8. HECTD2 and AKAP11 expression and thalamus-cortex meta-gene\nprograms estimated by geneNMF.\na) Gene expression of HECTD2 and AKAP11 across the cell nuclei clusters in the thalamus-\ncortex integration set. b) Spatial gene expression of HECTD2 and AKAP11 in the medial\nthalamus.\no Developmental origin of thalamic inhibitory neurons\no Only MGE, CGE, and LGE were considered, but it is not convincing that the observed thalamic\ninhibitory neurons were derived from LGE\no Existing literature points to midbrain-derived (SOX14+) inhibitory neurons\nWe apologize for the careless omission of discussion regarding thalamic inhibitory neurons\nduring word reduction. As the reviewer correctly noted, thalamic inhibitory neurons are\nthought to originate from the midbrain and are characterized by SOX14 and OTX2 expression,\nas reported by Kim et al. Our original intention was to discuss lateral ganglionic eminence (LGE)\nmarker expression alongside the existing literature on SOX14 and OTX2, as illustrated in the\noriginal Supplementary Fig. 1c. However, this discussion was inadvertently removed in the\ncourse of final manuscript preparation.\nUpon reflection, we recognize that we placed excessive emphasis on the LGE, which was based\nsolely on marker expression and remains highly speculative. In response, we have renamed\nthe cluster originally labeled as “InN_LGE” to “InN_Tha” to more appropriately reflect its\ncharacteristics. We have also revised the manuscript to include relevant arguments of SOX14\nand OTX2 expression.\nLine 91 in the revised manuscript\nInhibitory neurons were classified into three distinct classes based on their source tissues and\nganglionic eminence markers compiled from Schmitz et al.34. Thalamic inhibitory neurons\nexpressed lateral ganglionic eminence (LGE) markers including ZFHX4, FOXP2, MEIS2, and\nTSHZ134, as well as midbrain-derived inhibitory markers such as SOX14 and OTX2, consistent\nwith the findings from Kim et al.35 (Supplementary Fig. 1c). These neurons were annotated as\nInN_Tha without reference to their developmental origin, which, while intriguing, is beyond\nthe scope of the present study.\nAccordingly, we have revised “InN_LGE” to “InN_Tha” throughout the main text and figures\n(Fig. 1, Supplementary Figs. 1 and 2).\nc) Top: Expression of excitatory neuron (ExN) markers and selected inhibitory neuron (InN)\nmarkers within ExN clusters. Bottom: Expression of InN markers within InN clusters. NeuCyto\ncluster is retained as a technical reference.\no Comparison with Siletti et al. is unclear, esp. regarding sample subject variation\nIn responding to this comment, we realized that two of the HBCA datasets had been\nmisassigned: the correct samples should be 10X398-5 and 10X398-6 (epithalamus), rather\nthan 10X398-1 and 10X398-2 (hypothalamus). We sincerely apologize for this oversight and\nhave revised the relevant figures and corresponding text as described above. Thanks to this\ncorrection, the characteristics of PVT neurons become clearer. Please refer to the above\ncorrections for details.\nRegarding sample-to-sample variation, we have provided the metadata for the selected 20\nsamples in the revised Supplementary Data 11. In the original manuscript, the left panel of\nFig. 2a depicted the composition of PVT neurons across our dataset and the HBCA dataset;\nhowever, the color coding was unclear. Upon review, we recognized that the composition\nacross different samples was not informative. In the revision, we instead present the thalamic\ncellular composition for each sample in our dataset in Supplementary Data 11 and visualize\nthe case-control comparison of thalamic cellular composition in Supplementary Figs. 3b-c,\nwhich is more relevant to our study. As the HBCA dataset was not intended for PVT analysis,\nthe cellular composition of PVT neurons in the HBCA samples was highly variable. For this\nreason, we removed the HBCA panels from the original Fig. 2a.\nMain Figure\nREMOVED: Figure 2a right panel in the original manuscript\nLine 146 in the revised manuscript\nWe obtained 146,643 cell nuclei in 18 unique cell clusters and 2,139 droplets in two NeuCyto\nclusters after QC (Fig. 3a and Supplementary Fig. 3a). Two NeuCyto clusters were removed as\nnon-nuclei from the subsequent analysis. Tha_ExN_CALB in Fig. 3a corresponded to PVT\nneurons. All major cell classes were consistently detected across the 41 samples\n(Supplementary Fig. 3b). The Tha_ExN_CALB cluster was also robustly identified in all samples\namong thalamic excitatory neurons with two outliers, BD13 and CT10 (Supplementary Fig.\n3c).\nLine 852 in the revised manuscript\nWe deconvoluted the Visium spatial transcriptomic data using cell-type expression profiles\nobtained from snRNA-seq via CytoSPACE39 (v1.1.0) with default parameters to estimate the\nspatial distribution of each cell type. The reference cell types included excitatory/inhibitory\nneurons, astrocytes, OPCs, oligodendrocytes, and microglia, as defined in the integrated set\ncomprising the HBCA and our dataset (see Fig. 2a). The low representation of the\nExN_PVALB_2 cluster in Fig. 2d was consistent with its low/no proportion in the CT06 snRNA-\nseq dataset (Supplementary Data 11), supporting the validity of our deconvolution approach.\nAs a robustness check, we repeated the deconvolution using either our control thalamic\ndataset or CT06-derived data alone, both of which yielded consistent results (Supplementary\nFig. 13e-f).\nSupplementary Figure 3. Cell-type annotation in the thalamus.\nb) Proportions of major cell classes across individual thalamic samples. c) Proportions of\nexcitatory neuron subtypes across individual thalamic samples.\nSupplementary Data\n11.DataS11.HBCA.profiles.txt\nSample profiles of HBCA data used in this study, including cell cluster counts.\n88\nThis manuscript investigates cellular and molecular alterations in bipolar disorder (BD) through\nsingle-nucleus RNA-seq (snRNA-seq) on postmortem brain tissue from BD patients and controls,\nfocusing on the medial thalamus and frontal cortex. The authors first identify a subpopulation\nof thalamic excitatory neurons that represent neurons of the paraventricular thalamic nucleus\n(PVT), and provide a detailed molecular and spatial characterization of this cell population.\nSubsequent analyses reveal changes in cell type compositions and gene expression in BD.\nIntriguingly, these changes are more substantial in the thalamus compared to the cortex,\nhighlighting the importance of this region and PVT neurons in particular.\nOverall, this study makes important contributions to our understanding of both the\ncytoarchitecture of the human PVT and alterations in BD at the cellular compositional and\ntranscriptional levels. In addition, the analyses of the snRNA-seq data are thorough and I\ncommend the authors for their careful assessments of QC metrics and the robustness of their\nresults. There are several considerations that if addressed would strengthen the study.\nWe thank the reviewer for their thoughtful and constructive comments, as well as for\nrecognizing the novelty and thoroughness of our study. We appreciate the positive assessment\nof our analyses, QC measures, and the robustness of our results. We are encouraged that the\nreviewer finds our work to make important contributions to understanding the human PVT\nand its alterations in BD.\nMajor comments:\n(1) Dysregulation of circadian rhythms is a feature of BD. Were any associations observed\nbetween the transcriptional changes in BD and patients' time of death?\nWe appreciate the reviewer’s comment regarding the potential influence of circadian rhythm\non transcriptional states in BD. We agree that circadian rhythm dysregulation is a well-\nrecognized feature of BD. As previously reported in schizophrenia (Seney et al., Nature\nCommunications, 2019, 10:3355), gene expression patterns in BD may also differ from those\nin controls. While the possible association between transcriptional changes in BD and time of\ndeath is intriguing, the time of death in our dataset is only an estimate and is likely to be\ninaccurate; moreover, the discrepancy between the estimated and actual time of death cannot\nbe determined. This issue is compounded by the fact that our samples were derived from\nindividuals who died by suicide. Given that circadian rhythms can influence gene expression\nwithin a 3–4-hour window, reliance on imprecise estimates of time of death could be\nmisleading, rendering any analysis of its association with transcriptional changes in BD virtually\ninfeasible. This limitation also applies to postmortem interval (PMI), and we have now\nexplicitly acknowledged it in the manuscript.\nLine 620 in the revised manuscript\nBecause the exact time of death for the donors was unavailable, the PMI in our dataset was\nnot completely accurate. Given that circadian rhythms can influence both gene expression and\n90\ndisease-associated patterns79, precise determination of the time of death would facilitate the\ninvestigation of circadian rhythm-associated gene expression changes.\n(2) The gene set enrichment analyses in this study use terms from the Gene Ontology (GO)\ndatabase. Given that synaptic dysfunction is implicated in BD pathology, I recommend the\nauthors also perform these analyses using the curated gene sets from the SynGO database\n(Koopmans et al. 2019), which were also used in O'Connell et al. (2025).\nWe thank the reviewer for this valuable suggestion. We have now performed gene ontology\nenrichment analysis using the SynGO curated resource (https://www.syngoportal.org/,\nKoopmans et al., Neuron 2019) for differentially expressed genes (DEGs) in the relevant\nclusters, including PVT neurons. We used more relaxed threshold of FDR (FDR < 0.2) than\nconventional GO enrichment analysis (FDR < 0.05), because the number of synaptic DEGs in\nour DEG criteria was limited.\nWe found six GO terms that were enriched in the thalamic clusters and none in the cortical\nclusters. Among the six SynGO terms, GO:0099055 (integral component of postsynaptic\nmembrane), GO:0099056 (integral component of presynaptic membrane), and GO:0099061\n(integral component of postsynaptic density membrane) were enriched in the downregulated\ngenes in PVT neurons. GO:1905274 (regulation of modification of postsynaptic actin\ncytoskeleton) was enriched in the downregulated genes in microglia. This analysis further\nsupports our findings of synaptic dysregulation both in PVT neurons and microglia.\nThese results are included in a new Supplementary Table 5 and mentioned in the Results and\nDiscussion. In addition to SynGO analysis, we performed GSEA in response to Reviewer #1’s\ncomment.\nLine 204 in the revised manuscript\nGiven the enrichment of synaptic terms in the DEG-GOs, we additionally performed SynGO\nenrichment analysis45. The downregulated DEGs in the Tha_ExN_CALB cluster again showed\nthe strongest SynGO enrichment, including three terms related to synaptic structure: integral\ncomponent of postsynaptic membrane (GO:0099055), integral component of presynaptic\nmembrane (GO:0099056), and integral component of postsynaptic density membrane\n(GO:0099061) (Supplementary Table 5).\nComplementary GSEA based on fold changes supported the enrichment of neuronal\nand synaptic processes among downregulated genes in Tha_ExN_CALB (Supplementary Data\n7). GSEA identified potassium channel activity, including KCNQ3, KCND2, and KCNC2 as the top\nGO term for the downregulated genes in Tha_ExN_CALB. The overlap between DEG-GO and\nGSEA results (GSEA-GO) was more pronounced for downregulated than upregulated genes\n(9/24 vs. 6/58; two-sided Fisher’s exact test, P = 0.00931, odds ratio = 5.07) (Supplementary\nTable 5). Notably, the top three DEG-GO terms among downregulated genes—transporter\ncomplex (GO:1990351), monatomic ion channel complex (GO:0034702), and main axon\n(GO:0044304)—were also identified by GSEA. The downregulated genes in Tha_ExN_CALB\nexhibit greater biological coherence and interpretability across enrichment methods.\nSupplementary Table 6 summarizes the downregulated probable DEGs indicated by both DEG-\nGO and GSEA as biologically interpretable genes for BD. Among them, KCNQ3, CACNA1C,\nSHISA9, KCNC2, LRRC7, ROBO2, and CNTNAP4 were consistently enriched across all three\nanalytical approaches (DEG-GO, GSEA, and SynGO).\nLine 285 in the revised manuscript\nDisturbed interaction between PVT neurons and microglia in BD\nIn addition to Tha_ExN_CALB (PVT neurons), we identified downregulated robust DEG-GOs\nwithin Tha_Micro_1 (microglia). Notably, these two clusters were the only ones displaying\ndownregulated robust DEG-GOs (Fig. 3e-f). Furthermore, all five DEG-GOs identified within\nTha_Micro_1 were related to synaptic function (Supplementary Table 5). When we assigned\nthe DEGs in Tha_Micro_1 to co-expression network modules, six of ten downregulated\nprobable DEGs (60%) and 32 of the 55 downregulated possible DEGs (58.2%) were assigned\nto Module-11 of microglia (Micro-M11) (P = 8.5 × 10⁻⁴ and 1.20 × 10⁻¹⁴, FC = 4.36 and 4.22,\nhypergeometric test) (Supplementary Table 4). Among the probable DEGs, SYNDIG1, CYFIP1,\nand TIAM1 emerged as promising synaptic DEGs (Fig. 5a), each exhibiting robustness scores\nexceeding 0.95 and included within the five synaptic DEG-GOs. Additional GSEA identified\nsynaptic enrichment among the downregulated DEGs in Tha_Micro_1 (Supplementary Data\n7), again highlighting SYNDIG1 and CYFIP1. SynGO enrichment analysis suggested potential\nenrichment of GO:1905274 (regulation of modification of the postsynaptic actin cytoskeleton)\n(Supplementary Table 5).\nLine 353 in the revised manuscript\nAmong cortical DEG-GO, only GO:0045211 (postsynaptic membrane), an upregulated GO term\nin FrC_ExN_L6ITCar3, was identified by FC-based GSEA, and no SynGO enrichment was\nobserved. These findings suggest that cortical DEGs exhibit lower biological interpretability\ncompared to thalamic DEGs.\nLine 489 in the revised manuscript\nAmong all the thalamic and cortical cell types in our study, PVT neurons exhibited the highest\nnumber of DEGs in BD. Notably, downregulated DEGs in PVT neurons were more biologically\ninterpretable than upregulated ones and included putative BD risk genes. Among the PVT\nneuron DEGs, SHISA9, CACNA1C, and KCNQ3 stood out as particularly relevant to BD\npathophysiology, as they were centrally positioned in the PPI network and consistently\nhighlighted in GO, GSEA, and SynGO analyses.\nLine 1025 in the revised manuscript\nTo get insight into synaptic genes, SynGO45 enrichment analysis was added to the clusters with\nrobust DEG-GOs, using enricher function in clusterProfiler with the FDR threshold of 0.2 and\ngene set size of 5 – 100 genes. SynGO terms with at least two genes were retained to\naccommodate the smaller number of target genes compared to the above DEG-GO analysis.\nAs we set FDR < 0.2, the SynGO enrichment should be regarded as a potential one. We also\nperformed FC-based GSEA separately for the downregulated and upregulated gene sets to\nfurther evaluate the robustness of the identified biological signals, using gseGO function in\nclusterProfiler with the enrichmentScore threshold of 0.7.\nSupplementary Table 5 (abridged)\n(3) The cohort includes both BDI and BDII subtype cases. Is the dataset sufficiently powered to\nstratify analyses by subtype? I understand that the sample size of the BDII cases may preclude\nsuch analyses; if so, it would be worth mentioning this in the limitations section of the\nDiscussion (together with sex differences in lines 449-451).\nWe agree that differences between BD subtypes may be biologically meaningful. However, our\ncurrent cohort includes only five individuals with BDII, limiting the statistical power for\nsubtype-stratified analysis. We now explicitly acknowledge this in the revised Discussion (in\nthe limitations paragraph) and have included subtype information in Supplementary Table 1.\nIn addition to addressing the reviewers’ comments, we have revised Supplementary Table 1\nto better protect individual privacy. The original version included individual-level demographic\nand clinical data. Because Coroner data in Canada, where the brain tissue specimens were\ncollected, are nominal and publicly accessible, there is a theoretical risk of re-identification.\nIn\nthe revised manuscript, we now provide aggregated demographic and clinical information,\nincluding a summary of subtype information, instead of individual-level data.\nLine 616 in the revised manuscript\nWhile the differences between BD subtypes (BD-I and BD-II) are of interest, the limited number\nof BD-II cases (N = 5) in our cohort also precludes a reliable assessment of subtype-specific\neffects.\nSupplementary Table 1 (abridged)\n(4) In Fig. 8a, downregulated genes in Tha_ExN_CALB are enriched in SNP-based genomic risk\nof all 3 phenotypes tested. Is this observed for other psychiatric disorders (such as post-\ntraumatic stress disorder, attention deficit-hyperactivity disorder, borderline personality\ndisorder, and autism spectrum disorder)?\nWe appreciate the reviewer’s suggestion. We conducted additional MAGMA enrichment\nanalyses using publicly available GWAS summary statistics for attention-deficit/hyperactivity\ndisorder (ADHD), post-traumatic stress disorder (PTSD), obsessive-compulsive disorder (OCD),\nsubstance use disorder (SUD), eating disorder (ED), panic disorder (PD), and autism spectrum\ndisorder (ASD). These analyses revealed no significant enrichment of SNP-based genomic risk\nin the downregulated genes of PVT neurons for any of these additional disorders. In contrast,\nenrichment was observed for the GWAS data of bipolar disorder (BD), schizophrenia (SCZ), and\nmajor depressive disorder (MDD). While this could be partly derived from the sample size of\nthe original GWAS, where BD, SCZ, and MDD have larger sample sizes, notable enrichment of\nthese associations implicates a possible relevance of PVT neurons to these disorders. To\nillustrate this clearly, we have revised Fig. 8. The full association data are described in\nSupplementary Data 9.\nUnfortunately, GWAS summary statistics for borderline personality disorder (BPD) are not yet\npublicly available. According to the preprint by Streit et al. (medRxiv, 2024,\nhttps://www.medrxiv.org/content/10.1101/2024.11.12.24316957v1), the summary statistics\nare expected to be released upon official publication. We will incorporate the BPD analysis\nonce these data become available.\nLine 379 in the revised manuscript\nThis enrichment pattern extended beyond BD, as downregulated DEGs in PVT neurons also\nshowed the highest enrichment for MAGMA-based genomic risk associated with\nschizophrenia33 and major depressive disorder59, while showing weaker or nonsignificant\nenrichment for other psychiatric conditions: attention-deficit/hyperactivity disorder, post-\ntraumatic stress disorder, obsessive-compulsive disorder, substance use disorder, eating\ndisorder, panic disorder, and autism spectrum disorder (Fig. 8a, details in Supplementary Data\n9). These enrichment comparisons should be interpreted with caution, as differences in\nstatistical power across the original GWAS datasets may influence the observed patterns.\nLine 1038 in the revised manuscript\nTo assess specificity, we extended the analysis to GWAS summary statistics from other recent\nPsychiatric Genomics Consortium (PGC) datasets, including schizophrenia (SCZ)33, major\ndepressive disorder (MDD)59, attention-deficit/hyperactivity disorder (ADHD)88, post-\ntraumatic stress disorder (PTSD)89, obsessive-compulsive disorder (OCD)90, substance use\ndisorder (SUD)91, eating disorder (ED)92, panic disorder (PD)93, and autism spectrum disorder\n(ASD)94 (https://pgc.unc.edu/).\na) Left: MAGMA-based genomic risk enrichment in DEGs for BD, schizophrenia (SCZ), and\nmajor depressive disorder (MDD). Right: MAGMA-based genomic risk enrichment in\ndownregulated DEGs in PVT neurons for BD, SCZ, MDD, attention-deficit/hyperactivity disorder\n(ADHD), post-traumatic stress disorder (PTSD), substance use disorder (SUD), obsessive-\ncompulsive disorder (OCD), eating disorder (ED), panic disorder (PD), and autism spectrum\ndisorder (ASD).\nSupplementary Data\n09.DataS09.MAGMA_PsychiatricDisorders.txt\nMAGMA results for selected clusters and 10 psychiatric disorders.\n(5) In addition to the analyses of common variant-based BD genomic risk, it would also be\nhelpful to comment on how these findings relate to the rare loss-of-function variants recently\ndescribed in Thorgeirsson et al. (2025). For example, which cell types in thalamus and cortex\nexpress HECTD2 and AKAP11, and how does their expression in BD patients compare to\ncontrols?\nWe thank the reviewer for raising this important point. We examined the expression patterns\nof HECTD2 and AKAP11 across cell types in both the thalamus and cortex, as these two genes\nwere associated with BD in Thorgeirsson et al. (Nature Genetics 2025, 57:851-855). Both genes\nexhibited relatively high expression in glial populations in these regions. We did not observe\nconsistent differential expression between BD and control samples for either gene in any cell\n95\nclusters in this study. This lack of differential expression may reflect distinct mechanisms\nunderlying BD-associated genes identified through rare variants, as opposed to those\nimplicated by common variant analyses. We have added a brief description of these findings\nto the Results section and included new panels in Supplementary Fig. 8.\nLine 394 in the revised manuscript\nWe examined the expression patterns of BD-associated genes identified through rare variant\nanalyses, specifically HECTD2 and AKAP113,60. Neither gene was differentially expressed in any\ncell type clusters in the thalamus or cortex at an FDR threshold of < 0.1. The expression levels\nof HECTD2 and AKAP11 were not notably elevated in thalamic excitatory neurons, and their\nspatial expression patterns showed no specificity to the PVT region (Supplementary Fig. 8a-\nb). These findings may reflect distinct mechanisms underlying BD-associated genes identified\nthrough rare variants compared to those implicated by common variant analyses.\nSupplementary Figure 8. HECTD2 and AKAP11 expression and thalamus-cortex meta-gene\nprograms estimated by geneNMF.\na) Gene expression of HECTD2 and AKAP11 across the cell nuclei clusters in the thalamus-\ncortex integration set. b) Spatial gene expression of HECTD2 and AKAP11 in the medial\nthalamus.\n(6) While the data that support potential disruptions in PVT-microglia trans-cellular\ninteractions in BD are compelling (MP7, Figs. 5b-c), it was unclear to me why the authors\nfocused on SHISA9 in PVT neurons as a representative MP7 gene throughout the manuscript.\nAlthough SHISA9 is part of the core PPI network of BD-associated downregulated genes in PVT\nneurons (Fig. 5d), it appears to be expressed in only in a small fraction of Tha_ExN_CALB (PVT)\nneurons, while is expressed in most Tha_ExN_PVALB_1 neurons (Figs. 5d and 3a). I suggest\nthat the authors either characterize this SHISA9-positive subpopulation of PVT neurons further\n(for example, are these detected in all thalamic samples?), or reduce/remove references to\nSHISA9 in PVT neurons throughout the text. These include parts of the section on \"Disturbed\ninteraction between PVT neurons and microglia in BD\" (lines 251-284), the Discussion (lines\n415-417), as well as the following sentence in the abstract: \"This interaction was further\nexemplified by the concurrent downregulation of SYNDIG1 in microglia and SHISA9 in PVT\nneurons, both of which are involved in AMPA receptor signaling pathways\" (top of pg. 3, line\nnumbers not available).\nremove certain overstatements. The previously described “interaction between SHISA9 and\nSYNDIG1” represented an over-interpretation and has therefore been removed. Instead, we\nnow describe SHISA9 and SYNDIG1 more conservatively, noting their roles as auxiliary subunits\nof AMPA receptors (Kalashnikova et al. Neuron 2010, 65:80-93; Khodosevich et al. Neuron\n2014, 83:601-615).\nOur focus on SHISA9 among the downregulated genes in PVT neurons is primarily motivated\nby its identification as one of the 36 credible BD-associated genes reported in O’Connell et al.,\nNature (2025). In addition, we noticed that SYNDIG1, the top-ranked downregulated gene in\nmicroglia, shares the property of encoding auxiliary subunits of the AMPAR complex with\nSHISA9, in parallel with probable interactions between microglia and thalamic excitatory\nneurons, as estimated by GeneNMF approach. The biological interpretability of SHISA9 and\nSYNDIG1 was further reinforced by GSEA (as suggested by Reviewer #1) and SynGO analysis\n(as suggested by Reviewer #4). The DEGs as convergence of DEG-GO and GSEA are illustrated\nin a new Supplementary Table 6, which includes SHISA9 and SYNDIG1. In this revisional\nprocess, we have also noticed that TIAM1, one of downregulated genes in microglia as\nbiologically interpretable genes in Supplementary Table 6, might have roles in regulating\nsynaptic AMPAR function (Rao et al. J Neurosci 2019, 39:9306-9315).\nWe have also revised the visualization of SHISA9 expression. In the original version of Fig. 5d,\nSHISA9 and SYNDIG1 were displayed together. Due to limitations in Seurat's default color\nscaling, we were unable to independently adjust the color range for each gene, which led to a\nmisleading color annotation for SHISA9. In the revised manuscript, SHISA9 and SYNDIG1 are\npresented in separate panels (Fig. 5d). As shown in the updated figure, SHISA9 was broadly\nexpressed in Tha_ExN_CALB (PVT neurons), while the strongest expression was observed in\nTha_ExN_PVALB_1. Additionally, we included detailed expression profiles of SHISA9, alongside\nCACNA1C, to highlight that SHISA9 is highly expressed in thalamic excitatory neurons relative\nto cortical excitatory populations in Fig. 8c. This figure further illustrates high expression of\nSHISA9 in Tha_ExN_CALB (PVT neurons).\nWe selected SHISA9 and SYNDIG1 for the reasons outlined above. Although the meta-gene\nprogram encompassing thalamic excitatory neurons and microglia requires further\nexperimental validation, the concurrent downregulation of several genes involved in AMPAR\nregulation represents a noteworthy finding. In the original manuscript, the visualization of\nSHISA9 expression in PVT neurons was misleading; however, SHISA9 is in fact highly expressed\nin these neurons. We apologize for this earlier misrepresentation.\nAbstract in the original manuscript\nREMOVED: This interaction was further exemplified by the concurrent downregulation of\nSYNDIG1 in microglia and SHISA9 in PVT neurons, both of which are involved in AMPA receptor\nsignaling pathways.\nLine 279 in the original manuscript\nREMOVED: The coordinated downregulation of SYNDIG1 in microglia and SHISA9 in excitatory\nneurons exemplifies a probable disturbance in microglia-neuron interactions for synapses in\nthe BD thalamus, indicating a plausible pathological mechanism of SHISA9, one of the 36\nputative BD-causal genes1.\nLine 415 in the original manuscript\nREMOVED: The downregulation of SYNDIG1 in microglia and SHISA9 in PVT neurons might be\nan example of dysregulated molecular interactions within the Excit-Micro meta-gene programs\nassociated with BD.\nLine 219 in the revised manuscript\nSupplementary Table 6 summarizes the downregulated probable DEGs indicated by both DEG-\nGO and GSEA as biologically interpretable genes for BD. Among the downregulated probable\nDEGs, KCNQ3, CACNA1C, SHISA9, KCNC2, LRRC7, ROBO2, and CNTNAP4 were consistently\nenriched across all three analytical approaches: DEG-GO, SynGO, and GSEA.\nLine 311 in the revised manuscript\nSupplementary Table 6 summarizes the downregulated genes in Tha_ExN_CALB and\nTha_Micro_1 that were enriched in one or more GO terms in both DEG-GO and GSEA analyses.\nNotably, these genes belonged to the same co-expression modules in each cluster.\nLine 314 in the revised manuscript\nAt the gene level, SYNDIG1 and SHISA9 encode auxiliary subunits for AMPA receptor (AMPAR)\ncomplex56,57, with both genes’ robustness scores exceeding 0.7 (Supplementary Table 4).\nNotably, SYNDIG1 was prominently expressed in microglia, while SHISA9 was primarily\nexpressed in excitatory neurons (Fig. 5d). SHISA9 consisted of core PPI networks of BD-\nassociated downregulated genes in PVT neurons (Fig. 4d) and overlapped in the\ndownregulated genes in Tha_ExN_PVALB_1, as one of the 36 putative BD-causal genes1.\nBeyond microglia and excitatory neurons, thalamic T cells were included in this meta-gene\nprogram. Fig. 5e illustrates potential ligand-receptor interactions between neurons and\nmicroglia/T cells, as analyzed through LIANA58. Although we did not detect concurrently\ndownregulated DEGs as ligand–receptor pairs, KCND2 in Tha_ExN_CALB, LAMB1 in\nTha_ExN_PVALB_2, and KCNQ1 in Tha_Micro_1 were identified as downregulated DEGs\namong the LIANA-based ligand-receptor pairs. These genes encode ion channels or\nextracellular matrix proteins that may potentially contribute to Excit–Micro interactions.\nHowever, as this analysis is an inference, the precise ligand–receptor mechanisms remain to\nbe elucidated in future studies.\nLine 556 in the revised manuscript\nBeyond cell-type-based analysis, we identified trans-cellular meta-gene programs between\nexcitatory neurons and microglia (Excit-Micro meta-gene programs), likely reflecting\nphysiological synaptic interactions between these cell types. Microglia are known to modulate\nsynaptic plasticity and neuronal circuitry through dynamic interactions with neurons76.\nNotably, such meta-gene programs in the thalamus are disrupted in BD, further characterized\nby a significant decrease in synaptic gene expression in both PVT neurons and microglia. Ling\net al. reported downregulation of the synaptic neuron and astrocyte program in\nschizophrenia10; our findings highlight downregulation of the synaptic genes in the thalamic\nneuron and microglia in BD. The concurrent downregulation of SYNDIG1 in microglia and\nSHISA9 in excitatory neurons suggests disruptions of AMPAR signaling, as both genes encode\nauxiliary subunits of the AMPAR complex56,57. TIAM1, another downregulated gene in\nmicroglia, has also been implicated in regulating synaptic AMPAR function77. The preferential\nlocalization of SYNDIG1-positive microglia in cortical gray matter is consistent with this\ninterpretation. Nonetheless, these meta-gene programs remain data-driven constructs, and\nthe roles of SYNDIG1 in microglia—as well as their mechanistic links to SHISA9 in PVT\nneurons—are currently unknown. Future studies employing animal and in vitro cellular\nmodels will be crucial for elucidating these mechanisms.\nFigure 5. Disturbed interaction between microglia and PVT neurons in the thalamus.\nd) UMAP projection illustrating the expression patterns of SYNDIG1 and SHISA9.\nFigure 8. Comparative analysis between the thalamus and cortex.\nc) Relative expression of CACNA1C and SHISA9 across all clusters in the thalamus-cortex\nintegration dataset.\nbiological interpretations\nMinor comments:\n(1) Please provide details on the cohort in the Methods, including information on how\nindividuals were selected for inclusion and the approval of the study protocol. Details are in the\nReporting Summary, but must also be provided in the manuscript.\nWe have revised the Methods section to include additional information including cohort\nselection, inclusion criteria, and ethics approval.\nLine 634 in the revised manuscript\nPost-mortem brain samples\nPostmortem brain specimens were sourced from the Douglas–Bell Canada Brain Bank\n(www.douglasbrainbank.ca). Frozen samples were obtained from 21 BD patients and 21\nmatched controls (Supplementary Table 1), with written informed consent provided by the\nnext of kin. Due to the limited availability of postmortem thalamic tissue, no formal sample\nsize calculation was performed; instead, we included all eligible samples available at the time\nof collection, ensuring sex matching between BD cases and controls. Each group comprised 14\nmales and 7 females prior to quality control. The average age at death was significantly higher\nin the control group than in the BD group (52.5 vs. 44.3 years; two-sided t-test, P = 0.0472),\nreflecting our prioritization of paraventricular thalamic region availability and sex matching\nover age matching. The cause of death was significantly skewed toward suicide in the BD group\n(two-sided Fisher’s exact test, P = 2.33 × 10-8), consistent with the elevated suicide risk in BD\n100\nand the characteristics of this brain bank. Of the 42 donors, 41 were Caucasian and one Asian.\nDonors exhibited substantial heterogeneity in substances present at death, precluding\nstratification into subgroups of more than two individuals. Summary statistics of substances\npresent at death are presented in Supplementary Fig. 10a.\nLine 670 in the revised manuscript\nThis research was conducted in accordance with the Helsinki Declaration and received\napproval from the Research Ethics Committee of the Faculty of Medicine, Juntendo University\n(approval ID: M19-0278).\n(2) The final two paragraphs of the Discussion appear to be redundant.\nThank you for pointing this out. We have removed the final redundant paragraph of the\nDiscussion and streamlined the concluding section to avoid repetition.\nLine 452 in the original manuscript\nREMOVED: Despite the limitations of this study, our findings suggest that disturbances in PVT\nneurons constitute a central pathology in BD. Neuroimaging techniques targeting the PVT can\nenhance clinical assessments of BD. The presence of PVT-like neurons in thalamic organoids\nhighlights the feasibility of modeling PVT neurons from patient-derived peripheral samples.\nThis approach could facilitate the discovery of pharmacological interventions that target the\nmolecular mechanisms underlying BD. The integration of multi-modal methodologies,\nincluding cellular, molecular, and neuroimaging techniques, will offer a more comprehensive\nunderstanding of the pathological role of PVT neurons in BD. PVT neurons emerge as a\npromising target for the development of novel diagnostic and therapeutic strategies for\npatients with BD.\nWe have revised Supplementary Table 1 to better protect individual privacy. The original\nversion included individual-level demographic and clinical data. Because Coroner data in\nCanada, where the brain tissue specimens were collected, are nominal and publicly accessible,\nthere is a theoretical risk of re-identification. In the revised manuscript, we now provide\naggregated demographic and clinical information, including a summary of medication\nexposure, instead of individual-level data.\nSupplementary Table 1 (abridged)\nIn responding to Reviewer #1’s comment, we recognized that the original legend for Fig. 1d\nwas insufficient, and we have revised it to improve clarity.\nFigure 1. Characterization of thalamic excitatory neurons.\nd) Dot plot displaying selected neuronal markers across 20 neuronal cell clusters (see\nSupplementary Fig. 1 for details). Dot size represents the proportion of cells expressing the\nmarker, and color intensity indicates average expression. Colored bars indicate major cell\nclasses as in panel c.\n- Revising “hour” to “h” as SI style description throughout the main text.\n- Add company information such as “Illumina, San Diego, USA” and “10x Genomics,\nPleasanton, USA”.\n103"}]

Decision letter