Journal Logo

Cellular, Molecular and Developmental Neuroscience

Bioinformatics analysis and identification of genes and molecular pathways involved in Parkinson’s disease in patients with mutations in the glucocerebrosidase gene

Xu, Dan-Dan*,; Li, Guo-Qian*,; Wu, Zhi-Sheng; Liu, Xiao-Qiang; Yang, Xiao-Xia; Wang, Jie-Hua

Author Information
doi: 10.1097/WNR.0000000000001685



Parkinson’s disease is a relatively common degenerative disease of the central nervous system (CNS) [1]. Its incidence accounts for 2–3% of people over 65 years of age, causing a huge economic burden on families and society [2,3]. At present, it is generally believed that the onset of Parkinson’s disease may be caused by various factors such as aging, inheritance, genetic mutations and environmental toxins [1,4,5]. Among genetic factors, the glucocerebrosidase (GBA) gene is most commonly associated with Parkinson’s disease in the population. GBA is located in the 1q21 region; it is 7600 bp long and contains 11 exons. It encodes a protein called glucocerebrosidase, which catalyzes the breakdown of glucosylceramide into ceramide and glucose [6,7]. In recent years, scholars worldwide have made remarkable achievements in the study of the correlation between GBA gene mutations and Parkinson’s disease [8–10]. Many studies have confirmed that mutations in this gene, even for heterozygotes, increase the risk of developing the disease [11,12]. It has been determined that GBA gene mutations in patients with Parkinson’s disease are more common than in the general population [13,14]. An Israeli team first discovered a mutation in the GBA gene in German Jewish Parkinson patients [15]. GBA harmful mutations associated both with the onset of Gaucher disease and Parkinson’s disease in a heterozygous state include mutations N370S and L444P [16,17]. GBA mutations cause glucocerebrosidase deficiencies which may lead to abnormalities in lysosomal autophagy pathways and cause the accumulation of α-synuclein [18,19].

The aim of this study was to identify potential crucial genes and pathways associated with GBA gene mutations in patients with Parkinson’s disease to provide new insights into the role of GBA in the mechanisms underlying the development of this disease.


Dataset collection

