Transcriptome sequencing in a 6-hydroxydopamine rat model of Parkinson’s disease

ABSTRACT

The neurotoxin 6-hydroxydopamine (6-OHDA) has been widely exploited as a tool for modeling Parkinson’s disease (PD) in the rat. This study aimed to provide a comprehensive profile of the mRNAs and long noncoding RNAs (lncRNAs) in rats treated with 6-OHDA as a model of PD. Female SPF Wistar rats were randomly divided into two groups: a PD model group and a control group. The PD model was induced by 6-OHDA injection. RNA-seq analysis was performed on 6-OHDA-treated rats and corresponding controls. Novel lncRNAs were identified. Differentially expressed genes (DEGs) and differentially expressed lncRNAs were identified in the PD group compared with controls. Gene Ontology function and pathway enrichment analyses were conducted on the DEGs, followed by construction of a protein–protein interaction (PPI) network. In addition, prediction of lncRNA target genes and function prediction of lncRNAs were performed. Moreover, microRNAs (miRNAs) that interacted with the DEGs and differentially expressed lncRNAs were predicted to construct a miRNA–lncRNA–mRNA regulatory network. A total of 536 DEGs and 512 differentially expressed lncRNAs (44 up-regulated and 10 down-regulated known lncRNAs; 407 up-regulated and 51 down-regulated novel lncRNAs) were identified in the PD rat model compared with controls. The DEGs and target genes of lncRNAs were mainly associated with the innate immune response, 2′-5′-oligoadenylate synthetase activity, GTPase activity, GTP binding and the RIG-I-like receptor signaling pathway. IRF7 and ISG15 were hub proteins in the PPI network. Many mRNAs and lncRNAs interacted with other molecules in a competing endogenous RNA network, such as MAS1, TMPRSS2, NPTX1, XLOC_016191, XLOC_026924 and XLOC_005439. We conclude that IRF7, ISG15, MAS1, TMPRSS2, NPTX1, XLOC_016191, XLOC_026924 and XLOC_005439 may contribute critical roles in the pathogenesis of PD.

INTRODUCTION

Parkinson’s disease (PD) is a common neurodegenerative disease that impairs voluntary movement in middle-aged and elderly people, affecting 1–2% of the population over the age of 60 (Jungverdorben et al., 2017). It is characterized by dopamine deficiency in areas of the midbrain, resulting in motor dysfunctions such as slowness of movement, akinesia, rigidity, postural instability and resting tremor (Vo et al., 2017). There is still no method that can effectively prevent the occurrence and development of PD. Therefore, identifying the molecular mechanisms of the etiology and progression of the disease is a current research hotspot.

In recent years, an increasing number of studies have identified genes as likely candidates for involvement in PD. For instance, a network-based meta-analysis demonstrated that HNF4A and PTBP1 are longitudinally dynamic biomarkers for PD (Santiago and Potashkin, 2015). Yeo et al. showed that decreased expression of serum- and glucocorticoid-dependent kinase 1 promotes the expression of alpha-synuclein, which is related to dopaminergic cell death in chronic Parkinsonism mice (Yeo et al., 2018). Long noncoding RNAs (lncRNAs) are defined as transcripts longer than 200 nucleotides that are transcribed from DNA but not translated into proteins (Kopp and Mendell, 2018). Soreq et al. used an all-RNA transcriptome sequencing method to screen lncRNAs associated with PD and revealed, for example, a significant increase in the expression of lnc-frg1-3, which is associated with muscle rigidity in PD (Soreq et al., 2014). Another study found that the expression of a ubiquitin carboxyl-terminal hydrolase Uchl1 antisense strand lncRNA was regulated by the major transcription factor Nurr1, which is involved in the differentiation and maintenance of dopamine cells (Carrieri et al., 2015). However, there are substantially more genes and lncRNAs related to PD remaining to be explored.

In the current study, RNA-seq analysis was performed on a 6-hydroxydopamine (6-OHDA) rat model of PD and corresponding control animals. Novel lncRNAs were identified. Differentially expressed mRNAs and lncRNAs were identified in the PD group compared with controls. Gene Ontology (GO) function analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted on differentially expressed genes (DEGs), followed by construction of a protein–protein interaction (PPI) network. In addition, lncRNA target genes and function of lncRNAs were predicted. Moreover, microRNAs (miRNAs) that interacted with the DEGs and differentially expressed lncRNAs were predicted before the construction of a miRNA–lncRNA–mRNA regulatory network. This study aimed to provide a comprehensive profile of the mRNAs and lncRNAs in a 6-OHDA rat model of PD and to identify important new gene targets and lncRNA targets for PD, which should contribute to our understanding of the pathogenesis of PD.

