Revealing a novel contributing landscape of ferroptosis-related genes in Parkinson's disease.

Transcriptomics studies have yielded great insights into disease processes by detecting differentially expressed genes (DEGs). In this study, due to the high heritability of Parkinson's disease (PD), we performed bioinformatics analyses on nine transcriptomic datasets regarding substantia nigra from Gene Expression Omnibus database, including seven microarray datasets and two next-generation sequencing datasets. As a result, between age-matched PD patients and normal control, we identified 630 DEGs, of which 22 hub DEGs involved in PD or ferroptosis were found to be associated with each other at the transcriptional level and protein-protein interaction network, suggesting their high correlations among these hub genes. Moreover, 16 DEGs were singled out due to their comparable AUC (>0.6) in random forest classifiers, including seven PD-related genes (MAP4K4, LRP10, UCHL1, PAM, RIT2, SNCA, GCH1) and nine ferroptosis-related genes (GCH1, DDIT4, RGS4, MAPK9, CAV1, RELA, DUSP1, ATP6V1G2, ATF4 and ISCU). Furthermore, to probe the potential of those hub genes in predicting the PD progression and survival, we constructed a Cox model featured by an eight-gene signature, including four PD-related genes (SNCA, UCHL1, LRP10, and GCH1) and four ferroptosis-related genes (DDIT4, RGS4, RELA, and CAV1), and validated it successful in an independent dataset, indicating that it would be an effective tool for clinical research to predict PD progression. In conclusion, ferroptosis-related DEGs identified in this study were closely correlated with the known PD-related genes, revealing the involvement of ferroptosis in the development of PD. This study presented the potential of several ferroptosis-related genes as novel clinical biomarkers for PD.

Introduction

Parkinson’s disease (PD) is the second most common age-depended neurodegenerative disease [1], [2]. The prevalence of PD increases with age, and reaches up to 5 % in individuals older than 85 years [3]. With an aging society around the world, the number of PD patients is expected to increase more than 1.5 times by 2030, creating a huge social burden [4]. Clinically, PD patients are often characterized by several motor syndromes traits and/or non-motor symptoms traits. That is, resting tremor, rigidity, bradykinesia, and postural instability are the typical motor syndromes, while dementia, depression and sleep disorders belong to non-motor symptoms [1]. Pathological hallmarks of PD include the loss or pathological changes of dopamine neurons in the substantia nigra pars compacta associating with the nigrostriatal pathway [5].

In the recent 20 years, with the development of high-throughput sequencing technologies, many rare variants and candidate genes have been reported to contribute familial and sporadic cases of PD [6], [7]. Although several environmental factors have also been identified to contribute to the development of PD [8], more than 100 loci associated with PD risk were screened out at the genome-wide significance level by genome-wide association studies (GWAS), which could be the potential risk factors and biomarkers for PD [8], [9]. Our team have replicated more than nine PD-related genes in Chinese population [10]. Furthermore, we have also reported gene expression of PD associated genes and age at onset features for the first time [11].

As a newly discovered mode of cell death different from apoptosis and other known cell programmed death pathways, ferroptosis is characterized by iron-dependent lipid peroxidation accumulation [12], [13]. Simultaneously, several PD features have been well established as known key features and/or triggers in the ferroptosis cell death pathway, such as iron overload [14], reduced GSH levels [15], elevated lipid peroxidation [16], DJ-1 depletion [17], XCT downregulation [18], and CoQ10 reduction [19]. This strongly implicated that ferroptosis plays a role in the neurodegeneration observed in PD [20]. Indeed, PARK7 and PLA2G6, two known PD-deleterious genes linked to autosomal-recessive early onset PD, were reported to be associated with the regulation of ferroptosis, in which PARK7 was identified to act as a ferroptosis inhibitor by preserving the transsulfuration pathway [17] while PLA2G6 could be strongly associated with increased lipid peroxidation levels [21]. That indicated that there were some ferroptosis-related genes in the occurrence and development of PD. However, it remains unclear which ferroptosis-related genes are involved in the pathogenesis of PD.

Genome-wide expression profiling, or transcriptomics is a novel powerful approach for generation of new research method for the pathogenesis of PD [22]. High throughput transcriptomics study is a great tool for detecting DEGs from two or more samples [23]. This will be very helpful to improve the detection of the pathogenic genes and identify functional pathways among DEGs to reveal biological themes contributing to the disease process [22]. Due to high heritability in PD, transcriptomics study has the potential to reveal significant insights into disease processes.

