Identification of Hub Genes and Potential Molecular Pathogenesis in Substantia Nigra in Parkinson's Disease via Bioinformatics Analysis

Parkinson's disease (PD) is the second most common neurodegenerative disease, with significant socioeconomic burdens. One of the crucial pathological features of PD is the loss of dopaminergic neurons in the substantia nigra (SN). However, the exact pathogenesis remains unknown. Moreover, therapies to prevent neurodegenerative progress are still being explored. We performed bioinformatics analysis to identify candidate genes and molecular pathogenesis in the SN of patients with PD. We analyzed the expression profiles, GSE49036 and GSE7621, which included 31 SN tissues in PD samples and 17 SN tissues in healthy control samples, and identified 86 common differentially expressed genes (DEGs). Then, GO and KEGG pathway analyses of the identified DEGs were performed to understand the biological processes and significant pathways of PD. Subsequently, a protein-protein interaction network was established, with 15 hub genes and four key modules which were screened in this network. The expression profiles, GSE8397 and GSE42966, were used to verify these hub genes. We demonstrated a decrease in the expression levels of 14 hub genes in the SN tissues of PD samples. Our results indicated that, among the 14 hub genes, DRD2, SLC18A2, and SLC6A3 may participate in the pathogenesis of PD by influencing the function of the dopaminergic synapse. CACNA1E, KCNJ6, and KCNB1 may affect the function of the dopaminergic synapse by regulating ion transmembrane transport. Moreover, we identified eight microRNAs (miRNAs) that can regulate the hub genes and 339 transcription factors (TFs) targeting these hub genes and miRNAs. Subsequently, we established an mTF-miRNA-gene-gTF regulatory network. Together, the identification of DEGs, hub genes, miRNAs, and TFs could provide better insights into the pathogenesis of PD and contribute to the diagnosis and therapies.


Introduction
Parkinson's disease (PD) is the second most common neurodegenerative disease, with substantial socioeconomic burdens [1]. Te prevalence of global PD increased by 155.51% during 1990-2019 [2]. First described by James Parkinson in 1817, PD is characterized by a range of motor symptoms, including bradykinesia, resting tremor, rigidity, and posture instability. In addition, certain non-motor symptoms are also shown to be associated with PD, such as sleep disturbances, autonomic dysfunctions, cognitive and psychiatric dysfunctions, and sensory symptoms [3]. Currently, no available therapies can cause efective prevention of neurodegenerative progress. With the progression of the disease, symptoms that respond poorly to treatment may severely afect the quality of life.
Studies have established that the etiology of PD is multifactorial, which could include genetic factors, environmental factors, nervous system aging, and other factors [4,5]. Te selective loss of dopaminergic neurons in substantia nigra (SN) and the appearance of Lewy bodies in the cytoplasm of the remaining neurons are considered predominant pathological features of PD. In 1893, Bloq and Marinesco proposed the role of SN in the pathological development of PD. In 1912, Friedrich Heinrich Lewy identifed the presence of intraneuronal inclusions, now known as Lewy bodies, in the remaining neurons in patients of PD [6]. Studies have shown that the dopaminergic neurons in the ventrolateral tier of the SN are preferentially lost [7]. Moreover, 40-60% of dopaminergic neurons in the nigrostriatal system are lost before the frst appearance of the motor symptoms [8]. Te loss of dopaminergic neurons in SN could be regulated by several factors, such as oxidative stress, proteasome dysfunction, mitochondrial dysfunction, infammatory and immune response, apoptosis, and other mechanisms [9,10]. However, the exact pathogenesis of PD is still unclear. Tus, it is necessary and urgent to examine the pathogenesis of PD and fnd efective diagnosis and treatment strategies.
In recent years, RNA sequencing and microarray have become indispensable tools for identifying the expression of diferential genes, mRNAs, and non-coding RNAs [11]. Furthermore, several datasets of genes expressed in the SN of patients with PD can be downloaded from public databases, such as the Gene Expression Omnibus (GEO). Te development of bioinformatics analysis has made it possible to screen, compare, and analyze the existing data and identify the diferentially expressed genes (DEGs) that may be related to PD. Systematic studies on the relationship between DEGs could determine the biological processes involved in PD. Terefore, bioinformatics analysis helps explore the biological mechanism and potential biomarkers of PD. Te existing bioinformatics research related to PD shows several genes, mRNAs, and non-coding RNAs that may be associated with PD. For example, SLC6A3 is critical in maintaining the integrity of dopaminergic neurons, while SLC18A2 is essential for their survival by contracting intracellular toxicity [12]. ICAM1, also known as CD54, may increase neprilysin levels, essential to treat neurological diseases, including PD [13]. Studies have shown that HRAS may be related to L-DOPA-induced dyskinesia and cognitive impairment [14]. Similarly, miR-338 can decrease mitochondrial activity by reducing the cytochrome c oxidase IV, leading to neuronal damage in SN [15]. Tese fndings using bioinformatics analysis signifcantly contribute to our understanding of the causes and underlying molecular events of PD. However, additional studies are required to gain a more accurate understanding of PD pathogenesis.
Our study utilized bioinformatics analysis to explore the hub genes and potential molecular mechanisms in the SN of patients with PD, reveal the pathogenesis of PD, fnd diagnostic markers and therapeutic targets, and provide new perspectives and strategies for the diagnosis, treatment, and new drug development of PD. We identifed 86 common DEGs, constructed a protein-protein interaction (PPI) network, and selected 15 hub genes. To understand these genes better, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. Moreover, we identifed eight micro-RNAs (miRNAs) that could regulate the hub genes. By targeting these hub genes and miRNAs, we further identifed 339 transcription factors (TFs), thereby establishing an mTF-miRNA-gene-gTF regulatory network. Identifying DEGs, hub genes, miRNAs, and TFs could provide insights into the pathogenesis of PD and contribute to further diagnosis and therapies.