MATERIALS AND METHODS

Animals

Eighteen inbred female SPF Wistar rats weighing 200–220 g, 10 to 12 weeks of age, were used in the present study. Before experiments, the animals were kept in a controlled facility (temperature, 24 ± 2 ℃; humidity, 50 ± 5%) with ad libitum access to food and water under a 12/12 h light/dark cycle for one week. All animal procedures were carried out in compliance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals.

Experimental design and surgical procedure

The 18 experimental animals were randomly divided into two groups: control group and PD model group (n = 9 for each group). Rats were anesthetized with an intraperitoneal injection of 3% chloral hydrate (10 ml/kg). The animals were fixed into a stereotaxic instrument (Zhongshi Dichuang, Beijing) secured by ear bars adapted specifically for each rat. The top of the head was shaved. The skull was exposed by an incision in the midline and small holes were drilled in the skull to allow a 4-μl injection of 6-OHDA (Sigma, St. Louis) (dissolved in 0.2% ascorbic acid saline solution) into each of the substantia nigra (antero-posterior: +4.8 ± 0.1 mm, medio-lateral: −1.7 ± 0.1 mm, dorso-ventral: −7.5 ± 0.1 mm below dural surface) and the ventrolateral region of the caudate-putamen (antero-posterior: 4.6 ± 0.1 mm, medio-lateral: 0.9 ± 0.1 mm, V: 7.3 ± 0.1 mm from bregma) according to Paxinos and Watson (1997). Sham-operated rats, which received 0.02% ascorbic acid in sterile saline bilaterally at the same coordinates, served as controls in all experiments. 6-OHDA was injected at a rate of 1 μl/min, and the injection syringe was left in place for a further 10 min after injection to enable full absorption of the solution. Surgical clips were used for closing the incisions before rats were returned to the animal housing facility for recovery. To avoid infections, the rats received a daily injection of penicillin (50,000 units) for five days afterwards.

Rotational asymmetry analysis

At one week after surgery, behavioral changes and abnormal behavior of the 6-OHDA-lesioned PD rats were observed. Rats were intraperitoneally injected with 0.5 mg/kg of apomorphine (Sigma), and the number of rotations of rats thus induced was recorded manually for a 10-min period. Further inductions were performed at two, four and six weeks after surgery. If the rat turned steadily to the left, and the number of rotations was > 7/min in all of these measurements, it was considered a valid PD model and was used for further experiments (da Conceição et al., 2010).

RNA isolation and sequencing

Rats were killed by carbon dioxide asphyxiation and the brains were rapidly removed. The right and left striatum were isolated and rinsed with saline, snap-frozen in liquid nitrogen and kept at −80 ℃ until further use. Total RNA from the rat brain was isolated with RNAiso Plus reagent (Takara Bio, Kusatsu) according to the manufacturer’s instructions. RNA-seq libraries were prepared, followed by Illumina sequencing by Sangon Biotech (Shanghai).

Sequence read processing and genome mapping

