Transcriptomic profiling of differentially expressed genes and related pathways in different brain regions of Parkinson’s disease

Background (PD) and many its of how in the transcriptomes We generated a PD mouse model by injecting with MPTP solution. RNA sequencing was performed in the cerebral cortex, hippocampus, striatum, and cerebellum regions of the PD mouse. Compared with the control group, these four brain regions showed significant transcriptomic alterations, with the most differentially expressed genes (DEGs) found in the striatum region. The main DEGs were Lrrk2 , Mtor , Gxylt1 , C920006o11Rik , Vdac1 , Ano3 , Drd4 , and Ncan . DEGs were enriched using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes enrichment analysis methods, which identified significant GO and molecular pathways. In addition, we used network biology methods to analyze protein–protein relationships, which can accelerate the identification of new PD drugs. The results showed that LRRK2, DRD2, IGF-1, GNAI1, GNAI3, PRKACA, PPP2R5C, and PIK3R1 played a major role in protein regulation. Our analysis showed that these DEGs and proteins play an important role in the occurrence and development of PD. Our study also highlighted the potential use of this transcriptomic data for therapeutic strategies and treatment of PD. CC: cortex; HP: hippocampus; ST: striatum; CB: cerebellum; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DEGs: differentially expressed genes; PD: Parkinson's disease; DA: dopamine; LRRK2: leucine-rich repeat kinase 2; PPI: Protein–protein interaction; PCA: Principal component analysis; CN: control; AMPK: AMP-activated protein kinase; PI3K: phosphatidylinositide-3′-OH


5
in public databases [8][9][10]. A study by Alieva et al. (2017) showed significant differences in the expression levels of genes between the ST and substantia nigra in PD [11]. Moni et al. (2019) compared the DEGs in the brain and blood of PD to determine their relationship [12]. Booth et al. showed that leucine-rich repeat kinase 2 (LRRK2) G2019S mutation in the brain can reduce the neuroprotective capacity, thus causing the development of PD pathology [13]. Although PD has been extensively researched, no studies have shown that the disease occurs/develops complete pathways and key differentially expressed genes (DEGs). Therefore, under the current aging worldwide population, it is urgent to study the pathogenesis and control mechanism of PD.
In the present study, we used RNA-Seq to analyze the transcriptome of the cerebral cortex (CC), HP, ST, and CB on the PD. By analyzing these data, we identified the differential gene expression in each brain region, as well as the main functions and pathways for these differential gene enrichments. These data can provide biologically relevant insights into the occurrence/development of PD, and provide strategies for deciphering the potential diagnosis or prognosis of disease pathogenesis [14].

Animals and establishment of PD mouse model
C57Bl/6J mice (8 weeks old) were included in the present study. All mice were purchased from Shanghai Southern Model Biotechnology Co., Ltd. Mice were maintained at 22 °C and allowed free access to food and water in a 12 h dark/light cycle. Eight male mice were randomly divided into a control group and a model group, with four mice in each group. After 3 days of acclimation, the mice in the 6 model group were tested on pole test and rotarod test to collect baseline data.
Then, the mice in the model group were intraperitoneally injected with MPTP (Sigma, St. Louis, MO, USA) solution every day at a dose of 20 mg/kg for 10 days.
The control group mice were injected with normal saline. After 10 days, all the mice in the model group were tested for pole test and rotarod test, and the mice with obvious impairment in their sports ability were screened for subsequent experiments. The study was reviewed and approved by the Ethics Committee of Zhongda Hospital Southeast University.

Sample collection in different brain regions
The animals were anesthetized and killed immediately after the PD mouse model was established. The mice were anesthetized with tribromoethanol (500 mg/kg) (Sigma, Saint Louis, USA), and then the mice were killed by cervical dislocation. The control group and PD group were 3 mice respectively. Brain samples were dissected to isolate the CC, HP, ST, and CB regions. The mice brain was washed with 0.9% precooled saline. The long axis of the brain was placed perpendicular to the operator and the CC region was dissected along the middle of the herringbone. The crescentshaped HP region was clearly seen under the CC region. The left and right brain were gently opened with tweezers. The ST region was inside and close to the side of the CC region, with filamentous stripes, similar to serrations. The CB region was located at the bottom and middle of the brain. The separated mice brain was removed and placed in a 2 mL enzyme-free EP tube. The brain was immediately frozen in liquid nitrogen and stored at −80 °C until subsequent analysis.

