Novel Susceptibility Genes and Molecular Mechanisms Identified in Relation to Mitochondrial Dysfunctions in Parkinson’s Disease

Abstract

Background: Mitochondrial dysfunction is one of the major contributors in the pathogenesis of Parkinson’s Disease (PD). However, the mechanisms and genes involved in mitochondrial dysfunction in PD have still not been fully uncovered.

Methods: The PD associated datasets were accessed via the Gene Expression Omnibus (GEO) database, and we extracted 1870 Mitochondrial Dysfunction Related Genes (MDRGs) via the Genecards database. First, Differentially Expressed Genes (DEGs) between PD and control samples were screened out by differential expression analysis, and samples (PD and control) were considered as traits for Weighted Gene Co-Expression Network Analysis (WGCNA) to yield the key modules. Then, we took the intersection of DEGs and genes from key modules to acquire a set of intersected genes, which were then intersected with MDRGs to obtain candidate genes. Next, the Protein-Protein Interaction (PPI) analysis was implemented on candidate genes for identifying the hub genes, and they were subjected to single Gene Set Enrichment Analysis (GSEA). Finally, we constructed the drug-gene, Transcription Factor (TF)-mRNA-miRNA networks, and LASSO-logistic regression model.

Results: Between the PD and control samples, 945 DEGs were filtered out. Subsequently, MEblue and MEyellow modules all had a strong correlation with traits, therefore they were considered as the key modules. We also yielded 444 intersected genes and 31 candidate genes, respectively. In addition, 7 hub genes: DVL2, DNMT1, ABL1, RAF1, NOTCH1, RELA, and PDGFRB were identified via the protein-protein interaction analysis. These genes were found to be related to ‘Parkinson's disease’, ‘oxidative phosphorylation’ and other functional pathways. Likewise, drugs such as resveratrol and bortezomib were predicted according to the hub genes, and the TF-mRNA-miRNA network consisting of relationship pairs, such as SREBF1-ABL1, and ‘hsa-mir-23b-3p’-NOTCH1, was constructed. Ultimately, we found that the LASSO-Logistic regression model had an excellent ability to distinguish PD samples from control samples.

Conclusion: Through bioinformatic analyses, we identified novel mitochondrial dysfunction related hub genes including DVL2, DNMT1, ABL1, RAF1, NOTCH1, RELA, and PDGFRB in PD, and explored their mechanisms of action in the tricarboxylic acid cycle and neuroactive ligand receptor interaction signaling pathways, providing new insights for studying the pathogenesis of PD.

Keywords

Parkinson’s disease; Mitochondrial dysfunction; GEO; Bioinformatics

Introduction

Parkinson’s Disease (PD) is one of the most common neurodegenerative diseases affecting 2–3% of the population over the age of 65 [1]. Patients with PD are characterized by cardinal motor symptoms including bradykinesia, resting tremor, rigidity and posture and gait disturbances; as well as non-motor symptoms including but not limited to hyposmia, constipation, autonomic dysfunctions, impaired cognitive function, and rapid eye movement sleep behavior disorder [2]. These symptoms increase in severity as the disease gradually progresses and no current drugs have been found to stop or slow down the disease course. Although enormous efforts have been made to reveal the disease mechanisms, the pathogenesis of this devastating disease is still not fully understood. Aging, environment, and genetic influences play key roles in the pathogenesis of PD, with each of these factors varying in significance, interfering with key molecular pathways including mitochondrial homeostasis, endo-lysosomal pathway, neuroinflammation and oxidative stress, and causing gradual neuronal loss throughout the disease course [3].

Mitochondrial dysfunction has been found to be tightly associated with PD since 1989 when mitochondrial complex I deficiency was found in the post-mortem substantia nigra of PD patients [4]. As the powerhouse of cells, mitochondria not only produce adenosine triphosphate through oxidative phosphorylation, but also participate in other essential cellular activities including calcium homeostasis, programmed cell death, mitophagy, oxidative stress and neuroinflammation [5–7]. Mutations in several PD causative genes including SNCA, LRRK2, PARK2, PINK1, or DJ-1 have been found to result in mitochondrial dysregulation and contribute to the disease. For example, the role of PARK2/PINK1 pathway in mitophagy which clears damaged mitochondria has been extensively studied [8,9]. Moreover, several other mitochondrial function-related genes have been implicated in PD risk [10]. However, further identification of novel genes involved in PD mitochondrial dysfunction will shed light on the discovery of new disease mechanisms and the development of innovative disease-modifying therapeutics.