The raw sequencing data were first filtered to generate clean data using SeqPrep (Schmieder and Edwards, 2011) (https://github.com/jstjohn/SeqPrep). Reads containing adaptors were removed. Bases having a quality score of less than 20 at the 3′ end of the read were trimmed; if the remaining sequence still had a quality score of less than 10, the entire read was removed. Reads containing more than 10% of unknown nucleotides and reads with a minimum length of less than 30 bp were discarded. Clean reads were mapped to the rat reference genome (Rnor_6.0.88, http://www.ensembl.org/index.html) using tophat2 (Schmieder and Edwards, 2011) software (v2.1.0, http://ccb.jhu.edu/software/tophat/index.shtml) using the following parameters: read-mismatches 2, read-edit-dist 2, max-intron-length 5000000, num-threads 8, mate-inner-dist 40.

RNA sequence analysis and identification of differentially expressed transcripts

Transcript assembly was performed using Cufflinks (Mortazavi et al., 2008) software (v2.2.1, http://cufflinks.cbcb.umd.edu/) and the number of read-counts of the corresponding genes was then obtained using the HTSeq tool (Mortazavi et al., 2008) (V0.6.1p2) based on the rat genome annotation information. The results were normalized by the CPM (counts per million) method, and genes having a CPM value of less than 0.1 in five or more samples were defined as genes with low expression abundance. The obtained results were classified into lncRNA and mRNA according to the gene type in the annotation information.

A transcript of class code type iux in the assembled transcript (gene) is considered as an unknown transcript (i indicates that the predicted transcript is located entirely within the intron range, u indicates an unknown intergenic transcript, and x indicates exonic overlap with a reference on the opposite strand). Three software packages were used to predict the coding ability of transcript sequences, CPC (Kong et al., 2007) (v0.9-r2, http://cpc2.cbi.pku.edu.cn/), CNCI (Sun et al., 2013) (https://github.com/www-bioinfo-org/CNCI) and PfamScan (Mistry et al., 2007) (https://www.ebi.ac.uk/Tools/pfa/pfamscan/), with the following threshold settings: CPC_threshold = 0 (assembled transcripts > 0 was mRNA, or assembled transcripts < 0 was lncRNA), CNCI_threshold = 0 (assembled transcripts > 0 was mRNA, or assembled transcripts < 0 was lncRNA) and Pfam as a protein database. The results obtained from the intersection of t for all three judgment methods were used for subsequent analysis.

The edgeR (Lun et al., 2016) software package in R was used for differential expression analysis. Differentially expressed lncRNAs and mRNAs were identified in the PD model group compared with controls. Prior to analysis, genes with lower expression abundance were deleted and thresholds were set as |log₂ fold change (FC)| > 0.585 and P < 0.05. In addition, a two-dimensional clustering heatmap was generated.

GO and KEGG enrichment analysis

GO functions and KEGG pathways of DEGs were analyzed using the DAVID (Huang et al., 2008) (https://david.ncifcrf.gov/) tool. GO terms and pathways with gene count ≥ 2 were enriched. P values < 0.05 based on the hypergeometric test were considered significant.

PPI prediction

PPI analysis was performed on DEGs using the STRING (Huang et al., 2008) (v10.5, https://string-db.org/) online tool. Medium Confidence (combined score) > 0.4 was selected as the threshold. Cytoscape (Shannon et al., 2003) (v3.5.1, http://www.cytoscape.org/) was used for the construction of a PPI network. The topology of the network was analyzed, the connectivity of network nodes was ranked, and nodes with higher connectivity were considered as the important nodes.

Prediction of lncRNA target genes and function prediction of lncRNAs

For the resulting DEGs and differentially expressed lncRNAs, the correlation between them was calculated using Pearson’s correlation coefficient (Benesty et al., 2009), and the P value was corrected by the Benjamini and Hochberg method (Benjamini and Hochberg, 1995). A correlation coefficient greater than 0.9 and corrected P < 0.01 were chosen as the threshold. mRNAs that had a correlation with lncRNAs were potential target genes for lncRNAs. The lncRNA–mRNA network was built using cytoscape software (Shannon et al., 2003).

The functions of lncRNAs were predicted by performing GO function and KEGG pathway enrichment analysis on lncRNA target genes. We used the clusterProfiler (Yu et al., 2012) (http://www.bioconductor.org/) package to compare the KEGG pathways of lncRNA target genes. The screening threshold was set as P < 0.05.

miRNA prediction and ceRNA network construction

Competing endogenous RNA (ceRNA) regulatory networks are more elaborate and complex, involving more RNA molecules, including mRNA, pseudogene transcripts, lncRNAs and miRNAs (Yuan et al., 2017). We downloaded conserved miRNA-regulated mRNA data from miRanda (http://34.236.212.39/microrna/home.do), and the miRanda tool (Betel et al., 2010) was also applied to predict binding sites between miRNAs and lncRNAs (parameter: -sc 160, -en -20). miRNA–lncRNA, miRNA–mRNA and lncRNA–mRNA relationships were integrated to obtain the miRNA–lncRNA–mRNA regulatory network.

RESULTS

Overview of RNA-seq and mRNA and lncRNA identification

After the rotational behavior test, six PD rats were valid and were used for the subsequent sequencing. The six PD rats in the model group were named MW1-L, MW1-R, MW2-L, MW2-R, MW3-L and MW3-R; six rats in the control group were chosen for sequence analysis and named CW7-L, CW7-R, CW8-L, CW8-R, CW9-L and CW9-R. After sequence data preprocessing, more than 91% reads in each sample were aligned to the reference genome (Supplementary Table S1). The reads were mapped to the reference genome, resulting in 40,752 assembled transcripts, of which 13,613 genes had low expression abundance.

Under the threshold conditions (|log₂ FC| > 0.585 and P < 0.05), a total of 536 differentially expressed mRNAs (418 up-regulated and 118 down-regulated) and 512 differentially expressed lncRNAs (44 up-regulated known lncRNAs, 10 down-regulated known lncRNAs, 407 up-regulated novel lncRNAs and 51 down-regulated novel lncRNAs) were identified. Detailed information about differentially expressed mRNAs and known lncRNAs is presented in Supplementary Table S2. Figure 1 shows two-dimensional heatmaps of differentially expressed mRNAs and lncRNAs.

Fig. 1.

Heatmaps and functional enrichment analysis. (A) Heatmap of differentially expressed genes (DEGs). Red indicates up-regulated genes and blue indicates down-regulated genes. (B) Heatmap of differentially expressed lncRNAs. Red indicates up-regulated lncRNAs and blue indicates down-regulated lncRNAs. (C) The top five GO terms and pathways enriched in DEGs. Number of genes: the number of differentially expressed genes enriched in the indicated terms. The gray trend lines represent the –log₁₀ (P values), and the values (0–30) on the horizontal axis also apply to these data. BP, biological process; CC, cellular component; MF, molecular function.

GO and KEGG pathway enrichment analysis of DEGs

A total of 141 GO terms and 29 KEGG pathways were identified for DEGs. In the GO analysis, the top five enriched terms in the biological process (BP), cellular component (CC) and molecular function (MF) categories are shown in Fig. 1C. We found that the DEGs were mainly associated with the following BP terms: defense response to virus, cellular response to interferon-beta, innate immune response and negative regulation of viral genome replication, whereas with regard to the CC terms the DEGs appeared to be active in both the extracellular region and the extracellular matrix (ECM) (Fig. 1C). The main MF terms were 2′-5′-oligoadenylate synthetase activity, double-stranded RNA binding, CXCR3 chemokine receptor binding, GTPase activity and GTP binding. Additionally, the DEGs were mainly involved in pathways such as protein digestion and absorption, herpes simplex and influenza A infection, RIG-I-like receptor signaling pathway and ECM–receptor interaction.

PPI analysis

A PPI network for DEGs was constructed (Fig. 2), comprising 149 nodes and 434 interaction pairs. The key nodes in the network were identified. Table 1 lists the top 20 node genes with the highest connectivity in the network, such as interferon regulatory factor 7 (Irf7, degrees = 18), ISG15 ubiquitin-like modifier (Isg15, degrees = 13), ubiquitin-like modifier activating enzyme 7 (Uba7, degrees = 13), 2′-5′-oligoadenylate synthetase 1 (Oas1a, degrees = 12) and C-X-C motif chemokine ligand 10 (Cxcl10, degrees = 11).

Fig. 2.

Protein–protein interaction (PPI) network. Red and green circular nodes represent up-regulated and down-regulated genes, respectively.

Table 1. The top 20 genes in the protein-protein interaction network

Gene	Change	Degrees	Gene	Change	Degrees
IRF7	up	18	RGD1565117	up	8
ISG15	up	13	CXCL13	up	8
UBA7	up	13	HERC6	up	8
OAS1A	up	12	UBE2L6	up	7
CXCL10	up	11	DTX3L	up	7
OASL	up	10	MX1	up	7
USP18	up	9	NPY2R	up	7
RPS4X	down	9	NPY4R	up	7
CXCL11	up	8	CCR1L1	down	7
CXCL9	up	8	PTGDR2	up	7

Prediction of target genes for differentially expressed lncRNAs and functional analysis

Using the CPC, CNCI and PfamScan software to obtain intersection results, 7,506 novel lncRNAs were obtained (Fig. 3A), 458 of which were differentially expressed. Target genes for differentially expressed lncRNAs were predicted using Pearson’s correlation coefficient. lncRNAs with more than five target genes were analyzed. A lncRNA–mRNA interaction network was constructed (Fig. 3B). The top 20 lncRNAs with the highest connectivity are shown in Table 2, including XLOC_026924, XLOC_029397, XLOC_004631, XLOC_005439 and XLOC_028318.

Fig. 3.

Prediction of target genes for differentially expressed lncRNA and functional analysis. (A) Venn diagram of lncRNA prediction using three software packages. (B) lncRNA–mRNA interaction network. Red represents up-regulation, green represents down-regulation, circular nodes represent mRNA and triangular nodes represent lncRNA. (C) The top five GO terms and pathways enriched for target genes of lncRNAs.

Table 2. The top 20 lncRNAs with highest connectivity in the lncRNA-mRNA interaction network

Name	Change	Degrees	Name	Change	Degrees
XLOC_026924	up	41	XLOC_017775	down	15
XLOC_029397	up	40	XLOC_018598	down	15
XLOC_004631	up	34	XLOC_022926	up	15
XLOC_005439	up	26	Rn50_5_1638.1	down	14
XLOC_028318	down	19	XLOC_006399	down	13
XLOC_018657	up	18	AABR07027137.1	down	11
XLOC_037769	down	18	AABR07029901.1	up	10
XLOC_029657	down	17	XLOC_025867	up	10
XLOC_010572	down	16	XLOC_016202	up	9
XLOC_016191	up	15	XLOC_001547	down	8

To further understand the characteristics of these lncRNAs, functional enrichment analysis was performed on the target genes involved in the network. A total of 46, 15 and seven terms were enriched for the GO categories BP, MF and CC, respectively, and 10 KEGG pathways were enriched simultaneously. The top five most significant terms of BP, MF and CC as well as pathways are represented in Fig. 3C. The lncRNAs were associated with the defense response to virus, the cellular response to interferon-beta, the innate immune response, 2′-5′-oligoadenylate synthetase activity, GTPase activity, GTP binding, calcium channel activity and the RIG-I-like receptor signaling pathway.

miRNA prediction and ceRNA construction

We downloaded the conserved miRNAs for rat mRNA regulation from the miRanda website and screened for miRNAs that regulate the target genes of differentially expressed lncRNAs. In addition, miRanda software was used to predict binding sites between miRNAs and lncRNAs. Thus, the top 10 miRNAs that had binding sites with the differentially expressed lncRNAs are listed in Table 3, and are rno-miR-466b, rno-miR-349, rno-miR-34b, rno-miR-200b, rno-miR-483, rno-miR-27a, rno-miR-532-5p, rno-miR-34c, rno-miR-532-5p and rno-miR-466b. To integrate the miRNA–lncRNA, miRNA–mRNA and lncRNA–mRNA interactions, a miRNA–lncRNA–mRNA regulatory network was constructed (Fig. 4). Several clusters were extracted from the ceRNA network (Fig. 5). From the results, we found that DEGs such as neuronal pentraxin 1 (Nptx1), transient receptor potential cation channel subfamily V member 5 (Trpv5), MAS1 proto-oncogene, G protein-coupled receptor (Mgs1), MX dynamin like GTPase 1 (Mx1), calcium-binding protein 7 (Cabp7), transmembrane protease, serine 2 (Tmprss2), and radical S-adenosyl methionine domain containing 2 (Rsad2) could be regulated by many lncRNAs and miRNAs.

Table 3. The top 10 miRNAs that had binding sites with differentially expressed lncRNAs

miRNA	lncRNA	Max. score	Max. energy
rno-miR-466b	XLOC_022926	183.00	−26.10
rno-miR-349	XLOC_004631	183.00	−33.06
rno-miR-34b	Cd93	181.00	−26.90
rno-miR-200b	XLOC_036577	181.00	−20.54
rno-miR-483	XLOC_029657	178.00	−32.47
rno-miR-27a	Cd93	178.00	−24.35
rno-miR-532-5p	XLOC_029397	177.00	−28.28
rno-miR-34c	Cd93	177.00	−28.30
rno-miR-532-5p	XLOC_004631	176.00	−29.62
rno-miR-466b	AABR07006685.1	176.00	−29.69

Fig. 4.

miRNA–lncRNA–mRNA regulatory network. Red represents up-regulation; triangles are lncRNAs, diamonds are miRNAs and circles are mRNAs.

Fig. 5.

Several clusters extracted from the ceRNA network. Red circles represent up-regulated genes, red triangles represent up-regulated lncRNAs and diamonds represent miRNAs. Blue lines represent the miRNA–mRNA relationship, green lines represent the miRNA–lncRNA relationship and dashed lines represent the lncRNA–mRNA co-expression relationship.

DISCUSSION

In the current study, 536 DEGs (418 up-regulated and 118 down-regulated) and 512 differentially expressed lncRNAs (44 up-regulated known lncRNAs, 10 down-regulated known lncRNAs, 407 up-regulated novel lncRNAs and 51 down-regulated novel lncRNAs) were identified in a PD rat model compared with controls. The DEGs and target genes of lncRNAs were mainly associated with the innate immune response, 2′-5′-oligoadenylate synthetase activity, GTPase activity, GTP binding and the RIG-I-like receptor signaling pathway. IRF7 and ISG15 were hub proteins in the PPI network. Many mRNAs and lncRNAs interacted with other molecules in the ceRNA network such as Mas1, Tmprss2, Nptx1, XLOC_016191, XLOC_026924 and XLOC_005439.

Chronic neuroinflammation processes in PD have an important role that leads to protracted neuronal degeneration (Lee et al., 2009; Tansey and Goldberg, 2010). Evidence indicates that 2′-5′-oligoadenylate synthetase-like protein can be a useful biomarker for neurodegenerative disease diagnosis (Goldknopf et al., 2006). Several studies demonstrated the important roles of GTPase activity and GTP binding in neuronal toxicity and the development of PD (Guo et al., 2007; Ito et al., 2007; West et al., 2007). In addition, one study reported that cytosolic RIG-I-like helicases could be negative regulators of inflammation in the central nervous system (Dann et al., 2011). These studies provide an explanation for the results obtained by functional enrichment analysis, suggesting that our analysis is reliable.

Evidence demonstrates that type I interferon (IFN)-dependent immune responses are under the control of IRF7 (Honda et al., 2005). Recently, type I IFNs were reported to contribute to the early neuroinflammatory response and neuronal cell death in the mouse model of PD (Main et al., 2016). In addition, a previous genome-wide association study indicated that IRF7 was associated with systemic lupus erythematosus, which suggests a mechanistic relationship between systemic lupus erythematosus and PD via IFN signaling (Liu et al., 2015). The ubiquitin-like protein ISG15 is one of several hundred proteins induced by IFN-α/-β (Zhao et al., 2005). ISG15-mediated impairment of protein degradation was shown to contribute to progressive neurodegeneration (Wood et al., 2011). In this study, IRF7 and ISG15 were distinct hub proteins in the PPI network. Taken together, our results suggest that IRF7 and ISG15 play significant roles in the pathogenesis of PD, which warrants further experimental validation in future.

Neuronal pentraxin II (NPTX2) is highly up-regulated in PD (Moran et al., 2008). Qu and D’Mello showed that NPTX1 can be a target of histone deacetylase-3-mediated neurodegeneration (Qu and D’Mello, 2018). In this study, Nptx1 was an important player in the ceRNA network. Additionally, Mas1 and Tmprss2 were also shown to be critical in the ceRNA network. The appearance of Mx1 protein in reactive microglia may contribute to Alzheimer’s disease pathology (Wennström et al., 2015). Lahiry et al. showed a mutation in the serine protease TMPRSS4 in a novel neurodegenerative disorder (Lahiry et al., 2013). However, no study has reported the roles of XLOC_016191, XLOC_026924 or XLOC_005439. Nevertheless, we suggest that MAS1, TMPRSS2, NPTX1, XLOC_016191, XLOC_026924 and XLOC_005439 play roles in the development of PD. Future studies are required to address their functional significance.

CONCLUSIONS

In conclusion, the present study has provided a comprehensive analysis of the DEGs and differentially expressed lncRNAs in a rat PD model by RNA-seq and bioinformatics analysis. This study indicated significant functions for the immune response, 2′-5′-oligoadenylate synthetase activity, GTPase activity, GTP binding, and the RIG-I-like receptor signaling pathway. Moreover, IRF7, ISG15, MAS1, TMPRSS2, NPTX1, XLOC_016191, XLOC_026924 and XLOC_005439 may contribute critical roles in the pathogenesis of PD, which provides a foundation for further research into the molecular mechanisms of PD and may aid the identification of sensitive biomarkers and therapeutic targets for this disease.

ACKNOWLEDGMENTS

This work was supported by the Jilin Province Neurological Diseases Precision Medicine Science and Technology Innovation Center (No. 20170623006TC).

REFERENCES

•  Benesty,  J.,  Chen,  J.,  Huang,  Y., and  Cohen,  I. (2009) Pearson correlation coefficient. In: Noise Reduction in Speech Processing, vol. 2. pp. 37–40. Springer, Heidelberg.
•  Benjamini,  Y., and  Hochberg,  Y. (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statist. Soc. B (Methodological), 57, 289–300.
•  Betel,  D.,  Koppal,  A.,  Agius,  P.,  Sander,  C., and  Leslie,  C. (2010) Comprehensive modeling of microRNA targets predicts functional non-conserved and non-canonical sites. Genome Biol. 11, R90.
•  Carrieri,  C.,  Forrest,  A. R. R.,  Santoro,  C.,  Persichetti,  F.,  Carninci,  P.,  Zucchelli,  S., and  Gustincich,  S. (2015) Expression analysis of the long non-coding RNA antisense to Uchl1 (AS Uchl1) during dopaminergic cells’ differentiation in vitro and in neurochemical models of Parkinson’s disease. Front. Cell. Neurosci. 9, 114.
•  da Conceição,  F. S. L.,  Ngo-Abdalla,  S.,  Houzel,  J.-C., and  Rehen,  S. K. (2010) Murine model for Parkinson’s disease: from 6-OH dopamine lesion to behavioral test. J. Vis. Exp. 15, 1376.
•  Dann,  A.,  Poeck,  H.,  Croxford,  A. L.,  Gaupp,  S.,  Kierdorf,  K.,  Knust,  M.,  Pfeifer,  D.,  Maihoefer,  C.,  Endres,  S.,  Kalinke,  U., et al. (2011) Cytosolic RIG-I-like helicases act as negative regulators of sterile inflammation in the CNS. Nat. Neurosci. 15, 98–106.
•  Goldknopf,  I. L.,  Sheta,  E. A.,  Appel,  S. H.,  Simpson,  E. P., and  Yen,  A. A. (2006) 2’-5’-oligoadenylate synthetase like protein as a biomarker for neurodegenerative disease. U.S. Patent US20060115856A1. 2006-06-01.
•  Guo,  L.,  Gandhi,  P. N.,  Wang,  W.,  Petersen,  R. B.,  Wilson-Delfosse,  A. L., and  Chen,  S. G. (2007) The Parkinson’s disease-associated protein, leucine-rich repeat kinase 2 (LRRK2), is an authentic GTPase that stimulates kinase activity. Exp. Cell Res. 313, 3658–3670.
•  Honda,  K.,  Yanai,  H.,  Negishi,  H.,  Asagiri,  M.,  Sato,  M.,  Mizutani,  T.,  Shimada,  N.,  Ohba,  Y.,  Takaoka,  A.,  Yoshida,  N., et al. (2005) IRF-7 is the master regulator of type-I interferon-dependent immune responses. Nature 434, 772–777.
•  Huang,  D. W.,  Sherman,  B. T., and  Lempicki,  R. A. (2008) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57.
•  Ito,  G.,  Okai,  T.,  Fujino,  G.,  Takeda,  K.,  Ichijo,  H.,  Katada,  T., and  Iwatsubo,  T. (2007) GTP binding is essential to the protein kinase activity of LRRK2, a causative gene product for familial Parkinson’s disease. Biochemistry 46, 1380–1388.
•  Jungverdorben,  J.,  Till,  A., and  Brüstle,  O. (2017) Induced pluripotent stem cell-based modeling of neurodegenerative diseases: a focus on autophagy. J. Mol. Med. (Berl) 95, 705–718.
•  Kong,  L.,  Zhang,  Y.,  Ye,  Z.-Q.,  Liu,  X.-Q.,  Zhao,  S.-Q.,  Wei,  L., and  Gao,  G. (2007) CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349.
•  Kopp,  F., and  Mendell,  J. T. (2018) Functional classification and experimental dissection of long noncoding RNAs. Cell 172, 393–407.
•  Lahiry,  P.,  Racacho,  L.,  Wang,  J.,  Robinson,  J. F.,  Gloor,  G. B.,  Rupar,  C. A.,  Siu,  V. M.,  Bulman,  D. E., and  Hegele,  R. A. (2013) A mutation in the serine protease TMPRSS4 in a novel pediatric neurodegenerative disorder. Orphanet J. Rare Dis. 8, 126.
•  Lee,  J. K.,  Tran,  T., and  Tansey,  M. G. (2009) Neuroinflammation in Parkinson’s disease. J. Neuroimmune Pharmacol. 4, 419–429.
•  Liu,  F. C.,  Huang,  W. Y.,  Lin,  T. Y.,  Shen,  C. H.,  Chou,  Y. C.,  Lin,  C. L.,  Lin,  K. T., and  Kao,  C. H. (2015) Inverse sssociation of Parkinson disease with systemic lupus erythematosus: a nationwide population-based study. Medicine (Baltimore) 94, e2097.
•  Lun,  A. T. L.,  Chen,  Y., and  Smyth,  G. K. (2016) It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods Mol. Biol. 1418, 391–416.
•  Main,  B. S.,  Zhang,  M.,  Brody,  K. M.,  Ayton,  S.,  Frugier,  T.,  Steer,  D.,  Finkelstein,  D.,  Crack,  P. J., and  Taylor,  J. M. (2016) Type-1 interferons contribute to the neuroinflammatory response and disease progression of the MPTP mouse model of Parkinson’s disease. Glia 64, 1590–1604.
•  Mistry,  J.,  Bateman,  A., and  Finn,  R. D. (2007) Predicting active site residue annotations in the Pfam database. BMC Bioinformatics 8, 298.
•  Moran,  L. B.,  Hickey,  L.,  Michael,  G. J.,  Derkacs,  M.,  Christian,  L. M.,  Kalaitzakis,  M. E.,  Pearce,  R. K. B., and  Graeber,  M. B. (2008) Neuronal pentraxin II is highly upregulated in Parkinson’s disease and a novel component of Lewy bodies. Acta Neuropathol. 115, 471–478.
•  Mortazavi,  A.,  Williams,  B. A.,  McCue,  K.,  Schaeffer,  L., and  Wold,  B. (2008) Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628.
•  Paxinos,  G., and  Watson,  C. (1997) The Rat Brain in Stereotaxic Coordinates, 3rd Ed. Academic Press, San Diego, USA.
•  Qu,  Z., and  D’Mello,  S. R. (2018) Proteomic analysis identifies NPTX1 and HIP1R as potential targets of histone deacetylase-3-mediated neurodegeneration. Exp. Biol. Med. (Maywood) 243, 627–638.
•  Santiago,  J. A., and  Potashkin,  J. A. (2015) Network-based metaanalysis identifies HNF4A and PTBP1 as longitudinally dynamic biomarkers for Parkinson’s disease. Proc. Natl. Acad. Sci. USA 112, 2257–2262.
•  Schmieder,  R., and  Edwards,  R. (2011) Quality control and preprocessing of metagenomic datasets. Bioinformatics 27, 863–864.
•  Shannon,  P.,  Markiel,  A.,  Ozier,  O.,  Baliga,  N. S.,  Wang,  J. T.,  Ramage,  D.,  Amin,  N.,  Schwikowski,  B., and  Ideker,  T. (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504.
•  Soreq,  L.,  Guffanti,  A.,  Salomonis,  N.,  Simchovitz,  A.,  Israel,  Z.,  Bergman,  H., and  Soreq,  H. (2014) Long non-coding RNA and alternative splicing modulations in Parkinson’s leukocytes identified by RNA sequencing. PLoS Comput. Biol. 10, e1003517.
•  Sun,  L.,  Luo,  H.,  Bu,  D.,  Zhao,  G.,  Yu,  K.,  Zhang,  C.,  Liu,  Y.,  Chen,  R., and  Zhao,  Y. (2013) Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41, e166.
•  Tansey,  M. G., and  Goldberg,  M. S. (2010) Neuroinflammation in Parkinson’s disease: its role in neuronal death and implications for therapeutic intervention. Neurobiol. Dis. 37, 510–518.
•  Vo,  A.,  Sako,  W.,  Fujita,  K.,  Peng,  S.,  Mattis,  P. J.,  Skidmore,  F. M.,  Ma,  Y.,  Uluğ,  A. M., and  Eidelberg,  D. (2017) Parkinson’s disease-related network topographies characterized with resting state functional MRI. Hum. Brain Mapp. 38, 617–630.
•  Wennström,  M.,  Surova,  Y.,  Hall,  S.,  Nilsson,  C.,  Minthon,  L.,  Hansson,  O., and  Nielsen,  H. M. (2015) The inflammatory marker YKL-40 is elevated in cerebrospinal fluid from patients with Alzheimer’s but Not Parkinson’s disease or dementia with Lewy bodies. PLoS One 10, e0135458.
•  West,  A. B.,  Moore,  D. J.,  Choi,  C.,  Andrabi,  S. A.,  Li,  X.,  Dikeman,  D.,  Biskup,  S.,  Zhang,  Z.,  Lim,  K. L.,  Dawson,  V. L., et al. (2007) Parkinson’s disease-associated mutations in LRRK2 link enhanced GTP-binding and kinase activities to neuronal toxicity. Hum. Mol. Genet. 16, 223–232.
•  Wood,  L. M.,  Sankar,  S.,  Reed,  R. E.,  Haas,  A. L.,  Liu,  L. F.,  McKinnon,  P., and  Desai,  S. D. (2011) A novel role for ATM in regulating proteasome-mediated protein degradation through suppression of the ISG15 conjugation pathway. PLoS One 6, e16422.
•  Yeo,  S.,  Sung,  B.,  Hong,  Y. M.,  van den Noort,  M.,  Bosch,  P.,  Lee,  S. H.,  Song,  J.,  Park,  S. K., and  Lim,  S. (2018) Decreased expression of serum- and glucocorticoid-inducible kinase 1 (SGK1) promotes alpha-synuclein increase related with down-regulation of dopaminergic cell in the Substantia Nigra of chronic MPTP-induced Parkinsonism mice and in SH-SY5Y cells. Gene 661, 189–195.
•  Yu,  G.,  Wang,  L.-G.,  Han,  Y., and  He,  Q.-Y. (2012) clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287.
•  Yuan,  N.,  Zhang,  G.,  Bie,  F.,  Ma,  M.,  Ma,  Y.,  Jiang,  X.,  Wang,  Y., and  Hao,  X. (2017) Integrative analysis of lncRNAs and miRNAs with coding RNAs associated with ceRNA crosstalk network in triple negative breast cancer. OncoTargets Ther. 10, 5883–5897.
•  Zhao,  C.,  Denison,  C.,  Huibregtse,  J. M.,  Gygi,  S., and  Krug,  R. M. (2005) Human ISG15 conjugation targets both IFN-induced and constitutively expressed proteins functioning in diverse cellular pathways. Proc. Natl. Acad. Sci. USA 102, 10200–10205.