Microarray Data Analysis.
Te gene expression profles, GSE49036 and GSE7621, for the SN of patients with PD were obtained from the GEO database (https://www. ncbi.nlm.nih.gov/geo/), an international public repository. Microarray data of GSE49036 and GSE7621 were based on the GPL570 Platform (Afymetrix Human Genome U133 Plus 2.0 Array); GSE49036 included 15 SN tissues of PD samples and eight SN tissues of normal samples, while GSE7621 contained 16 SN tissues of PD samples and 9 SN tissues of normal samples. Additionally, datasets GSE8397 and GSE42966 were obtained from the GEO database to verify hub genes selected from GSE49036 and GSE7621. Microarray data for GSE8397 were based on the GPL96 Platform (Afymetrix Human Genome U133A Array) and GPL97 Platform (Afymetrix Human Genome U133B Array). A and B GeneChip of GSE8397 together contained 24 SN tissues of PD samples and 15 SN tissues of normal samples. Finally, the microarray data for GSE42966 were based on the GPL4133 Platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F) and included 9 SN tissues of PD samples and 6 SN tissues of normal samples (Table 1).

GO and KEGG Pathway Analyses of DEGs.
GO analysis is one of the most valuable methods for describing the features of genes comprehensively, which includes biological processes (BPs), molecular functions (MFs), and cellular components (CCs). Similarly, the KEGG database explores the functions and biological pathways of genes. GO and KEGG pathway analyses of the overlapping DEGs between GSE49036 and GSE7621 were performed using an online tool called DAVID (https://david.ncifcrf.gov/). Finally, the ggplot2 package in R Studio was applied to depict the bubble plots.

PPI Network Construction and Hub Gene Selection.
Te Search Tool for the Retrieval of Interacting Genes (STRING (https://string-db.org/)), together with Cytoscape software, was used to build a PPI network (medium confdence: 0.4). MCODE, a plugin in Cytoscape, determined the signifcant clusters, while the plugin Cyto-Hubba was used to screen the hub genes in the PPI network. GO and KEGG pathway analyses of hub genes and genes in clusters were also predicted by DAVID.

Prediction of Target miRNAs.
Two online miRNA databases, TargetScan (https://www.targetscan.org/vert_80/) and miRDB (http://www.mirdb.org/), were used to predict the miRNA targeting hub genes. Te intersection of miRNAs obtained above and the diferentially expressed miRNAs in SN between PD and normal individuals were obtained using the Venn online website. Tese miRNAs were considered target miRNAs in the study.

2.6.
mTF-miRNA-Gene-gTF Regulatory Network Construction. To further explore the functions of the abovefound hub genes and target miRNAs in PD pathogenesis, TFs related to hub genes (gTF) and TFs related to target miRNAs (mTF) were predicted by the online database RNAInter (http://www.rnainter.org/). Finally, an mTF-miRNA-gene-gTF regulatory network was established by Cytoscape software.

Identifcation of DEGs.
Te concise diagram of workfow is summarized in Figure 1. We selected two datasets, GSE49036 (15 SN tissues of PD samples and eight SN tissues of normal samples) and GSE7621 (16 SN tissues of PD samples and nine SN tissues of normal samples), for our study. Based on the criteria of P value <0.05 and |log FC| ≥ 1.0, we obtained 253 DEGs (29 upregulated genes and 224 downregulated genes) from GSE49036 and 1236 DEGs (732 upregulated genes and 504 downregulated genes) from GSE7621. We used volcano plots to visualize the DEGs in GSE49036 and GSE7621 (Figures 2(a) and 2(b)). Subsequently, we identifed 86 overlapping DEGs (2 upregulated genes and 84 downregulated genes) between GSE49036 and GSE7621 ( Figure 2(c), Supplementary Table 1).

GO and KEGG Pathway Analyses of DEGs.
For a comprehensive understanding of the DEGs, we performed GO and KEGG pathway analyses (P value <0.05) using DAVID. Te results of the GO analysis are presented in Table 2 and Figures 3(a)-3(c). For BP, these common DEGs were signifcantly enriched in chemical synaptic transmission, protein localization to the plasma membrane, homophilic cell adhesion via plasma membrane adhesion molecules, response to xenobiotic stimulus, axon guidance, regulation of ion transmembrane transport, dopaminergic neuron diferentiation, adult locomotory behavior, neurotransmitter transport, positive regulation of synapse assembly, and exocytosis. For CC, the common DEGs were mainly enriched in integral components of the plasma membrane, axon, dendrite, synapse, neuron projection, cell surface, glutamatergic synapse, and neuronal cell body. For MF, the common DEGs were enriched in calcium ion binding, protein N-terminus binding, ion channel binding, dopamine binding, monoamine transmembrane transporter activity, and high voltage-gated calcium channel activity. KEGG pathways were mainly enriched in the calcium signaling pathway, dopaminergic synapse, synaptic vesicle cycle, adrenergic signaling in cardiomyocytes, cocaine addiction, longevity regulating pathway-multiple species, retinol metabolism, and arrhythmogenic right ventricular cardiomyopathy (Table 3 and Figure 3(d)).

PPI Network Construction and Hub Gene Identifcation.
Using the data obtained above, we constructed a PPI network consisting of 49 nodes and 69 edges (Figure 4(a)). Te MCODE plugin in Cytoscape generated four modules (Figures 4(b)-4(e)); cluster 1, comprising 7 nodes and 19 edges, got the highest score (score: 6.333), while cluster 2 (7 nodes and 9 edges), cluster 3 (3 nodes and 3 edges) and cluster 4 (3 nodes and 3 edges) had the same score (score: 3). Te results of GO and KEGG analyses for these four modules suggest that DEGs are enriched mainly in transmembrane transport, chemical synaptic transmission, the biological process of locomotory behavior, and PD (P value <0.05; Supplementary Tables 2 and 3). However, we did not fnd any hits in the MF analysis for cluster 4 and KEGG pathway analysis for clusters 2, 3, and 4.

Validation of the Hub Genes.
To reinforce the reliability of hub genes in our study, we verifed the 15 hub genes in GSE8397 and GSE42966 and performed the box plots using GraphPad Prism 9 software ( Figure 6). Te expression levels of the 12 hub genes (SLC18A2, SLC6A3, KCNJ6, NR4A2, DRD2, RET, EN1, FGF13, SYNGR3, RIMBP2, KCNB1, and RAB3C) in SN tissues of PD samples were signifcantly decreased compared with those of normal samples (P value <0.01 and log FC ≤ −1.0). While the expression levels of FOXA2 and CACNA1E were also reduced in SN tissues of PD samples, they were not statistically signifcant (P value < 0.05 but −1< log FC < 0). On the other hand, the levels of UNC13C did not show any diference (P value >0.05). Together, our results indicated that the expression levels of 14 hub genes (SLC18A2, SLC6A3, KCNJ6, FOXA2, NR4A2, CACNA1E, DRD2, RET, EN1, FGF13, SYNGR3, RIMBP2, KCNB1, and RAB3C) were decreased in SN tissues of PD samples.

Prediction of Target TFs and Construction of mTF-miRNA-Gene-gTF Regulatory Network.
To better understand the 14 hub genes and eight target miRNAs found above, TFs targeting hub genes (gTF) and TFs targeting miRNAs (mTF) were identifed by the online database RNAInter. Moreover, we constructed an mTF-miRNAgene-gTF regulatory network using the software Cytoscape ( Figure 7). Tis network consisted of 208 nodes and 351 edges, involving six hub genes, eight target miRNAs, and 194 target TFs. Te top seven TFs of the network with the highest degrees (degree ≥5) were NHF4A, CDX2, FUS, E2F4, E2F6, ERG, and SUPT5H.

Discussion
Several eforts have been made to explore the pathogenesis of PD in recent years. However, the pathogenesis of PD is still unclear, and its efective therapy still needs more studies.
With the development of bioinformatics, microarray has become an indispensable tool to identify the expression of diferential genes, mRNAs, and non-coding RNAs. Moreover, several datasets of genes expressed in SN of PD have also been uploaded to the GEO database. In the past years,   Tis study explored the potential pathogenesis in SN of PD via bioinformatics analysis with diferent microarray datasets. We identifed 86 DEGs that were signifcantly enriched in various metabolism pathways. Further, we found that 14 hub genes, eight miRNAs, and seven TFs may play important roles in the pathogenesis of PD. GO and KEGG pathway analyses of hub genes suggest that the regulation of dopaminergic synaptic transmission might be involved in the pathogenesis of PD.
We found that three hub genes, DRD2, SLC18A2, and SLC6A3, in the PPI network were signifcantly enriched in the dopaminergic synapse, whereas SLC18A2 and SLC6A3 exhibited the highest degree of 9 and 8, respectively. DRD2 encodes the D2 subtype of the dopamine receptor. Recent studies show that dopamine agonists with high selectivity for DRD2 have already been used to improve symptoms in patients with PD [16]. DRD2 is related to peak-dose dyskinesias induced by levodopa in patients with PD [17]. Te DRD2 polymorphism, rs1076560 DRD2 G > T, might infuence gait function for patients with PD [18]. A clinical trial with 217 patients with PD on levodo pa therapy showed that DRD2 rs1799732 is an independent predictor of gastrointestinal symptoms associated with levodopa therapy [19]. SLC6A3 and SLC18A2 are the other two crucial candidate genes in sporadic PD. Encoded by SLC6A3, dopamine transporter (DAT) is the protein with the most selective expression of the most damaged dopaminergic neurons in patients with PD. Na + -K-ATPases on the plasma membrane can generate ion gradients. DAT reuptakes dopamine into presynaptic neurons from the synaptic cleft depending on the cotransport of Na+ and Cl− down the ions' concentration gradients [20]. Mainly present on the neuron terminals in SN, DAT is necessary for dopaminergic neurotransmission to control its intensity and duration [21]. A randomized trial by Moreau et al. showed that methylphenidate, an inhibitor of SLC6A3, can reduce the severity of gait hypokinesia and freezing in patients with advanced PD who received subthalamic nucleus stimulation [22]. SLC18A2 encodes vesicular monoamine transporter 2 (VMAT2), which can transport cytoplasmic monoamines into synaptic vesicles for storage, and then releases them extracellularly in the central nervous system, driven by H+ electrochemical force. Terefore, the concentrations of monoamine neurotransmitters in synaptic vesicles can maintain a high level, and those in the cytoplasm can maintain a low level. On the contrary, the decrease of     autopsies on six patients with PD and four healthy controls and gained dopamine storage vesicles from their striatum. Tey found that in patients with PD, the level of VMAT2 and synaptic vesicular dopamine uptake was signifcantly reduced, and dopamine storage impairment was located in the VMAT2 itself [25]. Interestingly, a decrease in vesicular function because of SLC18A2 mutation could lead to brain dopamine-serotonin vesicular transport disease, including infantile parkinsonism-dystonia-2. Several patients with brain dopamine-serotonin vesicular transport disease have been found in the world. Te homozygous c.710C > T (p.Pro237His) transition in SLC18A2 has been identifed in a 6-month-old male infant of China, two New Zealand siblings of European descent, and a 7-year-old female of Iraq [26][27][28]. Moreover, another variant, c.1160C⟶T in SLC18A2, has been observed in eight children in a Saudi Arabian family [29]. Tese observations are consistent with our results that the expression of DRD2, SLC18A2, and SLC6A3 in patients of PD is signifcantly reduced. So, we speculated that DRD2, SLC18A2, and SLC6A3 might participate in the pathogenesis of PD by infuencing the function of the dopaminergic synapse. GO analysis showed that CACNA1E, KCNJ6, and KCNB1 participated in the regulation of ion transmembrane   Figure 5: Interaction network of 15 hub genes. Te green color represents the downregulated genes screened based on fold change ≤ −1.0 and P value <0.05.  transport. Ion channels are proteins that generate and modulate electricity across biological membranes. CAC-NA1E encodes the high-voltage-activated Cav2.3 type R calcium channel. Voltage-gated Ca2+ channels, whose primary function is to initiate the synaptic transmission and neurotransmitter release, consist of 5 distinct subunits (α1, α2, β, c, and δ). Te α1 subunit can be divided into 3 subfamilies, namely, Cav1, Cav2, and Cav3 [30]. It has been reported that among all voltage-gated Ca2+ channel subtypes in adult SN dopaminergic neurons, Cav2.3 accounts for the most signifcant proportion. In SN dopaminergic neurons, the activity that generates oscillatory increases in free cytosolic Ca2+ levels, which are thought to impart mitochondrial stress and render these neurons more vulnerable to degeneration by PD stressors. In Cav2.3 knockout mice and Cav2.3 inhibitor SNX-482 using mice, the activityassociated nigral somatic Ca2+ signals reduced, which indicates that Cav2.3 contributes to neurodegeneration [31]. KCNJ6 encodes a potassium channel subunit called GIRK2, which belongs to the G-protein-gated inwardly rectifying potassium channel (GIRK) family. GIRK can be activated by ligand-stimulated G protein-coupled receptors (GPCRs), such as dopaminergic D2 receptors. Tus, the permeability of GIRK to K+ increases while the excitability of neurons decreases [32,33]. KCNB1 encodes an ion channel called the delayed rectifer voltage-gated K+ channel KCNB1 (Kv2.1).
Studies using the traumatic brain in mouse models indicated that the oxidation of Kv2.1 may cause neurodegeneration and cognitive impairment [34]. However, the specifc function and mechanism of KCNB1 in PD need further examination. Tese conclusions are consistent with our GO analysis results that CACNA1E and KCNJ6 are both enriched in the regulation of ion transmembrane transport and dopaminergic synapse, but KCNB1 is only enriched in the regulation of ion transmembrane transport. Terefore, our results suggest that CACNA1E, KCNJ6, and KCNB1 might infuence the function of the dopaminergic synapse by participating in the regulation of ion transmembrane transport. GO analysis in our study also showed that NR4A2, FOXA2, and EN1 are enriched in dopaminergic neuron diferentiation and adult locomotory behavior. NR4A2, also called Nurr1, is expressed predominantly in the central nervous system, especially in SN. NR4A2 is related to the diferentiation of the dopaminergic neurons in SN via activating the transcription of tyrosine hydroxylase and enhancing the expression of DAT [35][36][37]. Xu et al. performed heteroduplex analysis and sequencing analysis for the polymorphisms and mutations of NR4A2 in 225 patients with PD and 221 healthy control individuals. Tey found that a homozygous 7048G7049 polymorphism in intron 6 of the NR4A2 was higher in patients with PD than in healthy   [56] people. Teir analysis provides evidence for the association between NR4A2 and PD [38]. FOXA2 encodes a member of the forkhead class of DNA-binding proteins, which have been suggested to enhance the expression of the Nurr1induced DA phenotype [39]. EN1 encodes a homeodomain TF, which is necessary for the development of neurons in the midbrain, cerebellum, hindbrain, and spinal cord [40]. Simon et al. investigated the homeodomain TFs En-1 and En-2 in mice and found an increased En-1 in almost all dopaminergic neurons in SN and ventral tegmentum, but it was not required for their specifcation [41]. Tese studies further support the current fndings.
Additionally, our analysis showed that SYNGR3, RET, FGF13, RIMBP2, and RAB3C were downregulated in the SN of patients with PD and were the hub genes in the PPI network. SYNGR3 encodes an integral membrane protein called synaptogyrin-3, located on the synaptic vesicular membrane and involved in synaptic vesicular trafcking.
Te impairment of synaptic vesicular trafcking is one of the earliest pathological processes involved in PD [42]. In mice, synaptogyrin-3 interacts with DAT and increases its activity. Tis efect could be abolished in the presence of reserpine, a VMAT2 inhibitor [43]. Similar to our fndings, Simunovic et al. analyzed the gene expression profling of SN pars compacta in patients with PD and control individuals and found a reduction in the expression of SYNGR3 [44]. RET is a member of the cadherin superfamily. Glial cell line-derived neurotrophic factor (GDNF), which is one of the ligands of RET, can bind with the glycosylphosphatidylinositol-(GPI-) linked GDNF family receptor alpha 1 (GFRα1). Te combination of GDNF and GFRα1 can bind with RET to activate its intracellular tyrosine kinase activity [45]. Activated RET, in turn, activates the intracellular mitogen-activated protein kinase (MAPK), Akt (protein kinase B), and Src signaling cascades which can maintain the survival and regeneration of DA neurons [46]. Tus, GDNF is suggested to be one of  Figure 7: Te mTF-miRNA-gene-gTF regulatory network. Te square-shaped green nodes represent hub genes, the diamond-shaped yellow nodes represent miRNAs, the V-shaped blue nodes represent the TFs associated with hub genes (gTFs), and the triangle-shaped orange nodes represent the TFs associated with miRNAs (mTFs). the target-derived neurotrophic factors for the development of DA neurons and one of the factors that can maintain the survival of midbrain DA neurons [47,48]. Experiments in mice models confrmed that RET was crucial to the nigrostriatal DA system preservation. RET ablations in mice could lead to progressive loss of DA neurons and degeneration of DA nerve terminals in the striatum [49]. FGF13 encodes a protein that belongs to the fbroblast growth factor (FGF) family. Anatomic studies and electrophysiological recordings showed that FGF13 plays a powerful role in regulating excitability for hippocampal neurons [50]. RIMBP stands for Rab-interacting molecule-(RIM-) binding protein, whose interactions with Cav have been confrmed [51]. RIM-BP2 might have the most robust association with synaptic transmission [52]. RAB3C encodes a small GTPase. Mollard et al. concluded that similar to RAB3A, RAB3C is also localized on synaptic vesicles and is involved in vesicle trafcking in the nervous system [53]. However, there is no relevant report about the functions of FGF13, RIMBP2, and RAB3C in PD; therefore, further research is needed to understand their contribution to PD. miRNAs are small RNA molecules that can regulate gene expressions after transcription [54]. In our study, eight target miRNAs (hsa-miR-532-5p, hsa-miR-23b-3p, hsa-miR-198, hsa-miR-330-5p, hsa-miR-339-5p, hsa-miR-485-5p, hsa-miR-34a-5p, and hsa-miR-7-5p) were found to be associated with SN in individuals with PD. hsa-miR-532-5p has been shown to be associated with KCNB1. Briggs et al. studied SN in eight idiopathic patients with PD obtained from the Harvard Brain Tissue Resource Center and found a downregulation in hsa-miR-532-5p [55]. A meta-analysis showed that the miRNAs hsa-miR-23b-3p, hsa-miR-7-5p, and hsa-miR-34a-5p were associated with regulating the expression of SNCA [56]. On the other hand, the expression of hsa-miR-198, hsa-miR-330-5p, hsa-miR-339-5p, and hsa-miR-485-5p was signifcantly diferent between SN tissues in eight PD patients and four controls [57]. Consistently, these eight miRNAs were found to be associated with six hub genes (SLC18A2, NR4A2, CACNA1E, DRD2, EN1, and KCNB1) in our study. KCNB1 and CACNA1E are considered the most signifcant genes among these genes above and can be regulated by four miRNAs. Terefore, eight miRNAs might be associated with the pathogenesis of PD. However, future studies should pay more attention to these miRNAs.
Finally, we constructed an mTF-miRNA-gene-gTF regulatory network. We obtained 248 gTFs and 91 mTFs, which can regulate these six hub genes and eight miRNAs. EN1 and NR4A2 were found to be regulated by the highest number of TFs (90 and 61 TFs, respectively). Moreover, TFs, such as HNF4A, FUS, CDX2, SUPT5X, ERG, E2F4, and E2F6, could regulate the most signifcant number of genes and miRNAs in this network. HNF4A, hepatocyte nuclear factor 4 alpha, might infuence gluconeogenesis, diabetes, and lipid homeostasis [58,59]. Previous studies have highlighted the interaction between HNF4A and peroxisome proliferator activator receptor gamma (PPAR-c), a potential therapeutic target in PD [60]. Further, a meta-analysis has identifed HNF4A as the most signifcantly upregulated TF in the blood of patients with PD, while its relative abundance correlated with disease severity in patients with PD [61]. Further studies are required to examine the correlation of FUS, CDX2, SUPT5X, ERG, E2F4, and E2F6 TFs with PD.

Conclusions
We performed bioinformatics analysis on microarray datasets from studies on SN of patients with PD. Our study identifed 86 common DEGs and 14 hub genes in the PPI network. Among them, DRD2, SLC18A2, and SLC6A3 were shown to participate in the pathogenesis of PD by infuencing the function of the dopaminergic synapse. CAC-NA1E, KCNJ6, and KCNB1 might afect the function of the dopaminergic synapse by regulating ion transmembrane transport. Further, we predicted eight miRNAs (hsa-miR-532-5p, hsa-miR-198, hsa-miR-23b-3p, hsa-miR-339-5p, hsa-miR-330-5p, hsa-miR-485-5p, hsa-miR-34a-5p, and hsa-miR-7-5p) which were confrmed to be related to the SN of patients with PD. Using the obtained 248 gTFs and 91 mTFs, we further constructed an mTF-miRNA-gene-gTF regulatory network. Finally, TFs, such as HNF4A, FUS, CDX2, SUPT5X, ERG, E2F4, and E2F6, were found to regulate most number of genes and miRNAs. However, due to the lack of SN samples, additional experiments could not be conducted, and therefore, further studies are required to provide deeper insights into the pathogenesis of PD.

Data Availability
Tis study obtained the gene expression datasets from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). After a thorough review, we chose the GSE49036, GSE7621, GSE8397, and GSE42966. GSE49036 and GSE7621 were based on the GPL570 Platform (Afymetrix Human Genome U133 Plus 2.0 Array). GSE8397 was based on the GPL96 Platform (Afymetrix Human Genome U133A Array) and GPL97 Platform (Afymetrix Human Genome U133B Array), while GSE42966 was based on the GPL4133 Platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F). Te data are freely available online.

Conflicts of Interest
Te authors declare that they have no conficts of interest regarding the publication of this paper.

Authors' Contributions
Yunan Zhou, Zhihui Li, Chunling Chi, and Chunmei Li contributed equally to this work.

Supplementary Materials
Supplementary