Therefore, in this study, we accessed genetic datasets of PD patients and control cohorts from the public Gene Expression Omnibus (GEO) database, extracted 1870 Mitochondrial Dysfunction Related Genes (MDRGs) via the Genecards database. We identified seven mitochondrial dysfunction related hub genes that included DVL2, DNMT1, ABL1, RAF1, NOTCH1, RELA and PDGFRB through a series of bioinformatic methodologies including Differentially Expressed Genes (DEGs), Weighted Gene Co-Expression Network Analysis (WGCNA) and Protein-Protein Interaction (PPI) analysis. We have also constructed the drug-gene, Transcription Factor (TF)-mRNA-miRNA networks, and LASSO-Logistic regression model which confirmed the mechanisms of these hub genes in PD pathogenesis.

Materials and Methods

Acquisition of the data

The GSE26927 training set (chip data) (12 Parkinson’s Disease (PD) samples and 8 control samples) and GSE20163 (8 PD samples and 9 control samples), GSE7621 (16 PD samples and 9 control samples) validation sets were gained via Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds). Following this, 1,870 Mitochondrial Dysfunction Related Genes (MDRGs) were extracted via the Genecards database (https://www.genecards.org/) (Reference score >4).

Differential expression analysis and Weighted Gene Co-Expression Network Analysis (WGCNA)

Differentially Expressed Genes (DEGs) between PD and control samples from the GSE26927 training set were screened out by limma package (version 3.50.1) [11] setting adj. P<0.05. We considered samples (PD and control (from the GSE26927 training set)) as traits for WGCNA. First, samples were clustered to remove outliers, and a soft threshold was determined for the data. Gene modules were gained by constructing a co-expression matrix. Further, correlations of gene modules with traits were analyzed, and modules that had strong correlation with traits were selected as key modules.

Functional enrichment analysis of candidate genes

The DEGs and genes from key modules were intersected to yield a set of intersected genes, which were then intersected with MDRGs to acquire a set of candidate genes. Next, to explore the functional pathways involved by candidate genes, we implemented Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses by cluster Profiler package (version 4.0.5) [12] setting adjusted (adj.) P<0.05. Moreover, the location and expression pattern of candidate genes on chromosomes were analyzed utilizing the OmicCircos package (version 1.32.0) [13].

Screening for hub genes

First, we constructed a protein-protein interaction (PPI) network of candidate genes via the STRING database (http://string-db.org), and Cytohubba plug-in was applied to calculate the degree of each protein node for identifying the hub genes. After that, the expression of hub genes between PD and control samples was analyzed by setting adj. P<0.05 in the GSE26927 training set. We also analyzed the Pearson correlations among hub genes. The functional similarity analysis of hub genes was carried out utilizing the GOSemSim package (version 2.20.0) [14]. Finally, we implemented Gene Set Enrichment Analysis (GSEA) on hub genes to explore the functional pathways in which they were involved (background gene set: KEGG: c2.cp.kegg.v2022.1.Hs.entrez.gmt).

Construction of the drug-gene, Transcription Factor (TF)-mRNA-miRNA networks, and LASSO-Logistic regression model

We utilized the Drug-Gene Interaction database (DGIdb) (https://www.dgidb.org/) to predict the small molecule drugs of hub genes, and a drug-gene network was constructed based on the predicted results. Afterwards, the miRNAs and TFs of hub genes were respectively predicted via mirNet (https://www.mirnet.ca/) and JASPAR (https://jaspar.genereg.net/) databases, and a TF-mRNA-miRNA regulatory network was constructed based on the above factors. Finally, Least Absolute Shrinkage and Selection Operator (LASSO) analysis was implemented on hub genes via glmnet package (version 4.1-2) [15] to construct the LASSO-Logistic regression model, and confusion matrix heat map was drawn to analyze the ability of LASSO-Logistic regression model distinguishing PD from control samples. Moreover, we also plotted the Receiver Operating Characteristic (ROC) curve in the GSE26927 training set to assess the diagnostic ability of LASSO-Logistic regression model, and the evaluation result was validated by the same method in the GSE7621 and GSE20163 validation sets.

Results

Identification of DEGs and key modules

There were 945 DEGs between PD and control samples (Figure 1A), and the heat map revealed the |log2FC| top 10 up and down- regulated DEGs (Figure 1B). In addition, Figure 1C illustrated that the overall clustering of dataset samples was good, and there was no need to eliminate samples. According to the position of the red line in Figure 1D, the soft threshold was determined to be 7. At this point, the vertical coordinate R2 was approximately 0.85, and the mean value of the adjacency function was gradually approached to 0, indicating that the network was close to the scale-free distribution and showed a flat trend. Next, we filtered out 21 modules by constructing the co-expression matrix. Ultimately, we found that MEblue and MEyellow modules all had strong correlations with traits, therefore they were considered as the key modules (contained 3,127 genes in total) (Figure 1E).

Figure 1. The volcano plot of DEGs, red: upregulated genes; blue: downregulated genes; grey: genes not differentially expressed (A). The heatmap of DEGs, red: upregulated genes; blue: downregulated genes (B). Sample clustering for detecting outliers (C). Scale independence and mean connectivity analysis (D). Module-trait relationship (E).

Functional pathways involved in candidate genes

We finally yielded 444 intersected genes and 31 candidate genes, respectively (Figures 2A and 2B). Moreover, candidate genes were involved in GO entries such as ‘glial cell differentiation’, ‘mitochondrial inner membrane’, ‘apical part of cell’, ‘phosphatidylinositol 3−kinase binding’, and ‘positive regulation of cell development’ (Figure 2C). Meanwhile, the KEGG pathways in which they engaged including but not limited to ‘neurotrophin signaling pathway’, ‘notch signaling pathway’, ‘C-type lectin receptor signaling pathway’, ‘sphingolipid signaling pathway’, and ‘phospholipase D signaling pathway’ (Figure 3A). Finally, the heat map illustrated the location and expression pattern of candidate genes on the chromosomes (Figure 3B).

Figure 2. DEGs related to PD (A). DEGs associated with mitochondria (B). GO analysis of the modular genes, adjusted P value <0.05 was considered significant (C).

Figure 3. KEGG enrichment analysis of hub genes(A). Chromosomal location and expression pattern of differentially expressed genes in mitochondrial dysfunctions of PD. From inside to outside, first circle: Heat map of healthy controls; second circle: Heat map of Parkinson’s patients; third circle: Chromosomal map; fourth circle: The name and chromosomal location of each gene(B).

Acquisition of hub genes

The PPI network contained NOTCH1-SOX10, DVL2-VHL, PDGFRB-DNM2 and other reciprocal relationship pairs (Figure 4A). Additionally, we subsequently identified 7 hub genes (DVL2, DNMT1, ABL1, RAF1, NOTCH1, RELA, and PDGFRB) according to the degree of PPI network analysis (Figure 4B), and they were all significantly differently expressed between PD and control samples (Figure 4C). Thereafter, we found strong positive correlations among most of the hub genes, and NOTCH1 had the highest positive correlation with ABL1 (Figure 5A). The box plot revealed that PDGFRB had the highest similarity to other hub genes (Figure 5B). Nevertheless, hub genes jointly participated in ‘Parkinson's disease’ and ‘oxidative phosphorylation’ of these two KEGG pathways. Whereas ABL1, DNMT1, DVL2, NOTCH1, and RAF1 were all enriched to ‘cardiac muscle contraction’ and ‘neuroactive ligand receptor interaction’, while RELA and PDGFRB were jointly engaged in ‘citrate cycle TCA cycle’.

Figure 4. Protein-protein interaction network of differentially expressed genes in mitochondrial dysfunctions of PD (A). The identification of hub genes (B). Differential expression of hub genes (C).

Figure 5. Correlation analysis of the expression of hub genes (A). Functional similarities among hub genes, the higher the score of a gene is, the more functional relevance of this gene is to other genes (B).

The drug-gene, TF-mRNA-miRNA networks, and LASSO-Logistic regression model

The drug-gene network included resveratrol-RELA, bortezomib-NOTCH1, decitabine-DNMT1 and other relationship pairs (Figure 6A). Whereupon the TF-mRNA-miRNA network consisted of relationship pairs such as SREBF1-ABL1, PPARG-RELA, and ‘hsa-mir-23b-3p’-NOTCH1 (Figure 6B). Through LASSO logistic regression analysis, the optimal λ value was determined to be 0.01101454 (Figures 7A and 7B). Besides, the confusion matrix heat map revealed that the LASSO-Logistic regression model had an excellent ability to distinguish PD samples from control samples (Figure 7C). Concurrently, the Area Under the Curve (AUC) value of the LASSO-Logistic regression model was 1, indicating that it could diagnose PD patients with a high degree of accuracy (Figure 7D). The evaluation results in the GSE7621 and GSE20163 validation sets further demonstrated the predictive power of the LASSO-Logistic regression model.

Figure 6. Prediction of small molecule drugs through the drug-gene work analysis (A). Transcription factor-mRNA-miRNA regulation network analysis, pink oval: hub genes; yellow rectangle: miRNAs targeting at least two hub genes; green oval: transcription factors (B).

Figure 7. LASSO logistic regression analysis, Log(λ) curve (A); Lasso confidence interval (B). Confusion matrix heat map (C). Area under the Curve (AUC) (D).

Discussion

The More than 90 independent genetic risk signals have been found as common genetic components of PD and around 20 genes are identified as in monogenic PD forms [16]. Among these genes, several have been reported to participate in maintaining mitochondrial homeostasis. SNCA, the first gene described to cause PD in 1997, encodes the alpha-synuclein protein which has a mitochondrial targeting sequence at its N-terminal, and regulates mitochondrial morphology and complex I activity [17]. Accumulation of alpha-synuclein due to missense mutations or duplication/triplication of the SNCA gene was reported to increase oxidative stress and cytochrome C release, leading to dopaminergic neuron death [18]. The PARK2/PINK1 signaling pathway is the main regulator in mitochondrial quality control and maintains mitochondrial dynamics through mitophagy. Loss-of-function mutations in the PARK2/PINK1 pathway contribute to mitochondrial dysregulation and selective dopaminergic neuron loss in familial PD patients. However, in sporadic PD cases where there are no specific PARK2/PINK1 pathway defects, other genes may be involved. In this study, we have discovered 7 novel hub genes: DVL2, DNMT1, ABL1, RAF1, NOTCH1, RELA, and PDGFRB that are closely related to mitochondrial dysfunction. The identification of these new hub genes may help shed light on potential disease mechanisms, which could give rise to the development of novel disease-modifying treatments.

The release of cytochrome C, production of reactive oxygen species, activation of AMP-activated protein kinase, and release of mitochondrial DNA are the four main mechanisms by which mitochondria communicate with the rest of the cellular organelle to control cell fate and function [19-22]. The release of Tricarboxylic Acid (TCA) cycle metabolites has been recently reported to act as a fifth mitochondrial molecular mechanism. These metabolites enact functions through regulating DNA methylation, chromatin modifications, and post-translational protein modifications [19]. The significantly increased expression levels of RELA and PDGFRB in PD patients compared to controls, as demonstrated in this study for the first time, linked the release of TCA cycle metabolites to PD mitochondrial dysfunctions. Although the exact mechanisms of how RELA or PDGFRB regulate TCA cycle metabolites warrant further investigations, our drug-gene network analysis has shown an interaction between RELA and resveratrol, a small molecule that could potentially be repurposed as a disease-modifying treatment for PD. In addition, Peroxisome Proliferator-Activated Receptor (PPARG) that was reported to regulate mitochondrial biogenesis and neuronal survival elsewhere [23] was correlated with RELA expression, which is consistent with previous study [24]. These findings may imply potential indirect strategies to manipulate RELA gene and protein expression; in addition, RELA inhibitors or nucleic acid therapies targeting RELA may offer a more direct suppression approach. However, this will need to be validated further as a treatment strategy targeting PD mitochondrial dysfunctions.

Emerging evidence has shown a role for the neuroactive ligand-receptor interaction signaling pathway in modulating gene expression and regulating neuronal functions [25]. Furthermore, the function of this pathway was demonstrated to be affected by alpha-synuclein-induced microRNA dysregulation in vivo in an PD drosophila model [26]. Although little evidence has linked the neuroactive ligand-receptor interaction pathway with mitochondrial functions, rotenone, an inhibitor of mitochondrial complex I [27], has been shown to inhibit the transcriptional activity of key genes involved in this pathway [28]. Our KEGG analysis showed that five hub genes including ABL1, DNMT1, DVL2, RAF1, and NOTCH1 were enriched in the neuroactive ligand-receptor interaction pathway. Small molecules including bortezomib and decitabine, and microRNA including hsa-mir-23b-3p were predicted by the drug-gene and TF-mRNA-miRNA network analysis to interact with genes involved in the neuroactive ligand-receptor interaction pathway. These factors should be verified as potential therapeutic strategies to correct dysregulation in this pathway that has been related to impaired cognitive functions [29]. A common non-motor manifestation of PD throughout the disease course that can severely compromise quality-of-life [30].

Although enriched in the neuroactive ligand-receptor interaction pathway, NOTCH1 is required for the synaptic plasticity that contributes to memory functions and learning [31]; while, the Notch signaling pathway controls signal integration and the fate of glial and neural stem cells during development [32,33]. Overexpression of Notch signaling post-development, may however contribute to neuronal cell death in response to noxious stimuli, as is the case for hypoxic-ischemic injuries [34]. As is the case with increased expression level of NOTCH1 was found in human brain microvascular endothelial cells and iPSC-derived neuronal cells obtained from AD patients [35], our data reports a significant increase of NOTCH1 in postmortem brain tissues in PD patients compared to controls. Conversely, blocking the Notch pathway has been shown to provide significant protection against neuronal death [36]. Thus, suppressing NOTCH1 activities using the predicted hsa-mir-23b-3p, or other approaches including antisense oligonucleotides, may provide a strategy for reducing the gradual neuronal loss which occurs in PD.

Conclusion

In summary, this study applied a series of bioinformatic package tools and identified previously unreported hub genes that are involved in novel disease mechanisms of PD. The data from this research demonstrates for the first time the involvement of TCA cycle metabolites regulated by RELA and PDGFRB in PD mitochondrial dysfunctions. Although the number of samples analyzed in the training set and validation sets is limited and further validation required, the outcome of the study may serve to redirect the exploration of PD pathogenesis and drug development once follow-up investigations on this exploration of mitochondrial dysfunctions are completed.

Ethics Approval and Consent to Participate

Not applicable.

Consent for Publication

Not applicable.

Availability of data and Materials

All data generated from this study are available upon request.

Competing Interests

The authors declared no conflict of interests.

Funding

This research was funded by the Bryant Stokes Neurological Research Fund, National Natural Science Foundation of China (82101505), and the College Key Research Project in Henan Province (22A180030).

Author Constributions

DL and TW contributed to the conceptualization and design of the study. DL and FZ collected and analyzed the data and led the writing of the original draft. DL, SM, and CM edited, reviewed and revised the manuscript. All authors had full access to the data and approved the submission of the manuscript.

Acknowledgement

We would like to thank Professor Frank Mastaglia (Perron Institute for Neurological and Translational Science) for his helpful review of this manuscript. DL is very grateful to the generous salary support from the philanthropic Giumelli Family Foundation.