Total RNA isolation and cDNA library construction
Total RNA was isolated from the CC, HP, ST, and CB regions using Trizol reagent (Invitrogen, USA). An Agilent 2100 Bioanalyzer (Agilent Technologies, USA) determined the quality and quantity of RNA. Genomic DNA was removed with RNasefree DNaes I (Qiagen, German). Then, the RNA was purified with 1.8× magnetic beads (VAHTATM RNA Clean Beads, N412). NEB offered the NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB #E7490) and the NEBNext rRNA Depletion Kit (NEB #E6310) for the enrichment of non-ribosomal RNA. The cDNA library was constructed using a NEBNext® UltraTM II Directional RNA Library Prep Kit for Illumina®. The first step was to use reverse transcriptase and random primers to convert the interrupted RNA fragments into the first cDNA strand and the second cDNA strand was synthesized with DNA polymerase I and RNase H. The second step was the end repair process and adapter ligation step, and the adapter was then purified. The third step was the polymerase chain reaction (PCR) enrichment of adaptor-ligated DNA. PCR products were purified by NEBNext Sample Purification Beads. The cDNA library concentration was measured using a Qubit® 2.0 fluorometer and each sample was mixed to the same quality. Finally, an Illumina HiSeq 2500 (Illumina, USA) sequencer was used for sequencing.

Analysis of raw data
FasQC software was used to analyze the quality control of the sequencing data and to obtain the sequencing information of all samples. Then, all sample adapters were removed by Trimmomatic V36 software and clean sequence fragments were obtained. The clean reads of each sample were compared with the reference transcriptome of mice for sequence comparison. The quality of sequencing was assessed by analyzing the ratio of clean reads to reference genomes for each sample.

Differentially expressed genes (DEGs) analysis
The read counts of the model group were compared with those of the control 8 group. The gene expression between the model group and the control group was compared and analyzed using the DESeq2 R package (1.16.1). More than twice the differential expression and P < 0.05 were considered as DEGs.

Differentially expression genes enrichment analysis
Gene ontology (GO) enrichment analysis of DEGs was performed from the Metascape website (http://metascape.org/gp/index.html#/main/step1). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs was performed using the David 6.8 website (https://david.ncifcrf.gov/home.jsp). The GO term and KEGG pathways of different brain regions were explored with P < 0.05 as the significant enrichment. The functions of DEGs and major pathways involved were determined by analyzing the significant enrichment of GO term and KEGG pathways.

Protein-protein interaction (PPI) networks
In a PPI network, each protein found in an organism's proteome represents a graph node. If there is any interaction between two nodes, then these two nodes are connected. We used STRING online software (http://string-db.org/) to build a PPI network.

Statistical analysis
The experimental data were analyzed by SPSS 22.0 (SPSS Inc., Chicago, IL, USA) software. One-way ANOVA was used to compare the difference between PD group and control group. Statistical significance was defined as P<0.05.

MPTP-induced PD motor deficits
We evaluated the exercise capacity of MPTP-induced mice by pole test and rotarod test. Fig. 1A, B showed the time required for the movement of the mice in the pole test and rotarod test. In the pole test, MPTP-induced mice were used longer time than the control group (Fig. 1a). Compared with the control group, the MPTP-induced mice had poor coordination and the rotarod test time was significantly shortened ( Fig. 1b). This indicates that MPTP-induced PD mice were successfully constructed.

Summary of RNA-Seq data sets
To understand the transcriptome information of different brain regions in PD mice,  Table S1).
After quality control and filtering, clean reads were obtained. We compared the clean reads to the reference genome with Hisat2. The reads for each gene were calculated using HTSeq-count. Principal component analysis (PCA) was performed based on the DESeq2 method to analyze whether the same brain region samples were together. The results showed that the four brain regions of the control group were divided into different clusters (Fig. 2a). Compared with the CC and HP regions, the transcriptional characteristics of the CB region were significantly different. In addition, the transcriptome characteristics of the ST were significantly different from those of the CC and HP regions. However, the transcriptome characteristics of the CC and HP regions were not significantly different. The four brain regions of the model group were divided into different clusters (Fig. 2b). The difference between the CB region and the other three brain regions was the most obvious.

