To investigate the progression mechanism of PD and perform WGCNA, we first collected bulk RNA-seq data (GSE150646) from an SNCA overexpression rat model, including 20 samples (Hentrich et al., 2020). Quality control of single-cell samples is based on the number of genes detected, the number of total RNA molecules detected, and the percentage of mitochondria or ribosomes from each sample. Since they come from different experiments, we adopt different data filtering methods. Expressly, for GSE140231, this study first excludes the data of two cases of cerebral amyloid angiopathy and it only includes four cortex single cell data (Agarwal et al., 2020), and the data filtering threshold is 10,000 > nCounts_RNA > 1,000, nFeature >200, and the percentage of mitochondrial genes <10%. For GSE157783, all samples were included in the study and the data filtering threshold is nCounts_RNA > 1,000, nFeature >200, and the percentage of ribosomes genes <0.05% (Smajić et al., 2022). Moreover, we used a bulk RNA-seq dataset (GSE205450) of one human brain region for validation (Irmady et al., 2023).

This is an observational study. The Ethics Committee of Qingdao Agricultural University has confirmed that no ethical approval is required.

After obtaining high-quality single-cell data through quality control, the data were processed using the Seurat (R software package v4.4.0) (Hao et al., 2021). The two datasets were then merged, and dimension reduction was performed using the canonical correlation analysis algorithm. The dimensionality reduction and clustering results are visualized by the uniform manifold approximate projection (UMAP) method at a resolution of 0.3 and a dimension of 10 via plot1cell (R software package v0.0.0.9000) or Seurat (R software package v4.4.0) (Hao et al., 2021; Wu et al., 2022). Next, cluster annotation using classic marker genes, Oligodendrocytes [MOBP, MOG (Montague et al., 2006; Juryńczyk et al., 2019)], OPCs [VCAN (van Bruggen et al., 2017)], Astrocyte [AQP4, GFAP (Ikeshima-Kataoka, 2016; Jurga et al., 2021)], Neuronal [GAD1, GAD2 (Kodama et al., 2012)], Microglia [CD74 (Hwang et al., 2017)], Endothelial [EGFL7, CLDN5 (Jang et al., 2011)] and Ependymal [FOXJ1 (Shah et al., 2018)].

For bulk RNA-seq data, the DESeq2 (R software package v1.42.1) was used for the DEGs detection, and the input was the gene counts matrix (Love et al., 2014). Only genes satisfying padj < 0.05 & |log2FoldChange| > 0.5 were considered DEGs. For scRNA-seq data, FindMarkers() function of the Seurat was used for DEGs detection (Hao et al., 2021), and the threshold was set to logfc.threshold >0.1, min.pct > 0.1 and pvalue < 0.05. Unless otherwise specified, the default algorithms and parameters of the software were used.

For functional exploration of the gene sets, including DEGs, hub genes and candidate genes, functional enrichment analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), was performed. The clusterProfiler (R software package v4.10.1) and metascape (v3.5.20240101¹) were used for gene functional enrichment analysis (Wu et al., 2021; Zhou et al., 2019). The significance level of padj < 0.05 was considered as the cut-off threshold.

For bulk RNA-seq data, the WGCNA (R software package v1.72–5) was used for WGCNA, and the input was the gene counts matrix. First, missing values and outlier samples were checked (Langfelder and Horvath, 2008). The hclust() function was used to check outlier samples, and samples with obvious outliers were excluded from the analysis. Next, the pickSoftThreshold() function was used to calculate and pick a suitable soft threshold. The blockwiseModules() function was used to construct overexpression networks and module partitioning in one step, and the parameters we set are power = sft$powerEstimate, maxBlockSize = 6,000, TOMType = “unsigned,” minModuleSize = 300, mergeCutHeight = 0.3 and deepSplit = 2. Subsequently, the labeledHeatmap() function was used to visualize the relationship between modules and phenotypes. Note that since the genes in the gray modules did not participate in the clustering of any module, they were not included in the subsequent analysis. The module membership (MM) and gene significance (GS) algorithm was used to mine hub genes of WGCNA in bulk RNA-seq, which are often closely related to traits (Zhang et al., 2022; Zhang et al., 2023b). The cut-off threshold of hub genes was MM > 0.6 and GS > 0.4.

