Identification of key genes and pathways in Parkinson's disease through integrated analysis.

Parkinson's disease (PD) is a progressive, degene-rative neurological disease, typically characterized by tremors and muscle rigidity. The present study aimed to identify differe-ntially expressed genes (DEGs) between patients with PD and healthy patients, and clarify their association with additional biological processes that may regulate factors that lead to PD. An integrated analysis of publicly available Gene Expression Omnibus datasets of PD was performed. DEGs were identified between PD and normal blood samples. Gene Ontology enrichment and Kyoto Encyclopedia of Genes and Genomes pathway analyses, as well as protein‑protein interaction (PPI) networks were used to predict the functions of identified DEGs. Reverse transcription‑quantitative polymerase chain reaction (RT‑qPCR) was performed to validate the predicted expression levels of identified DEGs in whole blood samples obtained from patients with PD and normal healthy controls. A total of 292DEGs were identified between the PD and normal blood samples. Of these, 156 genes were significantly upregulated and 136 genes were significantly downregulated in PD samples following integrated analysis of four PD expression datasets. The 10 most upregulated and downregulated genes were used to construct a PPI network, where ubiquitin C‑terminal hydrolase L1 (UCHL1), 3‑phosphoinositide dependent protein kinase 1 (PDPK1) and protein kinase cAMP‑activated catalytic subunit β (PRKACB) demonstrated the highest connectivity in the network. DEGs were significantly enriched in amoebiasis, vascular smooth muscle contraction, and the Wnt and calcium signaling pathways. The expression levels of significant DEGs, UCHL1, PDPK1 and PRKACB were validated using RT‑qPCR analysis. The findings revealed that UCHL1 and PDPK1 were upregulated and PRKACB was downregulated in patients with PD when compared with normal healthy controls. In conclusion, the results indicate that the significant DEGs, including UCHL1, PDPK1 and PRKACB may be associated with the development of PD. In addition, these factors may be involved in various signaling pathways, including amoebiasis, vascular smooth muscle contraction and the Wnt and calcium signaling pathways.

Introduction

Parkinson's disease (PD) is a chronic, progressive, degenerative neurological disease, and is the second most common neurodegenerative disease following Alzheimer's disease (1). The clinical symptoms of PD may be motor-associated or patients may present with non-motor manifestations. Typical motor-associated manifestations include tremors, muscle rigidity, deterioration of voluntary movements and postural instability (2,3). The non-motor manifestations may include depression, autonomic dysfunction, cataracts and cognitive impairment, such as mild cognitive impairment and Parkinson's dementia (4–6).

The pathological features of PD include loss of midbrain dopaminergic neurons in the substantia nigra and intracellular inclusions of Lewy bodies predominantly composed of α-synuclein (7). Despite fairly typical clinical and pathological features, the etiology of PD remains to be elucidated. Numerous cellular mechanisms, including oxidative stress, mitochondrial dysfunction, neuroinflammatory abnormalities and cell apoptosis have been previously reported as potential causes (8).

Oxidative stress is a key pathogenic factor in PD. A previous study reported that the aggregation of α-synuclein forms to soluble oligomers may interact with metal ions to induce aberrant free radical production and neuronal toxicity, leading to oxidative stress and neuronal death in PD (9). In addition, mitochondrial dysfunction may be closely associated with increased reactive oxygen species formation in PD and affects various cellular pathways, leading to damage of intracellular components and cell death (10). α-synuclein is associated with the inner mitochondrial membrane, and may interfere with mitochondrial function by interacting with complex I and reducing mitochondrial complex I activity, which leads to increased free radical production and increased mitophagy in dopaminergic neurons (11,12).

Neuroinflammation is recognized as a pathophysiological contributor of PD. High levels of proinflammatory mediators, such as tumor necrosis factor-α, interleukin-1β and interferon-γ have been detected in the midbrain of patients with PD and contribute to the acceleration of nigral dopamine neuron degeneration (13,14).