Differentially expression gene in different brain regions
The transcriptome characteristics of the model group and control group mice were analyzed. DEG analysis explained the difference between the transcription in the 10 model group and control group. According to the values of fold change > 2 and P < 0.05, DEGs were obtained using the DESeq2 R package.
Cluster analysis was performed using the correlation distance and hierarchical algorithm to show gene expression patterns in different brain regions of the control group and model group ( Fig. 3a and b). In the same pathological state, the gene expression in the CB was significantly different from that in the CC, HP, and ST regions.
To analyze the effects of gene expression in different brain regions on PD, the DEGs in the model and control groups in the same brain region were compared. As shown shows the overlap of DEGs in the four brain regions. Of the four brain regions, only two genes (Gxylt1 and C920006O11 Rik) were co-expressed.

Differentially expressed genes in HP and ST regions of Parkinson mice
To further understand the transcript information in the HP and ST regions of mice with PD, RNA sequencing analysis was undertaken. Based on the results, unbiased hierarchical clustering analysis was performed for all genes in the HP region. The data of all genes were visualized using a heat map (Fig. 4a). For data quality control, we performed PCA (Fig. 4b). The results showed that the gene expression in the HP region of the model group was different from that in the HP region of the control group. A total of 179 DEGs were found in the HP region in the model group.

11
Additional file 1: Table S2 lists the top 10 upregulated genes and top 10 downregulated genes with the most significant differences in the HP region. In addition, hierarchical clustering analysis of the DEGs showed that there were differences in gene expression in the ST region of the PD model mice (Fig. 4c).
Additional file 1: Table S3 lists Table S4 and Table   S5 listed the top 10 up-regulated genes and the top 10 down-regulated genes in CC and CB regions, respectively. A heat map of the top 100 genes differentially expressed between the control group and model group in the CB region is provided in Additional file 2: Figure S1.

Gene ontology enrichment analysis of the DEGs in HP and ST regions
The biological functions of DEGs were identified by GO analysis. The GO analysis database can classify DEGs, including "biological processes," "cellular components," and "molecular functions." The GO term corrected by P < 0.05 was defined as significant enrichment of DEGs. In Fig. 5, the significant enrichment of GO terms in the three main categories (cellular component, biological process, and molecular function) is shown. In the HP region, the main five subcategories of biological process were "DNA methylation or demethylation," "posttranscriptional regulation of gene expression," "rhythmic process," "protein localization to organelle," and "regulation of DNA metabolic process," and the two subcategories for molecular function were "guanyl-nucleotide exchange factor activity" and "protein serine/threonine kinase activity" (Fig. 5a). In the ST region, the main five subcategories of cell component were "post synapse," "dendritic tree," "presynapse," "perikaryon", and "endoplasmic reticulum subcompartment." The five subcategories for biological processes were "positive regulation of nervous system development," "synapse organization," "regulation of transmembrane transport," "regulation of exocytosis," and "regulation of intracellular transport"; and five subcategories for molecular function, which are "protein domain specific binding," "protein kinase binding," "protein heterodimerization activity," "GABA receptor binding," and "guanyl-nucleotide binding" (Fig. 5d). Fig. 5b, 5c and 5e, 5f show the top 20 clusters of significantly enriched terms in the HP and ST regions, respectively.

Gene ontology enrichment analysis of the DEGs in CC and CB regions
GO enrichment analysis was performed on DEGs in the CC and CB regions to determine the effects of PD on these two brain regions. The results showed that 10 GO terms were enriched in the CC region (Fig. 6a), in which the cell components were significantly enriched in two GO terms ("axon" and "filopodium"), the biological process was significantly enriched in six GO terms ("negative regulation of protein kinase activity," "regulation of neuron death," "negative regulation of cell cycle," "ribosome biogenesis," "chromatin organization," and "organophosphate biosynthetic process"), and the molecular function was significantly enriched in two GO terms ("myosin binding" and "lipase activity"). A total of 14 GO terms were enriched in the CB region (Fig. 6b), in which the cellular components were significantly enriched in two GO terms, the biological process was significantly enriched in nine GO terms, and the molecular function was significantly enriched in three GO terms. In the CC region, the top five significantly enriched GO terms were 13 "axon," "associative learning," "myosin binding," "negative regulation of protein kinase activity," and "regulation of neuron death" (Fig. 6c). In the CB region, the top five significantly enriched GO terms were "negative regulation of transmembrane transport," "translation factor activity, RNA binding," "histone lysine methylation," "response to glucose," and "dioxygenase activity" (Fig. 6d).