For scRNA-seq data, the hdWGCNA (R software package v0.3.03) was used for WGCNA, and the input was the integrated Seurat object (Morabito et al., 2023). First, we filtered genes. We only included genes expressed in more than 5% of cells as subsequent genes. Merging multiple cell expression patterns (KNN algorithm) into one metacell can avoid the sparsity of single-cell data. Next, three cell types, Oligodendrocytes, Neuronal, and OPCs, were selected for WGCNA analysis. The TestSoftPowers() function was used for appropriate soft threshold screening and visualization. The ConstructNetwork() function was used for scale-free network construction and module division. The GetHubGenes() function obtained hub genes in different modules.

To investigate the transcriptional dynamics of genes in oligodendrocytes and OPCs during PD progression, we performed pseudo-time trajectories analysis by monocle (R software package v2.24.0), which constructs cell lineage development based on the changes in gene expression levels of different cell subsets over time (Trapnell et al., 2014). The reduceDimension() function is used to determine the trajectory, and the orderCells() function is used to sort cells. Since the algorithm trajectory does not conform to the actual biological law, we manually specified root_state = 2 in this step. The BEAM statistical analysis model was used to calculate cell fate trajectories before and after key fate nodes, using unsupervised clustering genes as markers.

To characterize the differences in signal transduction pathways between normal physiological conditions and PD, the CellChat (R software package v1.6.1) algorithm was used to analyze intercellular communication at the single-cell level (Jin et al., 2021). The Cell Communication Database uses CellChatDB.human as a reference. The global cell–cell communication network between normal and PD group was quantitatively and comparatively analyzed using the compareInteractions() function. The netVisual_bubble() function was used to visualize the differences in cellular communication between OPCs-neuronal and oligodendrocytes-neuronal.

Proteins interact with each other to form complex interaction networks to regulate various aspects of life processes, such as gene expression regulation, cell cycle regulation, etc. (Tomkins and Manzoni, 2021). The candidate genes were used as inputs for PPI construction using the String online database.² The default parameters were run, and we removed node proteins with no interactions during visualization.

ROC curve analysis was used to evaluate the diagnostic value of candidate genes for PD and to obtain the final biomarker. The pROC (R software package v1.18.5) was used to calculate ROC curve and visualization. Data from separate datasets (GSE205450) were used for input.

PD is a neurodegenerative disease that is closely related to age (Coukos and Krainc, 2024). However, the mechanism of PD progression has not been fully elucidated. Here, we first analyzed the data from the SNCA overexpression rat model and showed that SNCA overexpression could disrupt rat frontocortical gene transcription at an early stage. A total of 2,072 DEGs were obtained, including 1,433 up-regulated genes and 639 down-regulated genes (Supplementary Figure S1A). GO enrichment analysis showed that items related to learning and cognition were significantly enriched, such as “learning or memory,” “cognition,” “learning” (Supplementary Figure S1B). KEGG pathway results suggested that “Neuroactive ligand-receptor interaction” was the most enriched pathway (Supplementary Figure S1C). Interestingly, in the late SNCA overexpression stage (named PD_old group), only 211 DEGs were detected (Supplementary Figure S1D). The GO terms results showed that they were related to “glycosylceramide metabolic process,” “ensheathment of neurons,” and so on (Supplementary Figure S1E), and no KEGG pathway was enriched.

To gain insight into the pathogenesis of PD, especially in its early stages, we used the WGCNA approach. The hierarchical clustering tree showed that one sample had significant outlier performance, so we removed it (Figure 1A; Supplementary Figure S2A). We constructed a scale-free co-expression network for the remaining samples. First, we performed soft threshold selection. The results showed that it performed well and met the goal of building a scale-free network when β = 4 (Figure 1B). Fifteen modules were obtained, among which the turquoise module had the largest number of genes, the cyan module had the least number of genes, and the grey module contained no meaningless genes (Figure 1C). We next attempted to explore the relationship between these modules and PD progression, and the results indicated that the brown module was obviously positively correlated with PD progression (r = 0.83, p = 1e-05) and clearly negatively correlated with the normal phenotype (r = −0.83, p = 1e-05) (Figure 1D). Furthermore, we constructed a clustering matrix of PD phenotypes and modules, and the results showed that the brown module was closely associated with PD (Figure 1E).

