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

Glucocerebrosidase (GBA) mutations occur frequently in Parkinson's disease (PD) patients. This study aims to identify potential crucial genes and pathways associated with GBA mutations in patients with PD and to further analyze new molecular mechanisms related to the occurrence of gene mutations from the perspective of bioinformatics. Gene expression profiles of datasets GSE53424 and GSE99142 were acquired from the Gene Expression Ominibus database. Differentially expressed genes (DEGs) were detected, using the 'limma' package in R, comparing IDI-PD 1 (idiopathic PD patients) and GBA-PD 1 [PD patients with heterozygous GBA mutations (GBA N370S)] group samples. The functions of top modules were assessed using the DAVID, whereas gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed. Protein-protein interaction networks were assembled with Cytoscape software and separated into subnetworks using the Molecular Complex Detection Algorithm. Data from GSE53424 and GSE99142 were also extracted to verify our findings. There were 283 DEGs identified in PD patients heterozygous for GBA mutations. Module analysis revealed that GBA mutations in PD patients were associated with significant pathways, including Calcium signaling pathway, Rap1 signaling pathway and Cytokine-cytokine receptor interaction. Hub genes of the two modules were corticotropin-releasing hormone (CRH) and Melatonin receptor 1B (MTNR1B). The expression of CRH was downregulated, whereas that of MTNR1B was upregulated in PD patients with GBA mutations. The expression of CRH and MTNR1B has diagnostic value for PD patients with heterozygous GBA mutations. Novel DEGs and pathways identified herein might provide new insights into the underlying molecular mechanisms of heterozygous GBA mutations in PD patients.

Introduction

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.

Methods

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 (https://www.ncbi.nlm.nih.gov/geo/). 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 (https://david.ncifcrf.gov/) 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; http://string-db.org), 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; http://www.cytoscape.org/). 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.

Results

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

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).

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.

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
		Upregulated
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
		Downregulated
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

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

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, http://links.lww.com/WNR/A639. 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

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.

Discussion

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.

Limitations

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.

Acknowledgments

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.