KEGG enrichment analysis
KEGG pathway analysis of DEGs in the HP and ST regions of the model group helped identify the biological pathways related to DEGs. The P < 0.05 pathway was defined as the DEG enrichment pathway. A total of 179 and 521 DEGs in the HP and ST regions were enriched in 122 and 213 KEGG pathways, respectively. Table 1 shows the pathways and related genes that were significantly enriched in the HP region. In the ST region, the top 10 significantly enriched signaling pathways and related genes are shown in Table 2. There were three signaling pathways ("dopaminergic synapse," "glutamatergic synapse," and "adrenergic signaling in cardiomyocytes"), which was consistent with the study by Fu et al. (2019). However, only a few KEGG pathways were enriched in the CC and CB regions. Only the "focal adhesion" and "regulation of actin cytoskeleton" pathways were significantly enriched in the CC region (Table 3), and only the "adrenergic signaling in cardiomyocytes" and "RNA transport" pathways were significantly enriched in the CB region (Table 3).

PPI networks
There were 427 proteins in the PPI network and the nodes represented proteins (Fig.   7a). The most interacting proteins in the PPI network were LRRK2, DA receptor D2 (DRD2), protein kinase A catalytic subunits (PRKACA), PPP2R5C, insulin-like growth factor 1 (IGF-1), PIK3R1, guanine nucleotide binding protein alpha inhibiting 1 (GNAI1), and guanine nucleotide binding protein alpha inhibiting 3 (GNAI3). A total of 15 proteins that may be associated with PD were extracted from the PPI network and the protein network interactions were reconstructed (Fig. 7b).