The integration of the data from different studies would increase the statistical power when prioritizing disease related genes [24]. In order to more accurately study the important role of gene transcription in the pathogenesis of PD, we would integrate previously published transcriptomic data based on the substantia nigra of the brain, and re-analyze the differentially expressed genes (DEGs), with matching the basic situation of patients and controls in this study. This study provides a novel comprehensive landscape of ferroptosis-related genes for clinical diagnosis and therapy for PD.

Materials and methods

Data source and processing

We searched the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) with keywords “Parkinson” and “substantia nigra”. A total of nine transcriptomic datasets were obtained and presented in Table 1 and Table S1, including microarray and next-generation sequencing (NGS). Herein, the raw “CEL” files of those microarray datasets, i.e. GSE7621, GSE20141, GSE49036, GSE20164, GSE20163, GSE20292 and GSE34865, were downloaded, normalized, and log2 transform with R package affy (version 1.68.0). Then, we merged those datasets sequenced by microarray to remove batch effect with R package sva (version 3.38.0), including 66 Parkinson’s disease (PD) samples and 114 normal control (NC) samples. Similarly, the two NGS datasets, i.e. GSE166024 and GSE114517, were also merged to remove batch effect with R package sva, including 31 PD samples and 12 NC samples. Those PD-related genes were collected, including PD-deleterious genes from a recognized systematic review [25] and PD-risk genes from iPDGC Locus Browser [26] integrating with three GWAS researches.

Table 1

The datasets used in this study.

GEO datasets	Platform	Method	Tissue	PD samples	NC samples
GSE7621	GPL570	microarray	substantia nigra	16	9
GSE20141	GPL570	microarray	substantia nigra	10	8
GSE49036	GPL570	microarray	substantia nigra	15	8
GSE20164	GPL96	microarray	substantia nigra	6	5
GSE20163	GPL96	microarray	substantia nigra	8	9
GSE20292	GPL96	microarray	substantia nigra	11	18
GSE34865	GPL517	microarray	substantia nigra	0	57
GSE166024	GPL20301	NGS	substantia nigra	14	0
GSE114517	GPL18573	NGS	substantia nigra	17	12

Identification of DEGs and enrichment analysis

The differentially expressed genes (DEGs) between PD samples and NC samples were identified by using R package limma (version 3.46.0). We set thresholds at the P value of less than 0.05 and the absolute log2 fold change (log2FC) of greater than 0.4. Subsequently, the down-regulated DEGs and up-regulated DEGs were separately applied for GO enrichment and KEGG enrichment by using R package clusterProfiler (version 3.18.1).

Correlation analysis of DEGs and PPI network analysis