Fyn, a non-receptor tyrosine kinase, maybe involved in progressive neurodegenerative processes in PD by regulating microglial neuroinflammatory responses in PD in vitro and in vivo (15).

The pathogenesis of PD remains to de elucidated and currently available treatments are limited. Current treatment options for PD include levodopa, dopamine receptor agonists, acetylcholinesterase inhibitors, adenosine 2A receptor antagonists, catechol-O-methyl transferase inhibitors and monoamine oxidase B-inhibitors, however those options have limited efficacy for patients with PD (16). The present study performed bioinformatics analysis to integrate the mRNA expression data from blood samples obtained from patients with PD and normal healthy controls in the Gene Expression Omnibus (GEO) database, in order to identify differentially expressed genes (DEGs). The aim of the present study was to provide valuable information for the identification of novel therapeutic targets for PD, in order to facilitate the development of effective targeted treatments.

Materials and methods

Gene expression datasets

The raw gene expression profiling datasets submitted by four individual studies (GEO accession nos. GSE54536, GSE34287, GSE18838 and GSE6613; Table I) (17–20) were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo) (21). In order to identify differentially expressed genes in PD, four GEO datasets were downloaded and used for integrated analysis. Microarray studies investigating PD and normal blood samples were identified. Non-human studies, review articles and studies involving the integrated analysis of expression profiles were not used for the purposes of the present study.

Table I.

Details of the RNA expression profiles employed in the present study.

Author, year	GEO ID	Platform	Samples (N:P)	Refs.
Alieva et al, 2014	GSE54536	Illumina HumanHT-12 v4.0 Expression BeadChip	5:5	17
Potashkin et al, 2011	GSE34287	ExonHit Human Genome Wide SpliceArray 1.0	20:19	18
Shehadeh et al, 2010	GSE18838	Affymetrix Human Exon 1.0 ST Array	12:18	19
Scherzer et al, 2007	GSE6613	Affymetrix Human Genome U133A Array	23:50	20

GEO, Gene Expression Omnibus database; N, normal blood samples; P, blood samples from patients with Parkinson's disease.

Differential expression analysis of genes

