Eight different cell types detected in the SN and genes altered by PD.
For snRNA-seq analysis, a total of 33,293 nuclei were extracted from the SN of a total of 8 donors, with 4 donors from the CON group and 4 from the PD group, followed by FACS, and snRNA-seq procedure of the 10x Chromium platform (Fig. 1A and Supplementary Table 1). Visium spatial gene expression study was performed using two capture areas, namely the “medial” (A1) and “lateral” (D1) regions, from one midbrain FFPE section of a PD sample (Fig. 1A, sFig. 1A and Supplementary Table 1). Immunofluorescence image stained for phosphorylated α-Syn (green), neuronal marker (HuC/D and NeuN, red), and DAPI, was aligned with the capture areas to confirm the spatial locations of the Visium spots. The distribution of TH and Calretinin immuno-reactivity around the SN was also verified in adjacent sections (sFig. 1B).
First, in snRNA-seq data, eight different cell types formed clusters (Fig. 1B), and the integrated CON and PD samples overlapped well on the UMAP 2D space (Fig. 1C). There were no significant differences in proportions of cell types between CON and PD (Fig. 1D). However, fewer cells were obtained from the PD tissue compared to the CON tissue despite using a similar amount of tissue (total CON nuclei = 22631, total PD nuclei = 9129), and this pattern was consistently observed across all cell types (Fig. 1E).
The predominant cell types observed were oligodendrocytes, astrocytes, microglia and oligodendroglial precursor cells (OPCs). Additionally, a small population of MRC1-positive brain-resident macrophages was identified within a microglia cluster (Fig. 1B, F, Supplementary Table 2_1) [57, 58]. The neurons were organized into several distinct small clusters, and a cluster with characteristics of T cells was also detected, expressing high levels of CD44, known to be involved in cell trafficking into the central nervous system (CNS) [59, 60]. Vascular cells were clearly divided into clusters of endothelial cells and mural cells including pericytes, based on known markers such as CLDN5 and PDGFRB [61, 62].
We examined the DEGs upregulated or downregulated in PD compared to CON for each cell type (Fig. 1G and Supplementary Table 2_2). We found that the enriched gene ontology (GO) terms and KEGG pathways were well separated by both cell type and PD, indicating cell type-specific functions that may be compromised in disease circumstance or cell type-specific cellular responses (Fig. 1H). For example, downregulated genes in OPCs were linked to processes such as neuron projection development, neuron maturation, chondroitin sulfate proteoglycan metabolism contributing as a major component of the extracellular matrix (ECM), or axonogenesis-related genes implying axonopathy [63–65]. On the other hand, upregulated genes were associated with neurodegenerative environmental factors, including hypoxia and inflammatory cytokines.
We next investigated whether these transcriptomic signatures observed in our data were also associated with genetic risks previously studied in PD patients. Of the DEGs between CON and PD for each of the cell types in Fig. 1G, we compiled a list of PD-associated genes that had been identified in GWASs (including Korean GWAS), multi-omics, and genetic databases to see if they showed differential expression in our data [14, 16, 20, 40–46]. We found that many of the DEGs were associated with specific cells, with changes in OPCs, astrocytes and microglia being particularly prominent. This confirms that the DEGs we found recapitulate changes that are consistent and common in PD, particularly in glial cells (sFig. 2 and Supplementary Table 2_3).
To elucidate the spatial distribution of cell types in the SN of PD patients, unsupervised clustering on Visium spots resulted in 12 clusters (Fig. 1I and sFig. 3). Despite being initially clustered without spatial information, the Visium spots showed a clear spatial pattern that distinguished the SN and surrounding midbrain anatomical features. The top-ranked spatial DEGs included several cell type markers, allowing each cluster to be distinguished by enriched cell types (Fig. 1J, sFig. 3C, Supplementary Table 3).
Specific spatial clusters (SC) 0 and SC3 were identified as oligodendrocyte-rich regions, expected to be white matter. SC10, SC11, and SC6 were glial-rich regions representing progressive gliosis following neuronal loss. SC4 and SC9 are vascular cell-rich regions. SC7 and SC8 were neuron-rich clusters. Among the two neuron-rich clusters, SC7 had more pronounced expression of various neuronal markers than SC8, including SYT1, SLC18A2, GAD1, GAD2, CALB1 and TAC1 (also known as Substance P) (Fig. 1J and sFig 3C, 4A-D, 4I).
In particular, SC4 was a mixture of different cells, including oligodendrocytes, astrocytes, and vascular cells with a pronounced expression of PVALB (sFig. 3C, 4I), exactly representing the known features of the red nucleus (RN). All these results indicate that cell types in the SN of PD patients are organized in a spatially structured manner.
Neuronal Heterogeneity in the SN and Population Changes in PD
To observe changes that would be better detected in subpopulations of neurons in the SN, clustering was performed independently on 722 neuronal nuclei, resulting in a total of 11 neuronal subclusters (NC) (Fig. 2A). Examination of expressed neurotransmitters identified NC3 and NC8 as dopaminergic neurons expressing SLC18A2, SLC6A3, ALDH1A1 and TH. NC7, NC9, NC10, and NC11 were identified as GABAergic neurons predominantly expressing SLC6A1, GAD2, GAD1. NC1, NC2, and NC6 were found as glutamatergic neurons expressing SLC17A7 (Fig. 2B and Supplementary Table 4).
In a LIGER analysis [39] integrating our data with the transcriptional profile obtained by enrichment of human midbrain neurons recently reported by Kamath et al.[25], we confirmed that the neurons in our data retained the general features of their high-throughput profile of midbrain neuronal subtypes. Specifically, NC3 and NC8 matched well with dopaminergic neuronal subtypes expressing SOX6 and CALB1 (sFig. 5).
By dividing the proportion of identified neuronal subclusters between CON and PD groups in our dataset, we observed a reduction in dopaminergic (NC3 and NC8) and glutamatergic (NC1, NC2 and NC6) neurons, leading to an increased relative proportion of all GABAergic (NC7, NC10, and NC11) neurons, except for the NC9 subcluster, which co-expresses PVALB in the context of PD (Fig. 2C) [66].
To examine the proportional changes in neuronal subclusters observed in scRNA-seq result in the tissue, we performed RNAscope using the probes targeting SLC18A2 as a marker for dopaminergic neurons and GAD2 as a marker for GABAergic neurons (Fig. 2D). In the CON brain, SLC18A2 staining was strong in all the pigmented SN neurons while GAD2-positive staining was observed in non-pigmented neurons (Fig. 2D, top panel). These GABAergic neurons, known to be predominantly located in the substantia nigra pars reticulata (SNpr) in rodents [67, 68], were sporadically scattered and surrounding the SNpc in human samples. On the other hand, in advanced PD, where TH immunoreactivity was significantly reduced, SLC18A2-positive neurons had largely disappeared, except for a few remaining in the dorsomedial part of the SN, indicating the loss of dopaminergic neurons (Fig. 2D, middle panel and sFig. 6). GAD2-positive neurons were mainly preserved in number relative to SLC18A2-positive neurons, indicating the persistence of GABAergic neurons, which may explain the hypokinetic PD phenotypes (Supplementary Table 5) [69, 70].
In patients with confirmed synucleinopathy but no manifestation of clinical symptoms, namely preclinical PD, we found that there were still many neurons positive for SLC18A2 and GAD2, similar to those observed in CON. However, the intensity of immunostaining for TH, an enzyme critical for dopamine synthesis, was significantly decreased, which could serve as a good indication of early PD progression (Fig. 2D, bottom panel and sFig. 6).
Upon examining the enriched genes in each neuronal cluster, we were able to gain a more detailed understanding of the specific characteristics of neuronal subclusters (Fig. 2E and Supplementary Table 4). For example, among the dopaminergic neurons, NC3, which was highly reduced in number in PD, had high expression of TH, THY1 and RAB3A, while NC8, which was relatively less reduced in number in PD, had high expression of TTC6, KLHL1, and SEMA6D (Fig. 2E and sFig. 7). Different groups showed unique gene expressions. Among excitatory glutamatergic neurons, NC1 had high PAX6 expression, while in inhibitory GABAergic neurons, NC7 had high serotonin receptor 5-HT-2C (HTR2C) levels. Some groups, like NC4, shared features with other cell types, showing oligodendrocyte markers such as TF, ST18, and MOBP. NC5 specifically expressed MATN2, a gene known to be active in white matter neurons (Fig. 2E) [71].
To investigate the spatial distribution of neuronal subtypes, we deconvoluted the spatial data into neuronal subtype proportions using Cell2location (sFig. 8A) [48] followed by correlation analysis. This analysis revealed two spatially correlated groups of neuronal subtypes (Figs. 1I-J, 2F, and sFig. 7). The first group, including two glutamatergic subtypes (NC1 and NC6), two dopaminergic subtypes (NC3 and NC8), and three GABAergic subtypes (NC7, NC9, and NC11). The subtypes were mainly located in the region expected to be the dopaminergic nuclei in the midbrain, including particularly two neuron-rich SC7 and SC8 (Fig. 2F-G and sFig. 8).
The second group comprised a glutamatergic subtype (NC2), a GABAergic subtype (NC10), and two unknown neuron types (NC4 and NC5). These subtypes were primarily located in regions outside SN, particularly in SC4 (Fig. 2G and sFig. 8). This spatial separation of the two groups of neuronal subtypes suggests a transcriptomic distinction between the dopaminergic nuclei and other subregions of midbrain neurons, highlighting the unique molecular and spatial characteristics of different neuronal populations within the midbrain.
Next, we explored the genes showing expression differences between CON and PD samples within each neuronal cluster (Fig. 2H, sFig. 9A-D, and Supplementary Table 4). There was a greater number of DEGs between CON and PD within specific neuronal subtypes (sFig. 9A-B) compared to the total neuronal population (Fig. 1G), suggesting subtype-specific molecular changes associated with PD. With the exception of NC6 and NC9, where the cell number was not sufficient to perform DEG analysis due to the significant reduction in cell number caused by PD, we observed DEGs in the remaining clusters with diverse terms. DEGs related to protein processes in endoplasmic reticulum shown in NC7 and microtubule depolymerization in NC1, involved in different cellular processes (sFig. 9A-D) [72–74]. Among the most significant DEGs in each cluster (Fig. 2H), VDAC3, known to recruit Parkin and induce mitophagy [75], PLAAT3, a lipid-modifying enzyme involved in organelle degradation [76, 77], and ADAM22, involved in synaptic maturation [78], showed specific increases in dopaminergic NC3. In NC4, genes related to metal ion homeostasis such as SELENOP, a selenium transporter known for its protective role in Alzheimer's disease (AD) and PD [79, 80], and the copper transporter SLC31A2 [81] were downregulated. In NC7, SEMA3A and PLXNA1, involved in semaphorin signaling and previously implicated in neurodegenerative diseases [82, 83], were upregulated. In NC8, specific downregulation of genes such as RGS6, which is important for the maintenance of dopaminergic neurons in the ventral SN [84], and DCC, which plays a critical role in the development and survival of the ventral midbrain dopamine circuit [85–88], were observed. HSP90AA1, known as a stress sensor in neurons, showed a tendency to decrease in multiple clusters affected by PD, including NC2,3,4 and 8 [89, 90]. These subcluster-specific changes were consistent with the heterogeneous nature of PD effects on neurons, as indicated by previous ontology results (sFig. 9C, D).
We incorporated RNA velocity analysis to infer the dynamical changes in the neuronal transcriptome (sFig. 9E and Supplementary Table 6). Significant changes in the trajectory of cell development were observed when comparing healthy controls to PD patients. These changes were especially noticeable in the remaining dopamine-producing neurons (NC8). The dot product calculation between CON and PD velocity vectors for NC8 showed significant negative values (one-sided Wilcoxon signed-rank test, p = 2.64 x 10− 8), indicating alterations in both unspliced and spliced form of mRNA due to PD (sFig. 9F).
To further understand the regulators of these changes, SCENIC analysis was used to identify regulon activities in each neuronal subtype (Fig. 2I and Supplementary Table 7) [52]. Regulon activities were primarily clustered according to neuronal subtypes, but PD led to specific alterations, particularly in dopaminergic neurons NC3 and NC8 (Fig. 2I, dotted box). In NC3, regulons such as BHLHE41 and BHLHE40, circadian regulators, NFE2L1, which is important for dopaminergic neuronal survival, and TRIM28, which simultaneously regulates nuclear accumulation of α-Syn and tau, were significantly inactivated in PD [91, 92]. In NC8, PBX1, involved in protection against oxidative stress and is implicated in midbrain dopaminergic neuronal development and PD, was a specific regulon and significantly inactivated by PD [93]. SOX6, a marker gene for the division of SNpc into ventral and dorsal, was highly active in NC8 [94], suggesting distinct populations of dopaminergic neurons with differential susceptibility to PD. Taken together, despite the detected population-level changes in different neuronal subtypes, the transcriptomic changes in the neurons are likely driven by alterations in dopaminergic neurons. This finding highlights the critical role of dopaminergic neurons in the disease progression of PD as a pathological drive.
Key Transcriptomic Changes Driving PD-associated Astrocyte Trajectories
In astrocytes, we clustered cells with Seurat as in neurons and found that the clusters were not independently separated on UMAP and that transcriptomic features differed significantly between CON and PD rather than between clusters (sFig. 10). Therefore, we performed trajectory analysis on astrocytes to track pseudo-temporal expression changes as the disease progressed. The 6,128 astrocyte nuclei (Fig. 3A) were clustered into 8 astrocyte components (ACs) (Fig. 3B and Supplementary Table 8_1), resulted in two distinct astrocyte trajectories, leading to PD-enriched astrocyte components AC8 (Trajectory 1) and AC6 (Trajectory 2) (labels in Fig. 3B, Fig. 3C). Based on RNA velocity prediction [50] and GO enrichment analysis, we hypothesized that both trajectories start at AC3, which expressed genes associated with the most mature functions of normal astrocyte, including calcium signaling, cAMP signaling and long-term potentiation (Fig. 3D-E and sFig. 11) [95–98]. We visualized the expression changes and corresponding RNA velocity of the DEGs identified in the astrocytes cluster shown in Fig. 1G in a large pseudotime stream consisting of AC3-AC4-AC5-AC6-AC8 and observed that they showed a gradual change along a trajectory that was consistent in both analysis methods (Fig. 3F and Supplementary Table 8_1, 9). This also indicates that the pseudotime trajectory exactly represents the PD progression.
Among those genes, we validated expression changes in PD patient tissues, focusing on CHI3L1 and CD44, genes that have been reported to show expression changes associated with the neuronal injury, particularly the former in acute and the latter in chronic responses (Fig. 3G and sFig. 12A) [99–102]. It was observed that CHI3L1 was significantly upregulated in preclinical PD cases, particularly in SN, and cells displaying strong CHI3L1 expression were readily adjacent to pigment neurons. Additionally, CD44 was consistently co-expressed with CHI3L1 in preclinical PD cases, although its expression did not show a significant increase. Interestingly, in PD cases, we observed the expression of CD44 exceeded that of CHI3L1. These findings indicate that genes exhibiting changes in expression along the astrocyte trajectory could potentially serve as histological markers of reactive astrocytes not only in the end stage but also in the early response.
Next, to explore the differences between the two trajectories, we performed analysis to identify DEGs altered by pseudotime using a likelihood ratio test based on pseudotime spline fitting implemented in Monocle2 and observed transcriptomic changes associated with Trajectory 1 and Trajectory 2 (Fig. 3H and Supplementary Table 8_2). The result showed that 1,663 genes were shared among 1,950 DEGs from Trajectory 1 and 3,426 DEGs from Trajectory 2. Although approximately 85% of DEGs in Trajectory1 shared with Trajectory 2 across the disease trajectory, the specific terms still appeared (Fig. 3H, upper panel, and all gene lists for each Gene set are shown in Supplementary Table 8_3). Trajectory 1 was characterized by increased expression of genes related to metal ion homeostasis (Gene set 6) and glutamatergic excitatory synapses (Gene sets 7–8).
Trajectory 2 had 1,763 DEG genes that do not appear in Trajectory 1 (Fig. 3H, lower panel and Supplementary Table 8_3). This Trajectory 2 exhibited upregulation of genes involved in the p38 MAPK cascade and cellular response to general stress during the AC3 to AC4 transition (Gene set 5), potentially in response to extracellular tau or α-Syn accumulation [103–106]. Also, Trajectory 2 showed unique features, including upregulation of genes related to ribosome biology or coronavirus disease (Gene set 7) at the end state.
DEG analysis showed downregulated genes in AC6 compared to AC8, including MAGI2, CADM2 and NEBL (cell-cell interaction molecules), ERBB4, and SLC39A11 (metal ion transporter) (Fig. 3I and sFig. 12B). Genes upregulated in AC6 compared to AC8 included several ribosomal proteins (RPL34, RPS24, RPL41, RPS23, and RPS25) (Fig. 3J and sFig. 12C). Genes that are highly expressed in AC8 opposed to AC6 included cilia-related genes (DNAH5 and CABCOCO1) and the subunits of mitochondrial ATP synthase (ATP5F1E and ATP5ME) (Fig. 3K and sFig. 12C). In addition, the difference in FGFR3, known as a suppressor of GFAP expression along with astrocyte reactivity [107, 108], as well as the alteration in NFKBIA and TAB2, the NF-kb activator, suggest the major signaling pathways may involved (Fig. 3I-K).
Taken together, the two trajectories of astrocytes, although sharing many similarities, underwent specific genetic changes at the end, suggesting the presence of different types of reactivity in astrocytes or different circumstances in the damaged area.
Spatial correlation analysis revealed that astrocyte subtypes were primarily grouped into two categories based on their association with PD (Fig. 3L-M and sFig. 8B): CON-like (AC1 and AC3) and PD-associated (AC4, AC7, AC8, AC2, and AC6) (Fig. 3L, left panel), which are actually aligned with their pseudotime stages (Fig. 3B). The CON-like astrocytes were located in neuron-rich regions (SC7 and SC8) and the PVALB-high SC4 region. PD-associated subtypes were located in the glia-rich regions (SC10, SC11 and SC6), indicating ongoing reactive gliosis (Fig. 3M, left panel and Supplementary Table 3). Interestingly, in regions with PD-associated subtypes, there was a clear scarcity of CON-like astrocytes, indicating a mutually exclusive population distribution and possibly representing different states of the same cells.
Spatial correlation between late-stage subtypes AC6 and AC8 was high in the medial SN but low in the lateral SN, resulting in completely different groupings in hierarchical clustering (Fig. 3L). In this regard, the abundance score of AC6 was higher than AC8 in SC7 and SC8, but lower in SC0 and SC2, with this difference being more pronounced in the lateral SN (Fig. 3M, right panel). In the medial SN, AC4 and AC7, showed high spatial correlation with PD-derived astrocytes (AC8, AC6 and AC2), but in the lateral SN, these showed high correlation with CON-like astrocyte (AC3 and AC1) (Fig. 3M, middle and right panels). These results indicate the spatiotemporally different effects of PD exerted on medial and later SN [109, 110], which may be related to the different proportion and distribution of AC6 and AC8 astrocyte subtypes in these regions.
Identification of two distinct character transitions in PD microglia
5,391 microglia nuclei were sorted by pseudotime trajectory, displaying five microglia components (MCs) composed of one trunk and four symmetrical branches (Fig. 4A-B and Supplementary Table 10_1). Like astrocytes, microglia trajectory was also separated by CON and PD (Fig. 4C). Moreover, according to the hierarchical enrichment GO analysis, the MC1 corresponds to normal microglia, which regulate homeostasis through AMPK signaling and microtubule anchoring (sFig. 13). Three trajectories were identified: MC1-MC3-MC4 (Trajectory 1), MC1-MC3-MC5 (Trajectory 2), and MC1-MC2 (Trajectory 3) (Fig. 4D-E). Trajectory 1 and 2 (Fig. 4C and Supplementary Table 10_2) shared 788 pseudotime DEGs were categorized into seven Gene sets (Fig. 4F, 4f, and all gene lists for each Gene set are shown in Supplementary Table 10_3). Gene set 2 was associated with interferon-gamma response, macrophages, and IL-6/JAK/STAT3 signaling, indicating inflammatory response. The subsequent upregulation of Gene set 1 in MC5 suggested a full-blown stress response, including UPR and stress granules. A small Gene set 5, whose expression is transiently increased before and after MC3, was associated with dynamic changes in the nuclear lumen and the NF-kB signaling pathway. The largest Gene set 6 showed downregulation of autophagy-related genes.
Next, we noted the common features of downward trajectories 2 and 3. They shared 823 pseudotime DEGs divided into eight gene sets (Fig. 4F, 4f’, and Supplementary Table 10_3). Gene sets 3, 4, and 5 were associated with ribosome and related terms, cytoplasmic translation and coronavirus disease. Gene sets 3 and 4 showed a pattern of transiently increased expression in the middle of the pseudotime, followed by Gene set 5 including the additional term, protein secretion. There was a rebound in the genes of Gene set 6 was associated with ubiquitin protein transferase inhibitor activity, negative regulation of intrinsic apoptotic pathways, and mTORC1 signaling. Gene sets 7 and 8 were upregulated at the end of the pseudotime and were associated with microglial activation, maturation, immune function, and protein secretion.
These results suggest the distinct microglia trajectories associated with PD progression, highlighting the critical role of the MC3 state in determining the fate of microglia and the involvement of inflammatory responses, stress responses, and dysregulated cellular processes in disease-associated trajectories.
We then utilized the q-value of pseudotime DEGs to identify key genes for two phenotypes (Fig. 4G and Supplementary Table 10_2): one for the “Macrophage-like” phenotype (shared in Trajectory 2 and 3) and the other for the “PD-associated” phenotype (shared in Trajectory 1 and 2) (as labeled in Fig. 4E). The former phenotype was enriched with genes involved in ribosomal biogenesis or protein translation, such as RPS28, RPL18A, RPS18, RPS19, and RPLP1, consistent with our GO results (Fig. 4G, right panel). EEF1A1 (microglial polarization) [111] and macrophage markers TYROBP [112] were included. APOE [113], TREM2 [114–116] and DOCK8 [117], which have been implicated in phagocytic function, microglial activation, and relevance to neurodegenerative diseases, were also included. The later phenotype was enriched with genes linked to metal ion transport or ferroptosis such as TMEM163, SLC11A1, and HAMP, which showed high q-values [118–122] (Fig. 4G, left panel). Other top genes shared in all trajectory include CSGALNACT1 (an enzyme for the synthesis of chondroitin sulfate proteoglycan and major component of the glial scar) [123], NIBAN1 (stress response) [124], GAB1 (M2 polarization) [125, 126], and CD163 (macrophage marker) [127]. It also has BACH1 (defense against oxidative stress, mitochondrial dysfunction, and neuroinflammation) [128, 129], and HSPA5 (known as GPR78 or BIP, a regulator of UPR and ER stress) [130, 131].
To examine the spatial distribution of microglia, we examined the pairwise correlation between the abundances of five microglial components in the two Visium slides (Fig. 4H, left panel). Unlike neurons and astrocytes, the five microglial components were highly correlated indicating a similar distribution, which was more pronounced in lateral SN (sFig. 8C). MC1 and MC4 associated with Trajectory 1 were spatially correlated with oligodendrocytes (Fig. 4H, right panel). Consistent with this, MC1 and MC4 were enriched in the oligodendrocyte-rich regions SC11 and SC0 (Fig. 1J, 4I). In contrast, SC10 showed enrichment of the "macrophage-like" MC2, MC3 and MC5 rather than MC1 and MC4 (Fig. 4I). In addition, given that MC2 (Trajectory 3) had PD-independent features in snRNA-seq (sFig. 14C), MC5 was relatively enriched in SC9 (perivascular), SC2 and SC3 compared to MC2 (Fig. 1I). SC9 also showed high MC3 abundance.SC4 and SC5 had significantly fewer microglia.
We sought to confirm these analysis results in postmortem brain tissues using RNAscope (Fig. 4J and Supplementary Table 5). CSF1R is a well-known marker gene predominantly expressed in CNS microglia (Fig. 1F). In the snRNA-seq analysis, CSF1R showed significantly decreased expression specifically in the end of the Trajectory 2, which was associated with the “PD-associated” microglial phenotype (Supplementary Table 10_2). CD163, known as a marker of macrophages and previously associated with the anti-inflammatory M2 type, showed increased expression in MC2, MC3 and MC5, associated with the "macrophage-like" phenotype. However, CD163 relevance was negligible in MC4, which was the end product of Trajectory 1 (“PD-associated” phenotype) (sFig. 15A and Supplementary Table 10_2). In CON tissue, approximately 17% of microglia co-expressed CSF1R and CD163 with perivascular localization (Fig. 4J, top panel). In preclinical PD, the density of CSF1R-positive cells increased, with an increase in CSF1R/CD163 double-positive cells (around 36%) (Fig. 4J, middle panel), indicating early and rapid contribution of CD163 (potentially MC3) (sFig. 15 and Supplementary Table 10_1). In advanced PD, the number of cells positive for CD163, but lacking CSF1R was markedly increased, consistent with the snRNA-seq results.
Identification of four distinct Gene Regulatory Network (GRN) modules in PD glial cells
We constructed the GRN between transcription factors (TFs) and DEGs across the identified trajectories/phenotypes in astrocytes and microglia using TENET [53] and Leiden clustering [132] (Fig. 5A, 5B, Supplementary Table 11 ). Four distinct GRN modules were identified (Fig. 5A, C): Module 1 was a large module comprising all three microglial trajectories, enriched for immune system, cytokine secretion, and vesicle transport processes. Module 2 was a small module, centered on HSPA5, linked to UPR and ribosome biology changes in PD microglia. Module 3, associated with astrocyte trajectory 2, was enriched for mitosis, cell adhesion, potentially tau phosphorylation. Module 4 was associated with astrocyte trajectory 1, enriched for cilia-related function, primarily regulated by the TFs, MSRB3 or EAF2.
While the microglial trajectories clustered into one major module (module 1), distinct TFs were found to preferentially regulate specific trajectories within this module. For example, RREB1 and PPARD in Microglial Trajectory 1, POU2F2 and XRN2 in Microglial Trajectory 2, and ZNF710 and NFATC2 in Microglial Trajectory 3. Considering the prominent role of HSPA5 in module 2 (blue background circles, Fig. 5A) and its association with the UPR [133, 134], the diverse translation-related terms in GO suggest that changes in ribosome biology may be a major response of PD microglia to ER stress, presumably detrimental considering the PD progression in Microglia Trajectory 2 (Fig. 5C) [135, 136].
The astrocyte trajectories formed independent modules (Module 3, 4), suggesting distinct transcriptional regulation (Fig. 5A, C). In Module 3 (Astrocyte trajectory 2, grouped in the yellow circle in Fig. 5A), phosphorylation-related terms are likely the result of the contribution of STOX1, which has been implicated in tau protein phosphorylation (Fig. 5C) [137]. On the other hand, Module 4 (Astrocyte Trajectory 1, grouped in the green circle in Fig. 5A), primarily regulated by the transcription factors MSRB3 or EAF2, shows terms specific to cilia-related functions, consistent with the enriched cilia-related genes in AC8 (Fig. 5C, 3K).
Taken together, network analysis integrating TF-gene interactions helped further refine the distinct transcriptional programs and biological roles of the identified glial trajectories/phenotypes in PD pathogenesis.