The Pearson’s correlations between DEGs were calculated and visualized by R package Corrplot (version 0.92). In addition, the STRING database (https://string-db.org/) was downloaded and used to filter the protein–protein interactions (PPI) among DEGs. Then, the Cytoscape software (version 3.6.1) were used to visualize the PPI network. Therein, nodes represented genes, and their sizes denoted node degree in the network. Edges stood for the interactions between genes, and their sizes reflected the degree of relationships.

Construction of classification model

The R package randomForest (version 4.6–14) was applied to construct classification models, in which we set the parameters of ntree and mtry at 500 and 4, respectively. Here, we randomly selected 80 % samples as training set, and the rest 20 % samples as test set. The receiver operating characteristic (ROC) curves were plotted by R package pROC (version 1.18.0), and the area under curves (AUC) were used to evaluate the predictive power of classifiers.

Construction of prognostic model

Those important DEGs mentioned above were considered as features to perform multivariable Cox proportional hazards regression model by R package survival (version 3.2–13). Therein, Kaplan–Meier survival curve was implemented by R package survminer (version 0.4.9) and log-rank test was used to determine the significance. The risk score was fitted as follows:in which the i denotes the i-th feature, and the Feature and Coef represent its expressed value and coefficient in the fitted Cox model, respectively. And, the Coef is equal to the natural logarithm of hazard ratio (HR).

In dataset, the median of risk score was applied to stratify samples into high- and low-risk groups. And, between high- and low-risk groups, univariable Cox proportional hazards regression analysis was applied to determine the HR, 95 % confidence interval (CI), and P value.

Results

Identification of differentially expressed genes in PD

As shown in Table 1, seven datasets sequenced by microarray were obtained from GEO and merged with the removal of batch effect, including 66 PD samples and 114 normal control (NC) samples. When ignoring the effect of samples’ age, all samples derived from PD and NC groups were applied to perform differential analysis. However, only 182 DEGs were identified (Table S2A). As the age of most samples used in this study is not available, to accurately identify DEGs, the age-matched PD and NC samples were singled out, including 21 PD samples and 21 NC samples. Meanwhile, we identified 630 DEGs, including 161 up-regulated genes and 469 down-regulated genes (Fig. 1A, Table S2B). Interestingly, in addition to the 182 DEGs, more genes were identified in the age-matched PD and NC samples (Fig. S1A), which seems to be consistent with the fact that PD is an age-related disease.

Fig. 1

Identification and analysis of differentially expressed genes in Parkinson’s disease (PD). (A) Volcano plot shows the differentially expressed genes between age-matched PD samples and normal controls. (B) The significant KEGG pathways enriched by down-regulated genes in PD. (C) The significant GO biology process enriched by down-regulated genes in PD.

To further explore the involved biology pathways in DEGs identified in age-matched PD and NC samples, those down-regulated DEGs were separately applied for KEGG enrichment analysis and GO enrichment analysis, respectively (Table S3A, B). In the significant KEGG pathways enriched by those down-regulated genes, in addition to PD-related pathways, i.e. Parkinson disease and synaptic vesicle cycle, we also found several important pathways, such as oxidative phosphorylation and chemical carcinogenesis - ROS (Fig. 1B). Furthermore, similar biology process (BP) related to PD were also observed in GO enrichment analysis, such as axonogenesis, synapse assembly, synaptic vesicle endocytosis, and synaptic vesicle exocytosis (Fig. 1C). Simultaneously, the enrichment of cell component (CC) and molecular function (MF) was also found to involve synaptic vesicle, presynaptic membrane, dopaminergic synapse, and syntaxin-1 binding (Fig. 1C). Notably, several key BP were enriched as well, such as mitochondrial ATP synthesis coupled electron transport, transferrin transport, and iron ion transport. And, several key MF were enriched as well, such as ATPase-coupled transmembrane transporter activity and ferrous iron binding (Fig. 1C). Thus, we speculated that the ferroptosis played a role in patients with PD, leading to the occurrence and development of PD.

Screening of ferroptosis-related DEGs in PD

To explore the relationship between ferroptosis and PD, we firstly wonder whether there were ferroptosis-related genes in the 630 DEGs of PD samples. We obtained 259 ferroptosis-related genes from FerrDb [27], and 107 PD-related genes that we previously reviewed from literatures [11] (i.e. PD-deleterious genes or PD-risk genes). Subsequently, as shown in Fig. 2A, we found that nine DEGs belonged to PD-related genes (including three PD-deleterious genes [25], [28] (SNCA, UCHL1, LRP10) and six PD-risk genes [29] (GCH1, SH3GL2, MAP4K4, PAM, INPP5F, RIT2), and 14 DEGs were considered as ferroptosis-related genes (i.e. GCH1, GPX4, RGS4, RELA, MAPK9, ATP6V1G2, CAV1, ISCU, GOT1, ATF4, DUSP1, SLC7A5, DDIT4, VEGFA), suggesting that the 22 unique DEGs may act as the link between ferroptosis and PD (Fig. 2B, Fig. S1B). Especially, GCH1 was considered as a PD-related gene, and also belonged to a ferroptosis-related gene. The PD-related genes SNCA and PAM, as well as the ferroptosis-related genes MAPK9, ATP6V1G2, and GOT1, were clustered in BP of response to metal ion. Moreover, genes SNCA and ISCU were enriched in MF of ferrous iron binding, and genes ISCU and ATP6V1G2 were enriched in BP of iron ion transport. It suggested that PD-related genes were highly correlated with ferroptosis.

Fig. 2

Analysis of 22 hub differentially expressed genes in PD. (A) Venn diagram presents the genes intersection. (B) The abundance of 22 hub differentially expressed genes in age-matched PD samples and normal controls. (C) Pearson’s correlation at the transcriptomic level among the 22 hub differentially expressed genes. (D) Protein-protein interaction network among the 22 hub differentially expressed genes. Nodes in cyan, in pink, and in red represent the ferroptosis-related genes, PD-risk genes, and PD-deleterious genes, respectively. Edges denote the interaction, and its size reflect the intensity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Furthermore, the 22 hub DEGs mentioned above were further applied to investigate their Pearson’s correlations and protein–protein interactions (PPI), respectively. Primarily, as shown in Fig. 2C, we found that some DEGs were positively correlated with each other, in which MAPK9, GPX4, ATP6V1G2, RGS4, GOT1, and ISCU were considered as ferroptosis-related genes, while SNCA, PAM, INPP5F, UCHL1, SH3GL2, and RIT2 belonged to PD-related genes. These significant positive correlations indicated that they were associated at the transcriptional level.

Moreover, an interacted functional network containing 18 of the 22 hub DEGs were developed based on the STRING database (Fig. 2D). In particular, we observed that those ferroptosis-related genes were interacted with each other, and were interacted with PD-related genes, such as SNCA, UCHL1, GCH1, SH3GL2, RIT2, LRP10, and MAPK4. These seemed to indicate the relationship between ferroptosis and PD at the PPI network.

Construction of classifiers based on the key DEGs in PD

To probe the potential as characteristics to PD, these 22 identified hub DEGs were used to construct the random forest classifiers, respectively, and their AUCs were evaluated. Meanwhile, 16 classifiers, including MAP4K4, LRP10, UCHL1, PAM, RIT2, SNCA, GCH1, DDIT4, RGS4, MAPK9, CAV1, RELA, DUSP1, ATP6V1G2, ATF4 and ISCU, were found to achieve the AUC of greater than 0.6, indicating that their featured genes had a promising potential to characterize PD (Fig. 3). In particular, the classifiers featured by PD-related genes MAP4K4, LRP10 and UCHL1 showed the AUCs of greater than 0.7 (Fig. 3A-C), and the classifiers featured by ferroptosis-related genes DDIT4, RGS4, MAPK9 and RELA also showed the AUCs of greater than 0.7 (Fig. 3H-K). That indicated the potential of these genes as classification biomarkers for PD, and specifically suggested these ferroptosis-related genes may act as novel PD-related genes.

Fig. 3

Random forest classifiers with the AUC of greater than 0.6 constructed by the hub genes. (A-G) ROC curves of the PD-related genes, including MAP4K4, LRP10, UCHL1, PAM, RIT2, SNCA, and GCH1. (H—P) ROC curves of the ferroptosis-related genes, including DDIT4, RGS4, MAPK9, RELA, CAV1, DUSP1, ATP6V1G2, ATF4, and ISCU.

Survival analysis based on the key DEGs in PD

Since the substantia nigra samples used in this study were all sampled within a few hours after the death of the patients, thus the key DEGs mentioned above were further used to explore their potential for survival analysis in PD. Here, we selected those PD samples with age information from those microarray samples, and then regarded the age of patients as survival time to perform multivariable Cox proportional hazards regression analysis. As shown in Fig. 4A, when eight of these features were combined, the Cox model exhibited an excellent concordance index (0.79) with significance of less than 0.05. Notably, we observed that the ferroptosis-related gene DDIT4 presented a significant hazard ratio (HR) of greater than 1, indicating that the greater the expression of DDIT4, and the faster the PD pathogenic process. Indeed, DDIT4 was defined as a mark of ferroptosis in ferrDB, and Tian et al reported silencing DDIT4 partially reversed the erastin-induced ferroptosis in LUAD cell line, validating the role of DDIT4 as ferroptosis driver in vitro [30]. In addition, we found that the ferroptosis-related genes CAV1, RGS4 and RELA all presented significant HR of less than 1, indicating that the greater the expressed value, the lower the risk of poor survival. In ferrDB, CAV1 [31] and RELA [32] was defined as a suppressor of ferroptosis, while RGS4 was defined as marker. Overexpression of CAV1 was indeed found to inhibit the process of ferroptosis in head and neck squamous cell carcinoma and autoimmune hepatitis [31], [33].

Fig. 4

Survival analysis based on an eight-gene signature. (A) Forest plot of the Cox model. (B) Kaplan-Meier survival curve comparing the high- and low-risk groups in the training dataset. (C) Kaplan-Meier survival curve comparing the high- and low-risk groups in the independent dataset.

After fitting the risk score of each PD sample with the eight characteristic genes, we observed that the median of risk score was able to significantly stratify these samples into high- and low-risk groups (Fig. 4B, HR = 2.72, 95 % CI: 1.58–4.67, P = 9e-05), and the high-risk group exhibited a poor overall survival (age).

To validate the significance of the above prognostic model, an independent test dataset merged by two NGS datasets was used. Similarly, the risk score of each PD samples was fitted by the eight characteristic genes and their coefficients in the training model. Expectedly, the median risk score also remarkably divided samples into high- and low-risk groups, and high-risk group exhibited a poor overall survival (Fig. 4C, HR = 2.27, 95 % CI: 1.03–5.01, P = 0.03858), validating that the eight genes had good prognostic ability.

Discussion

Efficient high-throughput sequencing and bioinformatics analysis have contributed to our understanding of the molecular mechanisms of disease occurrence and development, which is necessary to identify novel potential clinical biomarkers. In this study, transcriptome datasets from substantia nigra were integrated and used to explore the (DEGs characterized to PD by bioinformatics analysis. Interestingly, in addition to some known PD-related signaling pathways, some DEGs were also found to enrich in ferroptosis-related biology pathways, such as iron ion transport, oxidative phosphorylation and chemical carcinogenesis - ROS. Particularly, the PD-related genes SNCA and PAM, as well as the ferroptosis-related genes MAPK9, ATP6V1G2, and GOT1 were enriched in the same biology pathway, i.e. response to metal ion. That suggested that ferroptosis may contribute to the occurrence and development of PD.

To further investigate the occurrence of ferroptosis in PD patients, 22 hub DEGs that have been reported to relate with PD or involve in ferroptosis were screened out for subsequent analysis. Therein, significant correlations at the transcriptional level and protein–protein interaction between PD-related genes and ferroptosis-related genes were observed, revealing the 22 unique DEGs may act as the link between ferroptosis and PD. It is worth pointing out that GCH1 appeared in both PD-related genes and ferroptosis-related genes. We have identified GCH1 as a PD-risk gene in Chinese PD patients [34], [35], and it was also reported to suppress ferroptosis by selectively preventing depletion of phospholipids with two polyunsaturated fatty acyl tails [36]. Besides, GPX4, a known key regulator of ferroptosis [37], was demonstrated to be up-regulated in neurons of substantia nigra and associated with dystrophic axons in striatum of Parkinson's brain [38]. Furthermore, 16 of the 22 hub genes were used to respectively construct classifiers to characterize PD, in which seven PD-related DEGs (MAP4K4, LRP10, UCHL1, PAM, RIT2, SNCA and GCH1) stood out. As we know, SNCA, one of the first identified PD pathogenic genes, encodes the α-synuclein protein that constitutes the major protein component of Lewy bodies [39]. However, the abnormal aggregation of α-synuclein plays a key role in the formation of the LBs and other α-synuclein pathological aggregates and is regarded as a critical step in the molecular pathogenesis of PD [40], [41]. Besides, UCHL1 is another gene that was first reported to cause PD two decades ago, while the pathogenicity of LRP10 is firstly identified in 2018 [42]. MAP4K4, PAM, RIT2 and GCH1 were report to be risk genes working to be associated with neuroprotection [43] and neurological function [44]. Besides, nine ferroptosis-related DEGs were found to classify PD samples from normal controls with comparable AUCs, suggesting the high association between ferroptosis and PD. Especially, four random forest classifiers featured by DDIT4, RGS4, MAPK9 or RELA all presented the AUC of greater than 0.7, suggesting that these ferroptosis-related DEGs may be used as diagnostic targets for PD. Moreover, a significant prognostic model characterized by an eight-gene signature was constructed to stratify PD samples into high- and low-risk groups, in which the high-risk group presented a poorer prognosis with successful validation in an independent data. In the eight-gene signature, in addition to four known PD-related DEGs, four ferroptosis-related genes were also included, i.e. DDIT4, RGS4, RELA and CAV1, suggesting that these ferroptosis-related DEGs may be used as therapeutic targets for PD. Indeed, DDIT4 was reported to be an inhibitor of mTOR pathway [45], while RGS4 is related with system x(c)(-)-mediated glutamatergic dysfunction [46]. Both the two signal pathways are closely related to the pathogenesis of PD. Above all, these hub genes are highly associated with PD. It follows that these ferroptosis-related DEGs were the potential PD risk genes.

In addition, there are three limitations in this present study. Firstly, the integrated dataset derived from several sequencing platforms was used to identified the DEGs, although we corrected the bias by removing the batch effect to maintain the reliability of research results as much as possible. Secondly, after strictly matching basic information such as age, the number of samples used in the study was small, which may reduce the statistical effectiveness of the research results. However, we believed that the improved statistical power by strict matching was more obvious and more significant. Thirdly, the hub DEGs in this study need verified the pathogenic mechanism with cell lines or animal experiments.

Conclusion

Taken together, based on the transcriptomics data of substantia nigra from PD patients, we performed functional analysis, such as enrichment and correlation analysis, classification model and prognostic model construction, to reveal the potential contribution of ferroptosis related genes to PD. Notably, we firstly reported that the risk score fitted by the eight genes was used to reflect the development of PD and predict the survival of patients. Meanwhile, several ferroptosis-related genes were, for the first time, uncovered to differentially express in PD, which could be used as novel potential clinical targets for PD.

CRediT authorship contribution statement

Xingxing Jian: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing. Guihu Zhao: Investigation, Methodology, Software, Writing – original draft. He Chen: Investigation. Yanhui Wang: Investigation. Jinchen Li: Investigation, Supervision, Validation, Writing – original draft. Lu Xie: Supervision. Bin Li: Conceptualization, Investigation, Software, Supervision, Writing – original draft, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.