Complete datasets from the GEO database were downloaded and preprocessed using log2 transformation and Z-score normalization in R (version 3.2.3; https://www.r-project.org/). The limma package (version 3.30.13; https://www.bioconductor.org/), which is a linear model for microarray data analysis (22), was used to identify DEGs by comparing the expression profiles in blood samples obtained from patients with PD with those of normal healthy controls in each individual dataset. Genes that demonstrated a P-value of <0.01 were considered as DEGs and were selected for further investigation.

Functional annotation of DEGs

The underlying function of DEGs was predicted using the Gene Ontology (GO) enrichment function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (23,24). A P-value of <0.05 and a false discovery rate (FDR) of <0.05 were set as the threshold values for selecting significantly enriched functional GO terms and KEGG pathways, respectively. The GENECODIS 3 online software (http://genecodis.cnb.csic.es/analysis) was used for this analysis (25–27).

Construction of protein-protein interaction (PPI) networks

Genes encoding closely interacting proteins are involved in mediating complex cellular functions. The Biological General Repository for Interaction Datasets (BioGRID, version 3.4), which is a database of known and predicted protein interactions (28), was used to screen interacting protein pairs from the identified DEGs. A PPI network involving the 10 most upregulated and downregulated DEGs was visualized using Cytoscape 3.1.0 software (http://cytoscape.org/) (29).

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

A total of 5 ml fresh whole blood samples were obtained from patients with PD (n=6; 3 male and 3 female; age range, 52–63 years) and normal healthy controls (n=4; 2 male and 2 female; age range, 54–64 years). Subjects were recruited from December 2015 to February 2016, in the Department of Neurology, Liaocheng People's Hospital, (Liaocheng, China). White cells were separated from whole blood samples with centrifugation at 200 × g for 10 min at room temperature, and total white cell RNA was extracted using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA). All patients provided written informed consent. The SuperScript III Reverse Transcription kit (Invitrogen; Thermo Fisher Scientific, Inc.) was used to synthesize cDNA from 1 µg RNA. RT-qPCR reactions were performed using Power SYBR Green PCR Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.), AmpliTaq Gold DNA polymerase (Applied Biosystems; Thermo Fisher Scientific, Inc.), and the Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). β-actin was used as an internal control for quantification of target gene expression. The relative expression of genes was calculated using the comparative 2−∆∆Cq method (30). The following PCR primers were used: Ubiquitin C-terminal hydrolase L1 (UCHL1), forward, 5′-GAAGGCCAATGTCGGGTAGAT-3′, and reverse, 5′-AAAGGCATTCGTCCATCAAGTT-3′; 3-phosphoinositide dependent protein kinase 1 (PDPK1), forward, 5′-AGAGGTCAGGCAGCAACATAGAG-3′ and reverse, 5′-TAAACCCTTCCGCTTATCCACT-3′; protein kinase cAMP-activated catalytic subunit β (PRKACB), forward, 5′-TCGAGTACCTCCATTCACTAGACCT-3′, and reverse, 5′-CTGCCTTATTGTAGCCCTTGC-3′; β-actin, forward, 5′-ACTTAGTTGCGTTACACCCTT-3′, and reverse, 5′-GTCACCTTCACCGTTCCA3′. The following thermocycling conditions were used: Polymerase activation step at 95°C for 10 min, followed by 45 cycles of 2-step qPCR (15 sec of 95°C denaturation, 60 sec of 60°C combined annealing and extension).

Statistical analysis

At least three independent experiments were performed for statistical evaluation. RT-qPCR experimental data are expressed as the mean ± standard deviation. Statistical significance was determined using a Student's t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

DEGs in PD

A total of 4 mRNA expression profiles were collected (GEO accession nos. GSE54536, GSE34287, GSE18838 and GSE6613), which included the blood samples obtained from 92 patients with PD and 60 normal healthy controls (Table I). Following normalization of the raw microarray data, 292 significant DEGs were identified inpatients with PD when compared with healthy controls. These consisted of 156 upregulated and 136 downregulated genes. Receptor interacting serine/threonine kinase 4 was the most significantly upregulated gene and mannose-P-dolichol utilization defect 1 was the most significantly downregulated gene (Table II). The 10 most significantly upregulated and downregulated genes are presented in Table II. The pattern of gene expression alterations for the 50 most significant DEGs is displayed as a heat map in Fig. 1.

Table II.

Significantly dysregulated genes in Parkinson's disease.

A, Top 10 upregulated genes
NCBI gene ID	Gene symbol	P-value
54101	RIPK4	6.22×10−7
6369	CCL24	2.42×10−5
5170	PDPK1	3.15×10−5
54328	GPR173	5.77×10−5
8111	GPR68	1.46×10−4
57214	KIAA1199	1.92×10−4
3699	ITIH3	2.05×10−4
55267	C22orf26	2.07×10−4
7345	UCHL1	2.47×10−4
58529	MYOZ1	2.71×10−4
B, Top 10 downregulated genes
NCBI gene ID	Gene symbol	P-value
9526	MPDU1	4.62×10−6
348995	NUP43	4.91×10−5
5567	PRKACB	5.90×10−5
79830	ZMYM1	8.01×10−5
10195	ALG3	1.05×10−4
10396	ATP8A1	1.30×10−4
2235	FECH	2.28×10−4
80169	C17orf68	2.59×10−4
200081	TXLNA	2.70×10−4
10178	ODZ1	2.84×10−4

NCBI, national centre for biotechnology information.

Figure 1.

A heat map depicting alterations in the expression patterns of the 50 most differentially expressed genes in patients with Parkinson's disease across four different datasets (Gene Expression Omnibus accession nos. GSE18838, GSE34287, GSE54536 and GSE6613).

GO classification of DEGs

DEGs in PD were annotated using GO database to gain insight into their putative biological roles. The threshold of GO terms was defined as P<0.05. Enrichment analysis of the most significant biological processes is presented in Table IIIA, and included asparaginyl-tRNA aminoacylation (GO: 0006421; P=2.32×10−7), cytotoxic T cell degranulation (GO: 0043316; P=2.32×10−7) and porphyrin-containing compound metabolic processes (GO: 0006778; P=9.88×10−6). Enrichment analysis of the most significant molecular functions is presented in Table IIIB, and included ferrochelatase activity (GO: 0004325; P=2.32×10−7), asparagine-tRNA ligase activity (GO: 0004816; P=1.83×10−7) and acid phosphatase activity (GO: 0003993; P=1.35×10−6). Enrichment analysis of the most significant cellular components is presented in Table IIIC, and included eukaryotic translation initiation factor 3 complex only (GO: 0005852; P=8.36×10−5).

Table III.

GO function enrichment analysis of differentially expressed genes in Parkinson's disease.

A, Biological process
GO ID	GO term	Count	P-value
0006778	Porphyrin-containing compound metabolic process	12	9.88×10−6
0006421	Asparaginyl-tRNAaminoacylation	8	2.32×10−7
0001913	T cell mediated cytotoxicity	6	5.58×10−5
0043316	Cytotoxic T cell degranulation	6	2.32×10−7
0002323	Natural killer cell activation involved in immune response	6	6.98×10−6
043320	Natural killer cell degranulation	6	6.98×10−6
0006611	Protein export from nucleus	11	1.45×10−6
0001732	Formation of translation initiation complex	6	1.45×10−5
0051084	‘De novo’ post-translational protein folding	11	6.37×10−5
0051085	Chaperone mediated protein folding requiring cofactor	11	4.49×10−8
0061077	Chaperone-mediated protein folding	11	1.04×10−6
B, Molecular function
GO ID	GO term	Count	P-value
0004325	Ferrochelatase activity	6	2.32×10−7
0004816	Asparagine-tRNA ligase activity	8	1.83×10−7
0003993	Acid phosphatase activity	10	1.35×10−6
0031489	Myosin V binding	6	2.98×10−5
0003989	Acetyl-CoA carboxylase activity	8	3.80×10−7
0016421	CoA carboxylase activity	8	2.44×10−6
0016885	Ligase activity, forming carbon-carbon bonds	8	9.88×10−6
0008113	Peptide-methionine (S)-S-oxide reductase activity	5	5.45×10−6
0004809	tRNA (guanine-N2-)-methyltransferase activity	7	1.00×10−7
0016423	tRNA (guanine) methyltransferase activity	7	3.05×10−6
0005540	Hyaluronic acid binding	12	4.03×10−7
C, Cellular component
GO ID	GO term	Count	P-value
0005852	Eukaryotic translation initiation factor 3 complex	11	8.36×10−5

GO, gene ontology.

Pathway enrichment analysis

KEGG pathway enrichment analysis was performed for the identified DEGs. FDR values of <0.05 was used as the threshold for pathway detection. The pathways that were most highly enriched were amoebiasis (FDR=6.92×10−4), vascular smooth muscle contraction (FDR=7.24×10−4) and the Wnt signaling pathway (FDR=8.17×10−4; Table IV).

Table IV.

Top 15 KEGG pathways obtained from enrichment analysis of differentially expressed genes in Parkinson's disease.

KEGG ID	KEGG term	Count	FDR	Genes
hsa05146	Amoebiasis	7	6.92×10−4	COL4A4, IL1R2, PLCB2, COL4A1, ARG1, NFKB1, PRKACB
hsa04270	Vascular smooth muscle contraction	4	7.24×10−4	CACNA1D, PLCB2, ADCY7, PRKACB
hsa04310	Wnt signaling pathway	8	8.17×10−4	PPP2R5B, AXIN1, WNT16, PLCB2, MYC, TBL1XR1, PPP3CB, PRKACB
hsa04020	Calcium signaling pathway	5	8.37×10−4	CACNA1D, PLCB2, ADCY7, PPP3CB, PRKACB
hsa04010	MAPK signaling pathway	10	8.44×10−4	NFKB2, IL1R2, CACNA1D, NFKB1, MYC, MAPKAPK3, PPP3CB, PRKACB, HSPA8, FLNB
hsa04724	Glutamatergic synapse	8	8.67×10−4	GRIA3, GRM8, CACNA1D, PLCB2, ADCY7, ADRBK2, PPP3CB, PRKACB
hsa05200	Pathways in cancer	11	9.51×10−4	COL4A4, NFKB2, AXIN1, WNT16, TPM3, RUNX1, CSF1R, COL4A1, PML, NFKB1, MYC
hsa04970	Salivary secretion	6	9.62×10−4	ATP1B1, TRPV6, PLCB2, ADCY7, PRKG2, PRKACB
hsa04640	Hematopoietic cell lineage	6	9.82×10−4	CD4, IL1R2, CSF1R, GP1BA, CD1C, CD44
hsa04540	Gap junction	4	1.02×10−3	PLCB2, ADCY7, PRKG2, PRKACB
hsa04912	GnRH signaling pathway	3	1.32×10−3	CACNA1D, ADCY7, PRKACB
hsa04062	Chemokine signaling pathway	8	1.39×10−3	CCL24, ITK, PLCB2, ADCY7, NFKB1, ADRBK2, CXCL12, PRKACB
hsa04920	Adipocytokine signaling pathway	5	1.57×10−3	G6PC2, AGRP, ACACB, NFKB1, PRKAG2
hsa04916	Melanogenesis	3	1.60×10−3	PLCB2, ADCY7, PRKACB
hsa04660	T cell receptor signaling pathway	6	1.70×10−3	CD4, ITK, NFKB1, GRAP2, PPP3CB, CD28

KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate.

PPI network construction

PPI networks of the 10 most upregulated and downregulated DEGs were constructed using Cytoscape software. The network consisted of 216 nodes and 231 edges. In the PPI network, the nodes with a high degree were classified as hub proteins, and degrees were defined to quantify the number of neighbors that a node demonstrated a direct connection with. As shown in Fig. 2, the red nodes represent upregulated genes and the green nodes represent downregulated genes. The most significant hub proteins were UCHL1, PDPK1 and PRKACB (Fig. 2).

Figure 2.

Protein-protein interaction networks of the 10 most significantly upregulated and downregulated DEGs. Nodes represent proteins and edges represent interactions between two proteins. Red nodes represent upregulated DEGs, whereas and green nodes indicate downregulated DEGs. Light-blue nodes represent protein products of genes predicted to interact with the DEGs. The solid lines indicate protein-protein interaction correlations, and dashed lines represent protein-protein co-localization. DEGs, differentially expressed genes.

RT-qPCR validation of differential expression mRNAs and target genes

In order to validate the microarray analysis data, the expression levels of significant DEGs (UCHL1, PDPK1, PRKACB) were quantified by RT-qPCR analysis in 6PD and 4 normal blood control samples. No significant difference between the expression levels of UCHL1 and PDPK1 in the PD and normal blood samples was observed (Fig. 3A and B); however, the expression levels of these genes were higher in the PD samples when compared with the controls. Similarly, no significant difference in the expression levels of PRKACB between PD and normal control samples was observed; however, PRKACB expression was lower in the PD group when compared with the controls (Fig. 3C).

Figure 3.

mRNA expression levels of (A) UCHL1 (B) PDPK1 and (C) PRKACB in PD and CON blood samples, as determined using reverse transcription-quantitative polymerase chain reaction analysis. UCHL1, ubiquitin C-terminal hydrolase L1; PDPK1, 3-phosphoinositide dependent protein kinase 1; PRKACB, protein kinase cAMP-activated catalytic subunit b; CON, normal blood samples; PD, blood samples from patients with Parkinson's disease.

Discussion

The results of the present study demonstrated that UCHL1 was upregulated in PD samples, and was one of the hub proteins in the PPI network. UCHL1 encodes the 223-amino acid ubiquitin carboxyl-terminal hydrolase-L1 enzyme, which is abundant in brain tissues. Variants of the UCHL1 gene are a risk factor for PD (31,32). UCHL1 is primarily expressed in neurons and has been implicated in various neurodegenerative diseases, such as PD (33). Additionally, the UCHL1 protein has been identified in Lewy bodies (34). A previous study revealed that UCHL1 may be essential for the pathogenesis of PD (35). UCHL1 was associated with α-synuclein accumulation in synaptic vesicles (36). The quantity of membrane-associated UCHL1 may affect the intracellular α-synuclein levels and neurotoxicity. Chemical inhibition of UCHL1 by farnesylation may reduce α-synuclein levels and improve neuronal cell viability in cellular models of α-synuclein-associated toxicity for PD (37). Notably, the UCHL1 S18Y polymorphism has been previously associated with increased susceptibility for PD (31,38).

The PDPK1 gene encodes the PDPK1 protein. Previous studies have reported that PDPK1 serves an important role in various neurodegenerative diseases, including Alzheimer's and Huntington's disease (39,40). The phosphatidylinositol-4, 5-bisphosphate 3-kinase/phosphoinositide-dependent kinase (PDK)/nuclear protein kinase C/Ras-related C3 botulinum toxin pathway axis is required by the 1–42 β-amyloid peptide to induce neuronal death, and may be involved in the development of Alzheimer's disease (41,42). Huntington's disease is caused by the dominant heritable expansion of atrinucleotide CAG repeat in the huntingtin gene (43). Increased levels of active PDK1 and reduced levels of extracellular-signal regulated kinase1/2 provides neuroprotection in striatal cells during Huntington's disease (44). However, the function of PDK1 in PD remains to be elucidated. The present study determined that PDK1 was upregulated in PD, and that PDK1 may therefore be involved in the development of PD.

PRKACB is a member of the serine/threonine protein kinase family and mediates signaling via cyclic adenosine monophosphate. The two catalytic subunit alleles of protein kinase A (PKA), PKACα and PKACβ, are encoded by the PRKACA and PRKACB genes, respectively (45). L-dopa is an effective treatment for patients with PD; however, long-term therapy is associated with complications, such as L-dopa-induced dyskinesia. It has been previously reported that levodopa/Benserazide may prevent L-dopa-induced dyskinesia by inactivating the DR1/PKA/P-τ pathway in a rat model of PD (46). In the present study, PRAKCB was observed to be upregulated in PD samples when compared with normal controls, which was in accordance with the integrated analysis of mRNA expression profiles. However, the biological role of PRAKCB in PD remains to be elucidated.

In the present study, the Wnt signaling pathway was identified to be one of the significant enrichment pathways of DEGs in PD. A previous study reported that the Wnt signaling pathway potentially regulates the expression levels of the key PD protein, α-synuclein (47). Dysregulated Wnt signaling serves an important role in the progression of PD through brain inflammation and impaired adult neurogenesis. Depletion of canonical Wnt signaling components has demonstrated a neuroprotective effect on midbrain dopaminergic neurons in a mouse model of PD (48).

In conclusion, the present study identified a total of 292DEGs in PD. The 10 most upregulated and downregulated genes were used for PPI network construction. In this network the UCHL1, PDPK1 and PRKACB genes were identified. The results suggest that these genes may serve an important role in PD via various signaling pathways, including amoebiasis, vascular smooth muscle contraction and the Wnt signaling cascade. These results may provide an important contribution to the identification of novel therapeutic targets for PD, as well as facilitate the development of effective targeted therapies. Further investigation into their biological functions may provide an additional insight into the role of the differentially regulated genes in the pathophysiology of PD.