Discussion
Analysis of DEGs in the different brain regions of PD by RNA-Seq can better explain the complexity of gene function and the biological pathways involved. This is the first comprehensive study that examines the gene expression changes in four brain regions. The present study aimed to investigate the changes of gene expression in pathological pathways in different brain regions of PD mice. DEGs screened from different brain regions can be used as biomarkers or therapeutic targets. Previous studies mainly focus on the transcriptome analysis of ST and substantia nigra [16], but less on HP, CC and CB regions. Emerging data suggest interactions between the dopaminergic systems and the hippocampus in synaptic plasticity, adaptive memory, and motivated behaviour. Therefore, it is important to study the pathological changes in the HP for the study of cognitive dysfunction in PD.
In the present study, the expression of genes in the CB region was significantly different from that in the HP, ST, and CC regions. Previous studies have shown that α-synuclein was found in the CB of PD and its expression level was significantly decreased [17][18][19][20]. Thus, this suggests that regulation of genes and proteins in the CB region may affect Parkinson's disease. Middleton and Strick showed that the CBthalamus-CC pathway affected motor and cognitive functions and was widely connected with the CC via specific pathways [21]. However, the specific mechanism requires further study. Ichinohe et al. also confirmed that there was a specific pathway between the CB and ST in rats, that is, the cerebello-thalamo-motor cortico-striatal pathway and the cerebello-thalamo-striatal pathway affect the function of the ST [22]. Therefore, it is necessary to study the differentially expressed genes in CC, HP, ST and CB brain region and analyze their biological functions.
In the present study, the main DEGs were Lrrk2, Mtor, Gxylt1, C920006o11Rik, Vdac1, Ano3, Drd4, and Ncan. Several genes have been linked to PD. However, Gxylt1 and C920006o11Rik as co-expressed genes in four brain regions have not been reported to date. In the present study, the DEGs were mainly related to DNA methylation or demethylation, brain development, and presynaptic and postsynaptic functions. Many studies have reported that these functions are related to the occurrence of PD [23][24][25][26][27], which is consistent with our results. Wen et al. showed that the brain and blood methylation levels of PD patients significantly changed compared with the control group [23]. This was mainly due to the change in the methylation level caused by TET1 gene mutation, which led to PD [24]. In the nervous system, DNA methylation plays a key role in regulating the cognitive function of the HP. However, the results of the present study showed that the changes in the DNA methylation status were caused by the changes in gene expression of Ctcf, Meg3, Usp9x, Apobec3, and Alkbh1. In addition, DEGs are also involved in the AMP-activated protein kinase (AMPK) signaling pathway, phosphatidylinositide-3′-OH kinase (PI3K)-Akt signaling pathway, mammalian target of rapamycin (mTOR) signaling pathway, gamma-amino butyric acid positive (GABAergic) synapse, glutamatergic synapse, and dopaminergic synapse. AMPK, PI3K-Akt, and mTOR signaling pathways are mutually regulated [28]. AMPK can be activated under oxidative stress, thereby inducing neuronal apoptosis and inhibiting the mTOR pathway, which plays a key role in neurodegenerative diseases [29][30][31][32]. The PI3K-Akt signaling pathway plays an important role in the development and function 16 of neurons, and is involved in the pathogenesis of PD, and directly affects the DA content in the substantia nigra [33]. Activated PI3K-Akt may positively regulate mTOR. The mTOR signaling pathway plays a regulatory role in the growth and development of neuronal cells, including the extension of axons and formation of dendrites [34,35]. Recent studies have shown that abnormal activity of Akt or abnormal expression of mTOR was associated with the occurrence and development of PD [36][37][38][39]. This is consistent with our results, i.e., that the expression of Akt and Mtor genes is significantly downregulated. Therefore, our RNA sequence results are reliable. Proper regulation of AMPK/PI3K-Akt/mTOR signaling may be a potential strategy for the prevention and treatment of PD. On the other hand, glutamate and GABA-mediated excitatory synaptic transmission and inhibitory synaptic transmission were the basis for maintaining the excitation/inhibition balance of the nervous system. Previous studies have shown that the excitation/inhibition imbalance is a key factor in neurological impairment [40,41], such as PD [42]. The D4 receptors in the cerebral cortex and hippocampal neurons have been shown to regulate GABA transmission [43]. However, in the present study, Drd4 gene expression was significantly increased in the cerebellum. This suggests that the cerebellum, cerebral cortex and hippocampus may affect PD by co-regulating the expression of genes. Therefore, this result once again shows that the differentially expressed genes in different brain regions may cause Parkinson's disease. In summary, the interaction of glutamatergic synapse and Glutamatergic synapse plays an important role in the pathogenesis of PD. Therefore, glutamate and GABA may help to develop more effective neuroprotective therapies and provide potential treatment strategies for PD patients.
To determine the molecular network regulators that cause PD, we used STRING software to analyze the correlation between the functions of proteins encoded by DEGs. We used the known PD-associated genes as nodes to reconstruct the PPI network map based on DEGs and related proteins reported in the literature. In the present study, the interactions between 427 DEGs encoded proteins were analyzed.  [45] and is involved in the regulation of receptor transport [46,47]. Beccano-Kelly et al. have shown that the overexpression of LRRK2 can increase the level of DRD2 protein, which played a major role in dopamine receptors [48]. Therefore, the pathogenesis of PD may be caused by alteration in synaptic vesicle trafficking. According to the literature, BAG5 protein has been confirmed to be associated with HSP70, which had been demonstrated to be present in the Lewy body and acts to prevent protein aggregation and neurodegenerative diseases in cells [49,50]. BAG5 protein can interact with a series of proteins to enhance the function of PARKIN to isolate protein aggregation and alleviate PARKIN protective protease function. In combination with these functions, BAG5 is a negative regulator of PARKIN and HSP70, which plays a role in promoting the death of dopaminergic neurons in PD.
Previous studies have been demonstrated that the COR domain of LRRK2 interacted with PARKIN in HEK293T cells, possibly mediated by a third party, and BAG5 may be one of them [51].
In addition, we found that the interaction of DRD2 and IGF-1 was mediated by GNAI1 in the present study. As a regulator of adenylate cyclase, GNAI1 plays an important role in the cholinergic synaptic pathway. Cholinergic M receptor agonists can induce DA release, while DA regulates acetylcholine release via DA receptors [52]. Alcantara et al. found that DA receptors were directly located in striatal cholinergic neurons, and activated DRD2 inhibited the release of acetylcholine [53]. Ebert et al. showed that the overexpression of IGF-1 was associated with the proliferation of dopaminergic neurons in the rat PD model, which may be due to the interaction of IGF-1 with DA receptors via cholinergic receptors, thus regulating the transport of DA neurotransmitters [54]. From the transcriptome data of the present study, it is speculated that the Prkaca, Ppp2r5c, Pik3r1, and Adcy1 genes may be related to the imbalance of DA and acetylcholine; however, the specific functions remain unclear and require further study and verification. A thorough understanding of these node proteins will contribute to the study of the pathogenesis of PD and the development of PD therapeutic drugs.

Ethics approval and consent to participate
All animal research protocols for this work were reviewed and approved by the Ethics Committee of Zhongda Hospital Southeast University.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Consent for publication
Not applicable.  Figure 1 Effect of PD on Pole test and Rotarod test. The network of differential expression in gene protein interaction.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.