The gene expression profiles GSE53424 and GSE99142 based on GPL10558 (Illumina HumanHT-12 V4.0 expression bead chip) were acquired from the Gene Expression Ominibus (GEO) database ( The GSE53424 dataset contained three Parkinson’s disease patients with heterozygous glucocerebrosidase mutations (GBA N370S) and two idiopathic Parkinson’s disease patients (IDI-PD). The GSE99142 dataset contained two Parkinson’s disease patients with heterozygous glucocerebrosidase mutations (GBA N370S) and three idiopathic Parkinson’s disease patients.

Identification of differentially expressed genes

We separated the data into IDI-PD 1 group and GBA-PD 1 group [Parkinson’s disease patients with heterozygous glucocerebrosidase mutations (GBA N370S) mutations] for further analysis. The GSE53424 and GSE99142 datasets were merged and normalized using the ‘sva’ R package. We identified differentially expressed genes (DEGs) between IDI-PD 1 and GBA-PD 1 groups in Parkinson’s disease using the ‘limma’ package in R. Values with P < 0.05 and |log2Fold change (logFC) | >2 were considered statistically significant.

Gene ontology enrichment and Kyoto Encyclopedia of Genes and Genomes pathway analysis

We used Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 ( to do the gene ontology enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to investigate the DEGs at the molecular and functional levels. Gene ontology enrichment analysis was used to annotate the cellular components, the biologic processes and the molecular functions of DEGs. KEGG pathway analysis was used to determine the route of gene cluster and related functions. When P < 0.05, we believe that statistics are significantly different.

Protein–protein interaction network analysis

Using the online database, STRING (STRING is part of the ELIXIR infrastructure: it is one of ELIXIR’s;, the protein–protein interactions (PPIs) of the DEGs were selected with score (median confidence) >0.4. The PPI network between IDI-PD 1 and GBA-PD 1 groups in Parkinson’s disease was then visualized using Cytoscape Software (Version 3.4.0; The results obtained by PPI network analysis were further module analyzed using Cytoscape Software. During the analysis, molecular complex detection (MCODE) was used to identify the most important module of the network map. The criteria for analysis were a degree cutoff of 2, MCODE scores >5, max depth of 100, k-score of 2 and node score cutoff of 0.2. Hub genes were excavated by setting the degrees. We used OmicShare to analyze hub gene clustering.

Statistical analysis

Statistical significance of the difference between two groups was determined using Student’s t-tests. The results were expressed as the means ± SD unless otherwise noted. The expression of hub genes for GBA mutations was compared for correlation analysis using Pearson rho (ρ) tests. The ability of hub genes to predict GBA mutations was determined using receiver operating characteristics (ROC) curves. All data were statistically analyzed using SPSS software, version 23.0 (IBM Corp, Armonk, New York, USA). Values with P < 0.05 were considered significantly different.


Identification of differentially expressed genes

Comparisons between the transcription profile datasets of the IDI-PD 1 and GBA-PD1 groups were obtained from National Center for Biotechnology Information GEO databases. DEGs from GSE53424 and GSE99142 series performed on GPL10558 platforms were screened, with P < 0.05 and |logFC| ≥2. For the analysis, five samples were included in the IDI-PD 1 group and five in the GBA-PD1 group. A total of 283 DEGs (130 upregulated and 153 downregulated) were identified after the comparison. Heatmap visualization of these DEGs was combined with a hierarchical cluster analysis; results suggested that both groups could be clearly distinguished (Fig. 1a). A DEG expression volcano plot of the upregulated and downregulated genes of the two expression groups is shown in Fig. 1b.

Fig. 1
Fig. 1:
Differential expressed genes (DEGs) analysis. (a) Heatmap of 283 DEGs. The diagram presents the result of a two-way hierarchical clustering of all the DEGs and samples. (b) Volcano map of DEGs. Magenta dots represent genes with |logFC | >2 and P value <0.05.

Enrichment analyses of differentially expressed genes

The gene ontology and KEGG databases were used for enrichment of biologic themes and pathways involving identified DEGs and used to feed the DAVID for further analysis. Gene ontology enrichment analysis predicted the functional roles of target host genes on the basis of three aspects, including biologic processes, cellular components and molecular functions. Enrichment analyses of the upregulated DEGs (Fig. 2a, Table 1) showed that the top three terms in biologic processes included chemical synaptic transmission, regulation of ion transmembrane transport, and defense response to bacterium. The cellular component groups were integral component of plasma membrane, voltage-gated calcium channel complex and myosin complex. The molecular functions of GBA mutations in Parkinson’s disease mainly focused on ion channel activity and leukotriene receptor activity. Pathway enrichment analysis showed that these DEGs were mainly involved in calcium signaling pathway and Rap1 signaling pathway (Fig. 2b, Table 2). Enrichment analyses of the downregulated DEGs (Fig. 2a) showed that the top three terms in biologic processes included multicellular organism development, G-protein coupled receptor signaling pathway and immune response. The cellular component groups were extracellular space, extracellular region and cell surface. The molecular functions of GBA mutations in Parkinson’s disease mainly focused on calcium ion binding, chemokine activity (Table 1). Pathway enrichment analysis showed that these DEGs were mainly involved in the cytokine–cytokine receptor interaction (Fig. 2b, Table 2).

Table 1 - The biologic process, cellular component and molecular functions in enriched analysis of differentially expressed genes between AAA and control (top three according to gene count)
Term Name Count P value
BP GO:0007268 Chemical synaptic transmission 5 3.4E−02
GO:0034765 Regulation of ion transmembrane transport 4 1.9E−02
GO:0042742 Defense response to bacterium 4 3.7E−02
CC GO:0005887 Integral component of plasma membrane 14 3.2E−02
GO:0005891 Voltage-gated calcium channel complex 3 1.0E−03
GO:0016459 Myosin complex 3 2.8E−02
MF GO:0005216 Ion channel activity 3 2.0E−02
GO:0004974 Leukotriene receptor activity 2 2.0E−02
BP GO:0007275 Multicellular organism development 11 6.1E−04
GO:0007186 G-protein coupled receptor signaling pathway 11 2.7E−02
GO:0006955 Immune response 9 2.3E−03
CC GO:0005615 Extracellular space 18 1.3E−03
GO:0005576 Extracellular region 16 3.4E−02
GO:0009986 Cell surface 8 3.2E−02
MF GO:0005509 Calcium ion binding 11 7.5E−03
GO:0008009 Chemokine activity 3 3.2E−02
BP, biologic processes; CC, cellular components;GO, gene ontology; MF, molecular functions.

Table 2 - The Kyoto Encyclopedia of Genes and Genomes pathway analysis of differentially expressed genes
Term Count P value Genes
UP Calcium signaling pathway 6 4.9E−03 P2RX5, CYSLTR2, LTB4R2, ADCY2, CACNA1E, CACNA1G
Rap1 signaling pathway 5 4.1E−02 FGF6, EFNA2, ITGB3, ADCY2, FYB
Down Cytokine–cytokine receptor interaction 7 1.5E−03 CCL25, IL1A, CCL24, CXCL11, TNFRSF18, LTB, IL17B

Fig. 2
Fig. 2:
Gene ontology and significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses of the differentially expressed genes (DEGs). (a) biological process, cellular component and molecular function between IDI-PD 1 and GBA-PD 1 groups.(b) Significantly enriched KEGG pathways of DEGs. BP, biologic processes; CC, cellular components;MF, molecular functions.

Protein–protein interaction network and modular analyses

A PPI network was constructed using the STRING database to study the relationship between DEGs and to identify hub genes (Fig. 3a,c). Additionally, we used the Cytoscape plug in MCODE to perform a module analysis which showed the top modules. The significant modules of the upregulated and downregulated DEGs based on module analysis of the PPI network contained 6 and 4 nodes, respectively (Fig. 3b.d). The node degrees in the top hub genes of the DEGs are shown in Fig. 3e. Functional enrichment analysis of the hub genes in the modules was mainly related to calcium signaling pathway, Rap1 signaling pathway and cytokine–cytokine receptor interaction.

Fig. 3
Fig. 3:
The protein-protein interaction (PPI) network constructed via the STRING (STRING is part of the ELIXIR infrastructure: it is one of ELIXIR’s) database for the differentially expressed genes (DEGs) with the colors of nodes representing the degree of gene interaction. (a) PPI network of the downregulated DEGs and the upregulated DEGs (c) identified from GSE53424 and GSE99142. The subnetworks were identified by Cytoscape molecular complex detection (MCODE) plug-in (b: upregulated, d: downregulated). (e) Genes from the two subnetworks with the degree of interaction in the PPI network.

Hub genes have clinical significance for Glucocerebrosidase mutations in Parkinson’s disease

The expression levels (after normalization) of the hub genes between Parkinson’s disease with GBA mutations and Parkinson’s disease with GBA no-mutations are shown in Supplementary Table, Supplemental digital content 1, The expression of the hub genes corticotropin-releasing hormone (CRH) was significantly downregulated, whereas that of melatonin receptor 1B (MTNR1B) was significantly upregulated in GBA mutations in Parkinson’s disease (Fig. 4a). The expression of CRH correlated negatively with MTNR1B (Fig. 4b). The area under the ROC curve indicated that the expression levels of CRH and MTNR1B have potential diagnostic value and might serve as biomarkers of GBA mutations in Parkinson’s disease (Fig. 4c).

Fig. 4
Fig. 4:
The relationship and their diagnostic value between DEGs in Parkinson’s disease with the GBA gene mutations. (a) The expression of corticotropin-releasing hormone (CRH) was upregulated and melatonin receptor 1B (MTNR1B) were downregulated in the Parkinson’s disease with the GBA gene mutations. (b) The expression levels of CRH and MTNR1B showed the negative correlation in the Parkinson’s disease with the GBA gene mutations. (c) Receiver operating characteristic analysis of the sensitivity and specificity of the predictive value of the CRH expression model and MTNR1B expression model.


In recent years, many studies have confirmed that GBA gene mutations are a very important and common risk factor for Parkinson’s disease [7]. Mutations in GBA cause reduced or null levels of the enzyme, resulting in an accumulation of glucocerebroside in macrophages (that cannot be hydrolyzed), that is ultimately stored in lysosomes. Consequently, abnormal lysosomal autophagy pathways occur, also causing an accumulation of α-synuclein [20]. Studying the underlying mechanism between GBA mutations and Parkinson’s disease will provide us with new insights into the treatment of Parkinson’s disease and might lead us to proposing new treatment strategies for Parkinson’s disease, either to alleviate or prevent the progression of this disease. To date, the molecular mechanism of the increased risk of Parkinson’s disease in GBA mutations carriers has not been fully elucidated. Although there are many studies on the mechanism underlying the effects of GBA mutations in Parkinson’s disease [11,21], several assumptions and models derived from them have great limitations. In our study, we identified the DEGs between IDI-PD1 and GBA-PD1 and then conducted the function annotation and the pathway enrichment analyses of these DEGs to explore for potentially crucial genes, either of pathogenetic or therapeutic relevance, in GBA-PD1. These analyses led to the identification of CRH and MTNR1B genes as key biomarkers that could be of mechanistic relevance for GBA-PD1 pathogenesis and progression.

CRH is synthesized and secreted by the hypothalamus, which can promote the release of adrenocorticotropic hormone from the pituitary gland. In addition, CRH also exists in other parts of the CNS. It was involved in various physiologic activities such as learning, memory, emotional response and neuroprotection. Studies have shown that the anxious behaviors of transgenic mice overexpressing CRH are significantly increased. Anxiety behavior in transgenic mice with upregulated CRH expression was significantly increased. The wild-type mice also developed anxiety signs similar to the transgenic mice with CRH overexpression [22]. Relevant research shows that CRH can increase the excitability of neurons. In addition, it can also participate in the establishment of long-term hippocampal enhancement and cerebellar long-term inhibition [23]. In-vitro studies have shown that CRH can resist a variety of oxidative stress (e.g. β-amyloid peptide) or excitatory amino acid (e.g. glutamate) damage to primary cultured or passaged neurons. In addition, CRH can also inhibit the neuronal apoptosis caused by LY294002 (PI3K inhibitor) [24].

Melatonin, also known as melanin, is an indole hormone mainly secreted by the pineal gland. It is widely distributed in the human body and has important physiologic effects. It can not only maintain the synchronization of biologic environment coordination, regulate biologic circadian rhythm, resist oxidation, antitumor activity, sexual maturation and reproductive immune response, but also play an important role in maintaining the normal nervous system function [22]. Melatonin receptors are widely distributed in the hypothalamus, pituitary, reproductive ovary, uterus, fallopian tube and testis. The melatonin that enters the nucleus binds to the receptor and then form a complex to regulate gene transcription in the target cell. After melatonin binds to specific receptors, it enters the cell, activates the corresponding second signal transduction system, and then passes through the third messenger of the nucleus, such as cAMP response element binding (CREB) protein, affects the expression of neurotrophic factor (BNDF), and then change the function of nerve cells, leading to sleep disorders, biologic rhythm disorders and depression. When melatonin binds to a specific receptor, it enters the cell and activates the corresponding second signal transduction system. Then, it passes through a third messenger in the nucleus, such as CREB protein, to affect the expression of neurotrophic factor (BNDF). In turn, it changes the function of nerve cells, leading to sleep disorders, biologic rhythm disorders and depression. The retina is the starting point of the central circadian rhythm system, and it is rich in dopamine and melatonin. Studies have indicated that visual abnormalities in Parkinson’s disease animal models and patients may be related to the loss of dopamine neurons and dopamine dysfunction in the retina [25,26]. As an important neuromodulatory hormone in the retina, melatonin exerts biologic effects by regulating the circadian rhythm, neuroprotection and retinal physiology and other corresponding receptors [27]. In recent years, some scholars have also believed that there is a mutual inhibition relationship between melatonin and dopamine in the retina. Melatonin may interfere with the synthesis and release of dopamine by the MT2 receptor [28]. Studies have found that the motor symptoms of Parkinson’s disease rats induced by 6-OHDA have a good correlation with the level of melatonin in the eye, suggesting that melatonin may be involved in the development of Parkinson’s disease [29].

We sequenced the GEO datasets GSE53424 and GSE99142 based on GPL10558 (Illumina HumanHT-12 V4.0 expression bead chip) with high consistency and reliability. We then merged and normalized the datasets, constructed a co-expression network, detected gene modules, and identified two hub genes (CRH and MTNR1B) and three signal pathways (calcium signaling pathway, Rap1 signaling pathway and cytokine–cytokine receptor interaction) that differed between IDI-PD 1 and GBA-PD 1 groups. These hub genes and signaling pathways might have crucial biologic functions in the pathogenesis of patients with Parkinson’s disease bearing GBA gene mutations. Further studies are needed to elucidate their potential role as diagnostic, prognostic or therapeutic biomarkers of GBA gene mutations in Parkinson’s disease.


The present study had some limitations. Our results were only obtained through bioinformatics analysis; they were not demonstrated by real-time PCR or tested in the animal experiment. In addition, the DEGs and key genes were identified for Parkinson’s disease patients with GBA gene mutations, a study using bioinformatics analysis is just the first step and there is still a long way until any of these findings can be translated into clinical applications. Even so, the results might provide a new insight into the molecular mechanism and potential biomarkers for Parkinson’s disease patients with GBA gene mutations.


The authors wish to thank Quanzhou First Hospital Affiliated to Fujian Medical University.

This study was funded by Health and Health Scientific Research Talent Training Project of Fujian Province (grant number 2019-1-92) and Science and Technology Plan Project of Quanzhou (grant number 2019N032S).

Conflicts of interest

There are no conflicts of interest.


1. Poewe W, Seppi K, Tanner CM, Halliday GM, Brundin P, Volkmann J, et al. Parkinson disease. Nat Rev Dis Primers. 2017; 3:17013
2. Vasconcellos LF, Pereira JS. Parkinson’s disease dementia: diagnostic criteria and risk factor review. J Clin Exp Neuropsychol. 2015; 37:988–993
3. Bäckström D, Linder J, Jakobson Mo S, Riklund K, Zetterberg H, Blennow K, et al. NfL as a biomarker for neurodegeneration and survival in Parkinson disease. Neurology. 2020; 95:e827–e838
4. Kyle K, Bronstein JM. Treatment of psychosis in Parkinson’s disease and dementia with Lewy bodies: a review. Parkinsonism Relat Disord. 2020; 75:55–62
5. Konovalova J, Gerasymchuk D, Parkkinen I, Chmielarz P, Domanskyi A. Interplay between MicroRNAs and oxidative stress in neurodegenerative diseases. Int J Mol Sci. 2019; 20:E6055
6. Tsuang D, Leverenz JB, Lopez OL, Hamilton RL, Bennett DA, Schneider JA, et al. GBA mutations increase risk for Lewy body disease with and without Alzheimer disease pathology. Neurology. 2012; 79:1944–1950
7. Do J, McKinney C, Sharma P, Sidransky E. Glucocerebrosidase and its relevance to Parkinson disease. Mol Neurodegener. 2019; 14:36
8. Huh YE, Chiang MSR, Locascio JJ, Liao Z, Liu G, Choudhury K, et al. β-Glucocerebrosidase activity in GBA-linked Parkinson disease: the type of mutation matters. Neurology. 2020; 95:e685–e696
9. Balestrino R, Schapira AHV. Glucocerebrosidase and Parkinson disease: molecular, clinical, and therapeutic implications. Neuroscientist. 2018; 24:540–559
10. Magalhaes J, Gegg ME, Migdalska-Richards A, Doherty MK, Whitfield PD, Schapira AH. Autophagic lysosome reformation dysfunction in glucocerebrosidase deficient cells: relevance to Parkinson disease. Hum Mol Genet. 2016; 25:3432–3445
11. Horowitz M, Pasmanik-Chor M, Ron I, Kolodny EH. The enigma of the E326K mutation in acid β-glucocerebrosidase. Mol Genet Metab. 2011; 104:35–38
12. Pankratz N, Beecham GW, DeStefano AL, Dawson TM, Doheny KF, Factor SA, et al.; PD GWAS Consortium. Meta-analysis of Parkinson’s disease: identification of a novel locus, RIT2. Ann Neurol. 2012; 71:370–384
13. Sidransky E, Samaddar T, Tayebi N. Mutations in GBA are associated with familial Parkinson disease susceptibility and age at onset. Neurology. 2009; 73:1424–1425
14. Gegg ME, Schapira AHV. The role of glucocerebrosidase in Parkinson disease pathogenesis. FEBS J. 2018; 285:3591–3603
15. Aharon-Peretz J, Rosenbaum H, Gershoni-Baruch R. Mutations in the glucocerebrosidase gene and Parkinson’s disease in Ashkenazi Jews. N Engl J Med. 2004; 351:1972–1977
16. Beutler E, Gelbart T, Scott CR. Hematologically important mutations: Gaucher disease. Blood Cells Mol Dis. 2005; 35:355–364
17. Lesage S, Anheim M, Condroyer C, Pollak P, Durif F, Dupuits C, et al.; French Parkinson’s Disease Genetics Study Group. Large-scale screening of the Gaucher’s disease-related glucocerebrosidase gene in Europeans with Parkinson’s disease. Hum Mol Genet. 2011; 20:202–210
18. Schöndorf DC, Aureli M, McAllister FE, Hindley CJ, Mayer F, Schmid B, et al. iPSC-derived neurons from GBA1-associated Parkinson’s disease patients show autophagic defects and impaired calcium homeostasis. Nat Commun. 2014; 5:4028
19. Aflaki E, Borger DK, Moaven N, Stubblefield BK, Rogers SA, Patnaik S, et al. A new glucocerebrosidase chaperone reduces α-synuclein and glycolipid levels in iPSC-derived dopaminergic neurons from patients with Gaucher disease and parkinsonism. J Neurosci. 2016; 36:7441–7452
20. Siebert M, Sidransky E, Westbroek W. Glucocerebrosidase is shaking up the synucleinopathies. Brain. 2014; 137:1304–1322
21. Lwin A, Orvisky E, Goker-Alpan O, LaMarca ME, Sidransky E. Glucocerebrosidase mutations in subjects with parkinsonism. Mol Genet Metab. 2004; 81:70–73
22. Kopustinskiene DM, Bernatoniene J. Molecular mechanisms of melatonin-mediated cell protection and signaling in health and disease. Pharmaceutics. 2021; 13:129
23. Druzin MY, Kurzina NP, Malinina EP, Kozlov AP. The effects of local application of D2 selective dopaminergic drugs into the medial prefrontal cortex of rats in a delayed spatial choice task. Behav Brain Res. 2000; 109:99–111
24. Facci L, Stevens DA, Pangallo M, Franceschini D, Skaper SD, Strijbos PJ. Corticotropin-releasing factor (CRF) and related peptides confer neuroprotection via type 1 CRF receptors. Neuropharmacology. 2003; 45:623–636
25. Brandies R, Yehuda S. The possible role of retinal dopaminergic system in visual performance. Neurosci Biobehav Rev. 2008; 32:611–656
26. Biehlmaier O, Alam M, Schmidt WJ. A rat model of Parkinsonism shows depletion of dopamine in the retina. Neurochem Int. 2007; 50:189–195
27. Tran T, Andersen R, Sherman SP, Pyle AD. Insights into skeletal muscle development and applications in regenerative medicine. Int Rev Cell Mol Biol. 2013; 300:51–83
28. Bartell PA, Miranda-Anaya M, McIvor W, Menaker M. Interactions between dopamine and melatonin organize circadian rhythmicity in the retina of the green iguana. J Biol Rhythms. 2007; 22:515–523
29. Meng T, Zheng ZH, Liu TT, Lin L. Contralateral retinal dopamine decrease and melatonin increase in progression of hemiparkinsonium rat. Neurochem Res. 2012; 37:1050–1056

corticotropin-releasing hormone; glucocerebrosidase gene; MTNR1B Parkinson’s disease

Supplemental Digital Content

Copyright © 2021 The Author(s). Published by Wolters Kluwer Health, Inc.