After identifying the key modules, we explored the functions of the genes within the brown modules. The GO annotation results showed that the genes in the module were involved in “oligodendrocyte differentiation,” “oligodendrocyte development,” “glial cell development,” “regulation of gliogenesis,” etc. Moreover, “oxidoreduction-driven active transmembrane transporter activity” was also enriched (Figure 2A). This result suggests that abnormalities in oligodendrocytes may drive the progression of PD. Next, the MM&GS algorithm identified the hub gene in the brown module (Figure 2B), and functional annotation of hub genes was also performed. The GO annotation results showed that the hub genes were related to “oligodendrocyte differentiation,” “oligodendrocyte development,” “glial cell development,” “ATP synthesis coupled electron transport,” etc. (Figure 2C). The KEGG pathway results showed that the hub genes were involved in “oxidative phosphorylation” (Figure 2D). This result suggests that mitochondrial dysfunction is also involved in the progression of PD. Subsequently, we compared the relationship between hub genes and DEGs. The results showed that hub genes had 162 overlapping DEGs in the two groups (Figure 2E). The GO annotation of 162 genes indicated that they were associated with “oligodendrocyte differentiation,” “mitochondrial respiratory chain complex I,” etc. (Figure 2F). The above results strongly suggested that oligodendrocyte abnormalities are significant for the progression of PD.

Considering that the above analysis was performed in a rat model, to verify our hypothesis that oligodendrocyte abnormalities are involved in PD progression, we reanalyzed two human brain tissue scRNA-seq data sets (see detail in Supplementary Table S1). The UMAP map showed that 47,562 cells, including 28,621 cells in the normal group and 18,941 cells in the PD group, were obtained by dividing into 17 clusters, and a total of 7 cell types were annotated (Figure 3A; Supplementary Figure S3A). The dot plot showed the expression characteristics of related cell type marker genes in different groups, Oligodendrocytes [MOBP, MOG (Montague et al., 2006; Juryńczyk et al., 2019)], OPCs [VCAN (van Bruggen et al., 2017)], Astrocyte [AQP4, GFAP (Ikeshima-Kataoka, 2016; Jurga et al., 2021)], Neuronal [GAD1, GAD2 (Kodama et al., 2012)], Microglia [CD74 (Hwang et al., 2017)], Endothelial [EGFL7, CLDN5 (Jang et al., 2011)] and Ependymal [FOXJ1 (Shah et al., 2018)] (Figure 3B). Moreover, we used UMAP to display the maps of different groups and different samples (Supplementary Figures S3B,C).

Taking the WGCNA results into account, OPCs were also included in the analysis because they are precursors of oligodendrocytes. First, we extracted oligodendrocyte populations from the integrated data and re-clustered and visualized them. The results showed that there were significant differences in OPCs between the PD and Normal groups (Figure 3C). 1,435 DEGs were detected in OPCs (Supplementary Table S2). This result seems to suggest that PD progression was driven by OPCs. Interestingly, the top GO term of DEGs in OPCs was “regulation of nervous system development” (Figure 3E). The KEGG pathway results showed that “Chemical carcinogenesis - ROS” was most enriched (Supplementary Figure S3D). Next, the same analysis was performed on oligodendrocytes, and the UMAP showed the differences between PD and Normal groups (Figure 3D). Further, we performed differential analysis with a total of 851 DEGs (Supplementary Table S3). The GO annotation showed that the 851 DEGs were participated in “protein folding,” etc. (Figure 3F). The KEGG pathway results suggested that they were related to neurodegenerative diseases, such as “Parkinson disease,” “Huntington disease,” “Alzheimer disease” (Supplementary Figure S3E).

To dissect the fate decisions of OPCs/oligodendrocytes throughout the study period, they were sorted according to their gene expression patterns in a pseudo-temporal trajectory. We defined the stage in which OPCs were as the initial stage of the trajectory, and then a total of three cell stages were obtained, which were classified into two different cell fates (Figures 4A,B). Subsequently, we attempted to explore the transcriptional regulatory program of OPCs differentiation into oligodendrocytes, and the results suggested that PI3K/AKT/mTOR signaling may play a key role (Supplementary Figures S4A,B). Next, we observed the proportion of cells in different groups at three different cell stages, and we found that the proportion of PD cells was higher in state_3. We speculated that the stage where state_3 was located may be driven by the development of PD (Figure 4C). Next, we used the BAEM algorithm to try to parse the driving factors that lead to different cell fate decisions. Interestingly, we found that the emergence of cell fate_1 was driven by two distinct gene expression patterns. The SNCA gene was highly expressed in cell fate 1 and lowly expressed in cell fate_2, which is consistent with the fact that high SNCA expression promotes PD progression (Hentrich et al., 2020; Figures 4A,D). We divided the genes driving cell fate decisions into four gene sets. C2 was highly expressed in cell fate_1(PD), and functional enrichment analysis showed that “regulation of autophagy” and “Parkinson disease” were significantly enriched (Figure 4E). C3 was highly expressed in cell fate_2 (Normal), and functional enrichment analysis showed that “regulation of TOR signaling” and “cell cycle” were significantly enriched (Figure 4F).

To investigate the dynamics of cell–cell communications between oligodendrocytes/OPCs and neurons, we performed cell–cell communication analysis using CellChat (R software package v1.6.1). The results found that the number and intensity of cell communication in the PD group were higher than those in the normal group (Supplementary Figure S5A). Interestingly, cellular communication between oligodendrocytes/OPCs and neuronal cells was enhanced in PD group (Supplementary Figure S5B). Next, the signals of cell communication in different samples were aligned, and we found two signaling networks related to nervous system development have strong information flows, including NRG, NRGX signaling network (Supplementary Figure S5C). Interestingly, our cell communication networks observed that neurons communicated strongly with OPCs in both NRG and NRXN signaling (Figures 5A,B). Specifically, NRG3-ERBB4 was only present between OPC-neurons (Figure 5C). This result suggests that PD may be driven by OPCs because the communication between OPCs-neurons is closer than that between oligodendrocytes (Maldonado and Angulo, 2015; Xiao and Czopka, 2023; Boulanger and Messier, 2017; Orduz et al., 2015).

To construct the gene co-expression network in OPCs/Oligodendrocytes/Neuronal and identify hub genes involved in PD, hdWGCNA analysis was performed. The hdWGCNA results suggested that when the soft threshold β = 4, scale-free network construction can be performed (Figure 6A). The hierarchical clustering tree showed the relationship between modules and genes (Figure 6B). A total of 8 valuable modules were identified (Figure 6C). The UMAP plot was used to display the OPCs/Oligodendrocytes/Neuronal gene co-expression network and highlighted the hub genes of each module involved in nervous system development (Figure 6D). The dot plot showed that the turquoise module is highly expressed in PD group, and the co-expression network of the top 25 hub genes was displayed (Figures 6E,F). Gene function annotation results showed that the hub genes of the turquoise module were involved in “negative regulation of neurogenesis,” “JNK kinase binding,” and so on (Figure 6G).

We compared the relationship between hub genes and DEGs from scRNA-seq data, and the results showed that hub genes had 45 overlapping DEGs in OPCs groups. In contrast, it had no overlap with oligodendrocytes, which suggesting that PD progression may be more sensitive to OPCs. In other words, studying OPCs may be more critical (Figure 7A). Only four terms were enriched, including “regulation of RNA splicing,” “cell junction assembly,” “negative regulation of neuron projection development” and “protein polyubiquitination” (Figure 7B). Then, the PPI network was built, and showed PPP2R2B, DNM3 and ATXN1 may play a significant role (Figure 7C). Importantly, we sought to identify markers that drive PD and performed ROC binary analysis model in an independent dataset (Supplementary Table S1). The results showed that AGPAT4, DNM3, PPP1R12B, PPP2R2B and LINC00486 levels could distinguish PD patients from healthy controls. The area under the ROC curve was followed by 0.764, 0.745, 0.725, 0.721, and 0.711, which suggesting these five gene may serve as the potential biomarkers for predicting PD (Figure 7D).