Antidyskinetic Effects of MEK Inhibitor Are Associated with Multiple Neurochemical Alterations in the Striatum of Hemiparkinsonian Rats

L-DOPA-induced dyskinesia (LID) represents one of the major problems of the long-term therapy of patients with Parkinson's disease (PD). Although, the pathophysiologic mechanisms underlying LID are not completely understood, activation of the extracellular signal regulated kinase (ERK) is recognized to play a key role. ERK is phosphorylated by mitogen-activated protein kinase kinase (MEK), and thus MEK inhibitor can prevent ERK activation. Here the effect of the MEK inhibitor PD98059 on LID and the associated molecular changes were examined. Rats with unilateral 6-OHDA lesions of the nigrostriatal pathway received daily L-DOPA treatment for 3 weeks, and abnormal involuntary movements (AIMs) were assessed every other day. PD98059 was injected in the lateral ventricle daily for 12 days starting from day 10 of L-DOPA treatment. Striatal molecular markers of LID were analyzed together with gene regulation using microarray. The administration of PD98059 significantly reduced AIMs. In addition, ERK activation and other associated molecular changes including ΔFosB were reversed in rats treated with the MEK inhibitor. PD98059 induced significant up-regulation of 418 transcripts and down-regulation of 378 transcripts in the striatum. Tyrosine hydroxylase (Th) and aryl hydrocarbon receptor nuclear translocator (Arnt) genes were down-regulated in lesioned animals and up-regulated in L-DOPA-treated animals. Analysis of protein levels showed that PD98059 reduced the striatal TH. These results support the association of p-ERK1/2, ΔFosB, p-H3 to the regulation of TH and ARNT in the mechanisms of LID, and pinpoint other gene regulatory changes, thus providing clues for identifying new targets for LID therapy.

Other established molecular hallmark of LID is FosB, a truncated isoform of FosB, that is a chronic transcription factor.
We aimed at reducing the phosphorylation of ERK1/2 to examine its role in the regulation of other biomarkers of LID and associated transcriptional changes. We used the MEK inhibitor PD98059 and whole transcriptome analysis using microarray analysis in rats with unilateral 6-OHDA lesions and chronic exposure to L-DOPA treatment. In line with previous work (Santini et al., 2007), PD98059 prevented the phosphorylation of ERK1/2 and histone H3, and reduced the abnormal accumulation of FosB. Also, the inhibitor showed a powerful effect of counteracting AIMs. Importantly, numerous genes were regulated in association with AIMs. We focused on the function of the following genes, tyrosine hydroxylase (Th) and aryl hydrocarbon receptor nuclear translocator (Arnt), which regulates Th (Teh et al., 2007). The administration of PD98059 neutralized the L-DOPA-induced changes in Arnt and Th levels. These results provide evidence for a significant interplay between LID biomarkers (p-ERK1/2, FosB, p-H3) and the regulation of striatal Th and Arnt genes in the pathophysiology of LID development. The analysis of the transcriptome may also provide clues for identifying further genes that are involved in LID mechanisms.

Animals
Adult, male Sprague-Dawley rats (Beijing HFK Bioscience Co., Ltd., China) weighing 230-250 g (7 weeks old) were housed with free access to food and water, 12 h light/dark cycle, constant temperature and humidity. Animal use and care were conformed to the Guidelines of Laboratory Animals Ethics of Tongji Medical College, Huazhong University of Science and Technology. The protocol was approved by the Ethics Committee of Huazhong University of Science and Technology.

6-OHDA Lesion and Cannula Implantation
Rats were deeply anesthetized with isoflurane (induction 3%, maintenance 1.5%) in oxygen and mounted on the stereotaxic frame. A temperature controller system was used to maintain body temperature at 37 • C. 6-OHDA (2 µg/µl, 4 µl) was injected into the right medial forebrain bundle through a 10-µl microsyringe at a rate of 0.5 µl/min at the following stereotactic coordinates: relative to bregma, −4.4 mm anterior (A), +1.5 mm lateral (R), and 7.8 mm deepness (D) (Lindgren et al., 2009). In 33 rats, a guide cannula connected to a microinjection system (RWD Life Science Co., Ltd) was implanted to the right lateral ventricle (A −0.7 mm, R +1.5 mm, D 4.0 mm) according to atlas of Paxinos and Watson (2005) and secured to the skull with dental cement. The efficacy of the dopaminergic lesion was tested by measuring contralateral turning behavior with an acute subthreshold dose of apomorphine (0.05 mg/kg s.c.) 2 weeks post-surgery. Only rats exhibiting more than 200 turns contralateral to the lesion side in 30 min were considered to be compatible with the model of full lesion, and were chosen for further study (Boldry et al., 1995). TH staining revealed also a complete lesion of substantia nigra (SN) ipsilateral to the lesion (see Figure 1).

Drugs Treatment and Behavioral Assessment
Three days after apomorphine screening (Figure 2A), 45 successful hemiparkinsonian rats received intraperitoneal injections of L-DOPA plus benserazide (12/6 mg/kg) once daily. L-DOPA-induced AIMs were recorded every other day. Rats were observed for 1 min every 35 min intervals for a total of 140 min following L-DOPA treatment. AIMs were evaluated using the validated AIMs scale. Orofacial, limb, and axial dyskinesia were graded from score 0 to 4: 0 = absent; 1 = occasional, present during less than half min; 2 = frequent, present during more than half min; 3 = continuous but interrupted by strong sensory distraction; 4 = continuous, not interrupted by strong sensory distraction (Winkler et al., 2002;Lundblad et al., 2004). The total of axial, limb and orofacial dyskinesia was also named ALO dyskinesia and the maximum ALO dyskinesia score in each session was 48. Rotation (contralateral turns) test was performed as before (Breger et al., 2013) at day 1, 4, 12, 16, and 20, and only turns of completed 360 • were counted.
From day 10, 42 rats with complete dyskinesia (ALO AIMs score = 12 at 70 min after L-DOPA) were divided into 3 groups: LID group, LID + PD98059 group, and LID + vehicle group. There was cannula implantation for each rat in the latter two groups. A volume of 3 µl PD98059 solution (0.4 µg/µl) or vehicle was infused into the right lateral ventricle over 6 min through the cannula connected to the microinjection system half hour before L-DOPA administration from day 10 to 21 (Figure 2A). The infusion needle was left in place for 2 min at the end of infusion for drug diffusion. Normal rats with sham operation and hemiparkinsonian rats were chosen as normal group and PD group, respectively. There were 14 rats in each group.
The stepping test was carried out 15 min before (pre-) and after (post-) L-DOPA treatment as previously described (Pinna et al., 2007(Pinna et al., , 2010 on days 2, 10, and 18. Rats were moved on the surface of the table for 0.7 m in 4 s by the experimenter. The number of adjusting steps of left and right forelimbs in the forward directions was counted.
All behavioral tests were performed by a blinded examiner.

Tissue Preparation
At day 21, 4 h after the last L-DOPA treatment, 4 rats of each group (normal, PD, LID, LID + Vehicle, and LID + PD98059) were anesthetized with an overdose of pentobarbital and then perfused transcardially with saline followed by 4% icecold paraformaldehyde (PFA) in phosphate buffer (pH 7.4). Brain tissues were removed, post-fixed with 4% PFA overnight and 25% sucrose for 24 h. Then the fixed samples were embedded in paraffin for immunohistochemistry and immunofluorescence. Another 10 rats of each group were sacrificed by decapitation. Brain tissues were rapidly removed, dissected for striatum on powdered dry ice. Striatum from 4 of the 10 rats were immediately stored at −80 • C until protein extraction for Western blotting analysis, and 6 of the 10 were incubated in RNAlater (Qiagen) at 4 • C overnight and then stored at −80 • C for further Microarray analysis (3 rats) and real-time quantitative PCR (3 rats).

Immunohistochemistry and Immunofluorescence
The paraffin-embedded brain tissues were sectioned at thickness of 4 µm. SN and striatum sections were mounted on glass slides, deparaffinized by xylene for 15 min (2 times), dehydrated in graded ethanol solutions, baked in the basic antigen retrieval buffer (pH = 6.0), and washed with phosphate buffer (pH 7.4) for 5 min (3 times). After washing, sections were blocked with 3% Bovine Serum Albumin for 30 min at room temperature (RT), then incubated with diluted primary antibody in a humidified chamber overnight at 4 • C overnight.
The following primary antibodies were used: rabbit polyclonal antibody against FosB (1:100, sc-48, Santa Cruz, detecting FosB-FosB) for the immunohistochemistry staining in striatum, rabbit polyclonal antibody against TH (1:750, ab112, Abcam) for the immunohistochemistry staining in SN. Then all sections were washed with phosphate buffer (pH 7.4) for 5 min (3 times), subsequently incubated with biotinylated goat anti-rabbit IgG at 37 • C for 50 min, washed again as above, incubated with Horseradish peroxidase labeled streptavidin fluid at 37 • C for 30 min, washed, followed by DAB solutions for 5 min, washed, counterstained with Harris hematoxylin for 3 min, dehydrated in graded ethanol solutions, and eventually cover slipped. Images were collected through an Olympus camera connected to the microscope at the same light intensity, and analyzed using Image-Pro Plus software by an independent experimenter blinded to the sections by counting the number of positive cells in the SN or dorsal striatum of the lesioned side (relative to bregma, A −4.5 mm for SN and A 0.6 mm for striatum sections, N = 4 sections for counting). TH staining in SN was used to estimate the extent of the dopaminergic lesion (shown in Figure 1).
Immunofluorescence staining shared a same procedure with immunohistochemistry staining before secondary antibodies incubation. Primary antibodies used in the double staining contained rabbit polyclonal antibody against FosB-FosB (1:100, sc-48, Santa Cruz, detecting FosB-FosB) and mouse monoclonal against p-H3 (1:200, Ab14955, Abcam). Primary antibodies used in the single staining contained mouse monoclonal antibody against HIF-1β (1:200, ab2771, Abcam, detecting ARNT) and rabbit polyclonal antibody against TH (1:750, ab112, Abcam). After being washed, sections were incubated in dark with an appropriately diluted Alexa 488-or Cyanine 3-coupled secondary antibodies for 50 min followed by DAPI (4 ′ , 6-diamidino-2-phenylindole) dyeing nucleus for 10 min. Images were collected using laser confocal microscopy marked with image manipulation software, and analyzed using Image-Pro Plus software by an independent experimenter blinded to the sections by counting the number of positive cells in the dorsal striatum of the lesioned side (relative to bregma, A 0.6 mm for striatum sections, N = 4 sections for counting).

Microarray
Samples (striatum of lesioned side) were sent to Shanghai biotechnology Corporation for whole transcriptome analysis using microarray. Total mRNA was extracted using TRIZOL Reagent (Cat#15596-018, Life technologies, Carlsbad, CA, US) according to the manufacturer's instructions. RNA integrity was evaluated with the Agilent Bioanalyzer 2,100 and RNA 6,000 Nano/Pico Kit (Agilent Technologies, Santa Clara, CA, US). Qualified total RNA was further purified by RNeasy micro kit (Cat#74004, QIAGEN, GmBH, Germany) and RNase-Free DNase Set (Cat#79254, QIAGEN, GmBH, Germany).
Concentrations of extracted RNA were assessed with the Nanodrop spectrophotometer (Nanodrop Technologies). A total of 100-150 ng RNA per sample was reverse transcribed to double stranded cDNA and then transcribed into cRNA using the Genechip WT Expression Kit (Affymetrix). Second cycle was carried out following generation of cRNA in order to transform the cRNA into single-strand cDNA. The cDNA was fragmented and the Genechip WT Terminal Labeling Kit (Affymetrix) was used to label the single-stranded DNA with biotin. Samples were hybridized to an Affymetrix Genechip Rat Gene 2.0 ST Array Platform. Array hybridization and wash was performed using GeneChip R Hybridization, Wash and Stain Kit (Cat#900720, Affymetrix, Santa Clara, CA, US) in Hybridization Oven 645 (Cat#00-0331-220V, Affymetrix, Santa Clara, CA, US), and Fluidics Station 450 (Cat#00-0079, Affymetrix, Santa Clara, CA, US) according to the manufacturer's instructions. Slides were scanned by GeneChip R Scanner 3000 (Cat#00-00212, Affymetrix, Santa Clara, CA, US).

Data Analysis
The Expression Console software (Affymetrix, Santa Clara, CA, US) was used to format the raw microarray data. Data pre-processing, including background adjustment, log fold transformation and normalization was completed using the "exon level" option in the software. Moreover, the dataset was normalized using the Robust Multi Array Average (RMA) method to control the inter-array variability. Normalized signal intensities of probes which belong to one transcript of each sample were processed by median for further data analysis. Gene chip and RNA quality were assessed by examining total mRNA expression for each striatum. We have submitted the microarray data to the GEO repository. The GEO accession numbers is GSE93695.
R software was used to screen the differential gene expression transcripts (DETs) among these samples. The adjusted P < 0.05 and fold change (FC) ≥ 1.2 or ≤ 0.8333 were used as the cut-off criteria. After getting the DETs, we proceeded with the functional analysis, searching the function and signaling pathway of the genes. In this step, we used the known databases mainly, including DAVID Gene Ontology (http://david.abcc. ncifcrf.gov/). Considering the large amount and complex branch structure of Gene ontology (GO) biological processes, we used a significance threshold P < 0.05 for biological process terms.
The co-expression analysis starts by constructing a matrix of pairwise correlations between all pairs of transcripts across samples of LID + PD98059, LID + Vehicle and LID group. We built an unsigned co-expression network with all the 37,177 transcripts in microarray using the Weighted Gene Correlational Network Analysis (WGCNA) package [PMID: 19114008]. GO term enrichment tests were performed for individual gene coexpression modules compared to a background set of all genes expressed in these brain samples using the R packages GOstats (version 2.26.0), biomaRt version (2.14.0), AnnotationDbi (version 1.20.7), and org.Hs.eg.db (version 2.8.0).
To specialize LID + PD98059 from other treatment, unsupervised hierarchical clustering of the candidates in every trait module was performed by bootstrapping analysis using MeV software (http://www.tm4.org/). Bootstrapping analysis provides confidence values for the stability of each cluster derived by hierarchical clustering.
RT-qPCR data were expressed as fold changes in relative gene expression compared with the Normal group using beta-actin levels as an endogenous control. The significance level was set at p < 0.05. Data are presented as mean ± SEM.

Statistics
Data were analyzed using one-way analysis of variance (ANOVA) followed by Tukey HSD or LSD post-hoc tests for multiple comparisons between groups, and Student's t-test for comparing TH positive neurons in the intact and lesioned side of SN. All statistical analyses were performed in SPSS 21.0 software. The significance level was set at p < 0.05. Data are presented as mean ± SEM for behavioral assessments, immunohistochemical count, and blot quantifications.
The Stepping test showed that L-DOPA could significantly improve motor function of the contralateral forelimb in rats with unilateral nigrostriatal 6-OHDA-lesion (Figures 2F,G). PD98059 did not affect the stepping test scores (Figures 2F,G).

Genes Associated with LID and Regulated by MEK Inhibitor
Gene expression was compared between the lesioned striatum of rats without L-DOPA treatment (PD group) and that of rats with AIMs after repeated administrations of L-DOPA for 21 days (LID group; ALO AIMs score = 12 at 70 min after L-DOPA). Compared with PD group, LID group exhibits 606 upregulated transcripts and 932 down-regulated transcripts (p < 0.05 and fold change ≥ 1.2 or ≤ 0.8333; Figure 5A). These  changes in gene expression are caused by L-DOPA treatment and are mainly related to the following functions, response to endogenous stimulus, enzyme linked receptor protein signaling pathway and several others described in Figure 5B (see also Supplementary Table 2 for details).
In comparison with LID + Vehicle group, LID+PD98059 group exhibits 418 up-regulated transcripts and 378 downregulated transcripts (P < 0.05 and fold change ≥ 1.2 or ≤ 0.8333; Figure 5A). Thus, PD98059 treatment caused large changes in gene expression mainly related to the following functions, negative regulation of multicellular organismal process and several others described in Figure 5B (see also Supplementary Table 4 for details).
Changes in gene regulation between PD and Normal groups and between LID + Vehicle and LID groups were also detected (see Figure 5 and details in Supplementary Tables 1, 3).
Changes in Th, Arnt, and FosB gene expression detected with microarray analysis were verified with quantitative RT-PCR (Figures 6C-E). Similar to microarray results, Th and Arnt were up-regulated by L-DOPA and normalized by PD98059 (P < 0.05). Only one of the four probes detecting FosB gene (Probe Set ID: 17630237) presented a trend in line with expectations shown in the quantitative RT-PCR result (Figure 6E), namely up-regulated by L-DOPA and down-regulated by PD98059. The gene symbol is set as LOC100360880 in the data, and probe details are shown in the Supplementary Table 5.

Identification of PD98059-Specific Coexpression Modules
To further identify a hierarchical network view of co-expressed genes across LID+PD98059 group, LID + Vehicle group and LID group subtypes, we applied WGCNA to a dataset containing 3 specimens in each group. All candidates (including protein coding genes and lncRNA) in this dataset were hierarchically clustered under unsupervised average linkage and classified into 45 modules (Figures 7A,B) labeled by color. Each module was comprised of mutually exclusive co-expressed candidates. Candidates with no distinct module assignment were grouped in a gray module by WGCNA. Two of these modules, containing steelblue and red, were identified using any pre-assigned phenotype (GS P < 0.05, Figure 7C). The module membership (MM) vs. gene significance (GS) plots for these modules ( Figure 7C) showed that MM and GS are highly correlated, indicating that the candidates most significantly associated with the trait are often also the most important (central) elements of the respective module. Following the unsupervised module generation, individual candidate correlations to a specific treatment were quantified by GS. The average GS of all candidates within each module is summarized in Figure 7D. This analysis unveiled positive or negative correlation of certain modules with PD98059 treatment. The steelblue module contained candidates negatively correlated to the PD98059 treatment ( Figure 7D). To validate the robustness of the co-expression network as a specific classifier, it was first applied by unsupervised hierarchical clustering bootstrap analysis to the expression value of each sample in a test dataset from which steelblue module were derived, making it a "PD98059-treatment module, " including protein metabolic process and RNA modification terms with GO enrichment analysis (Figure 7E). Similarly, the red module contained candidates positively correlated to the PD98059 treatment ( Figure 7D). And the genes in red module were involved in protein deacetylation and ubiquitin-dependent SMAD protein catabolic process (Figure 7E).

Identification and Validation of Hub Networks of Coexpressed Genes Closely Related to PD98059 and LID
Four genes in the steelblue module, namely LOC365559, Runx1, Far1, and Hspa4, and 9 genes in the red module, namely Cadm2, Kif6, Ephb2, Map2k6, Zfp316, Hook2, Tshz3, RGD1564887, and Th, were also common specific DETs between LID + PD98059 vs. LID + Vehicle and LID vs. PD (red and blue font, respectively; Figure 6B). This overlap suggests that the genes are crucially regulated by PD98059 treatment and LID. We further analyzed these genes in hub networks of co-expressed genes across transcriptome platforms. There were more than 60 DETs in the steelblue module closely related with Hspa4 ( Figure 8A). There were also many DETs in the hub network of the red module related with Map2k6, Hook2, Zfp316, Kif6, and Th ( Figure 8B). The regulation of these genes may thus be related to LID development.

Effects of MEK Inhibitor on the Striatal Expression of TH and ARNT
The expression of TH and ARNT was also analyzed using Immunofluorescence (Figure 8). Chronic L-DOPA administration increased the number of TH and ARNT positive neurons in the striatum of 6-OHDA-lesioned rats (Figure 9). This effect was reversed by the inhibitor of MEK PD98059 [ Figure 9; TH, F (3, 12) = 34.268; HIF-1 β (ARNT), F (3, 12) = 28.696; one-way ANOVA followed by Tukey HSD and LSD post-hoc test, N = 4/group]. No difference was found between LID + Vehicle and LID groups.

DISCUSSION
In this study, we showed that lateral ventricle administration of PD98059, a selective MEK inhibitor, could reduce LID in 6-OHDA-lesioned rats. Axial, limb, and orofacial AIMs scores   were significantly decreased since the fourth day of PD98059 administration. In addition, rotational behavior had little changes with the addition of the inhibitor indicating that this agent did not neutralize the antiparkinsonian effect of L-DOPA. These findings are in line with previous results showing that SL327, other MEK inhibitor, could significantly reduce LID in 6-OHDAlesioned mice (Santini et al., 2007) and definitely confirm the key role of ERK1/2 in the development of rodent LID. As classical inhibitors of MEK, PD98059 was widely used in a variety of experiments. It was proved to be highly effective in preventing the phosphorylation of ERK1/2 that could disrupt drug-paired contextual cue memories (Miller and Marshall, 2005). In this study, the increased striatal phosphorylation of ERK1/2 in hemiparkinsonian rats chronically exposed to L-DOPA could also be largely reversed by PD98059.
One of the most significant results of this study is that the antagonistic effect of PD98059 on LID and p-ERK1/2 was accompanied by other specific molecular changes, particularly the reduction of FosB expression. We have confirmed for the first time that inhibiting the over-activation of ERK1/2 could decrease the accumulation of FosB protein in the striatum. The FosB gene is ranked in the first place among 28 genes blocked by SL327, other MEK inhibitor with antidyskinetic effects in mice, which is consistent with the present results showing that MEK inhibitor could reduce the transcription of FosB (Charbonnier-Beaupel et al., 2015). However, we found in the microarray results that probes for detecting the FosB gene (mRNA Accession: NM_001256509) were probably not perfectly designed, although Affymetrix Genechip Rat Gene 2.0 ST was one of the most advanced whole transcriptome gene chip when we started the study. Only one of the four probes (Probe Set ID: 17630237) presented a trend in line with expectations, namely up-regulated by L-DOPA and down-regulated by PD98059. The gene symbol is set as LOC100360880 in the data and the probe details are shown in the Supplementary Table 5. Furthermore, we carried out realtime quantitative PCR to detect the FosB levels in samples. The change pattern among groups confirmed the result of Probe 17630237.
The transcription factor FosB is a recognized hallmark of LID (Andersson et al., 1999;Tekumalla et al., 2001;Cenci, 2002;Cao et al., 2010). Striatal overexpression of FosB could result in abnormal neuronal electrical properties leading not only to dyskinesia but also to an insensitive response to L-DOPA that can be responsive to selective blockade of the expression of this transcription factor (Engeln et al., 2016). The kinetic profile of FosB likely explains the delayed antidyskinetic effect of PD98059 to the fourth day of treatment. This transcription factor possesses unique stability properties compared with all other Fos family members (McClung et al., 2004;Alibhai et al., 2007;Nestler, 2008Nestler, , 2015, and this is owned to two contributing mechanisms. First, FosB lacks two critical degron domains, which are associated with ubiquitination and degradation in the C-terminus of full-length FosB and all other Fos family proteins. Second, FosB is modified by some protein kinases at the Nterminus, adding further stability to the protein (Ulery et al., 2006;Gajewski et al., 2009;Ulery-Reynolds et al., 2009;Cates et al., 2014). Therefore, FosB proteins tend to accumulate in neurons, and could persist for several weeks after withdrawal of relevant drug exposure. In this study, AIMs were significantly reduced starting several days after initiating MEK inhibitor treatment, and the FosB level was significantly decreased at the end of the experiment. Further studies are necessary to explore the role of certain kinases that mediate FosB catabolism and the contribution of ERK1/2 in the process.
Our microarray data showed that changes in the expression of a considerable number of genes occur with LID development and can be changed by MEK inhibition. Compared with the normal state, there were a total of 961 transcripts changed as a result of dopaminergic lesion. The regulated genes are related directly to synaptic transmission, nerve impulse conduction, behavior and other functions implying that dopamine depletion causes extensive transcriptional alterations. Furthermore, changes in 1,538 transcripts were involved in the development of LID after chronic L-DOPA treatment, and changes in 796 transcripts were associated with PD98059 treatment. The genes regulated with LID and PD98059 widely participate in transmembrane receptor protein tyrosine kinase signaling pathway, synaptic transmission, synaptic plasticity regulation, Ras protein signal transduction regulation, and other relevant functions. Notably, L-DOPA and PD98059 could affect some of the same GO biological processes. The preliminary analysis of GO biological processes, including enzyme linked receptor protein signaling pathway, synaptic transmission, cell-cell signaling, and learning and memory support that these genes, especially Hspa4, Map2k6, Hook2, Zfp316, Th, and likely other genes in the hub networks play important roles in the pathophysiology of LID. We focused on the most significant changes in gene expression occurred with PD98059 treatment. According to the in-depth analysis of the data, we found 13 genes that are up-regulated by L-DOPA and could be neutralized by PD98059, and 17 genes that are down-regulated by L-DOPA and could be normalized by PD98059.
It is important to note that among the 13 genes up-regulated by L-DOPA and neutralized by PD98059, one was the abnormal Th expression. Several studies have previously demonstrated that striatal Th-positive neurons are correlated with severity of LID and FosB expression in hemiparkinsonian mice (Darmopil et al., 2008;Charbonnier-Beaupel et al., 2015;Keber et al., 2015). Studies have shown that striatal Th neurons may cooperate with serotonergic terminals synthesizing dopamine and producing supraphysiological synaptic DA concentrations, a mechanism thought to contribute to LID (Keber et al., 2015). However, it was also reported that the striatal Th-positive neurons co-expressed with dynorphin and enkephalin, suggesting that they are medium spiny neurons of the direct and indirect striatal output pathways (Darmopil et al., 2008). Clearly, Th colocalization with dynorphin is coherent with our data and the relation to LID because there is sufficient evidence in support of the direct pathway role in LID mechanisms. Instead, Th colocalization with enkephalin is at odds with the specific association of Th expression with LID mechanisms. It is important to consider that the function of newly developed Th expression in striatal neurons in models of PD is still under poorly understood, and it is possible that Th expression serves different functions including a compensatory mechanism to enhance dopamine actions in both striatal pathways. On the other hand, the use of D1 or D2 receptor knock-out mice demonstrated that D1R, but not the D2R is necessary for L-DOPA-induced expression of striatal Th-positive neurons (Espadas et al., 2012). Although, transgenic models of dopamine receptor KO may undergo different gene regulations with L-DOPA treatment, and here we used wild type animals, the present data also support the link of striatal Th expression to LID. In our study, PD98059 significantly neutralized the abnormal Th expression, which suggests a connection between Th and MEK, and possibly p-ERK1/2. Of interest, p-ERK1/2 could indirectly influence Th expression through regulation of the orphan nuclear receptor Nur-related factor 1(NURR1) because putative ERK1/2 phosphorylation sites were found proximal to the N-terminal AF-1 region of NURR1, and NURR1 directly induces transcription of Th gene in midbrain dopamine neurons of the substantia nigra (Jacobsen et al., 2008). However, we have not found a distinct change in the Nurr1 gene with microarray among the studied groups of rats. Therefore, further work is needed to determine whether the connected regulation of p-ERK1/2-NURR1-Th plays a mechanistic role in LID.
Another regulated gene in relation to LID and PD98059 was Arnt, whose official full name is aryl hydrocarbon receptor nuclear translocator. ARNT protein is required for activity of the Ah (dioxin) receptor and the ligand-binding subunit to translocate from the cytosol to the nucleus after ligand binding. Then, the complex initiates transcription of genes. Its heterodimer with HIF1A acts as a transcriptional regulator in response to hypoxia (Mannello et al., 2011;Dela Cruz et al., 2014). It was revealed that ARNT, in concert with neuronal PAS domain protein 1 (NPAS1), negatively modulates the expression of Th and that this regulation occurs with NPAS1 directly binding on the Th promoter (Teh et al., 2007). However, this is conflicting with our data showing that Arnt gene expression correlates with Th expression changes associated with L-DOPA and PD98059 treatment, upregulation and suppression, respectively. Alternatively, Th expression may be induced by other transcription factors, for instance NURR1, and Arnt regulation could result from adaptive changes to counteract high Th expression.
Ubash3b and Sik1 genes were regulated in parallel to Th and Arnt in this study. The official full name of Ubash3b is ubiquitin associated and SH3 domain containing B. This gene encodes a protein that promotes accumulation of activated target receptors on the cell surface, exhibits tyrosine phosphatase activity, and down-regulates proteins that are dually modified by both protein tyrosine phosphorylation and ubiquitination. One possibility would be that UBASH3B contribute to the accumulation of FosB due to UBASH3B mediated accumulation of activated target receptors. The official full name of Sik1 is salt-inducible kinase 1. Sik1 is one of the CREB-target genes. SIK1 induction is thought to act as a negative feedback signal preventing persistent CREB/TORC1-dependent transcription in situations of long-lasting neuronal activity (Hu et al., 2015). It has been shown that Sik1 plays a key role in cocaine addiction (Dietrich et al., 2012), which also shares a mechanism associated with FosB upregulation. Map2k6 and Ephb2 are among the 17 genes down-regulated by L-DOPA and up-regulated by PD98059. The official full name of Map2k6 is mitogenactivated protein kinase kinase 6. MAP2K6 is the activator of p38 MAPK and is involved in synaptic plasticity including cue-induced relapse to heroin seeking, learning and memory (Bolshakov et al., 2000;Fanous et al., 2013). The official full name of Ephb2 is Eph receptor B2. EphB2 is a member of the EphB family of receptor tyrosine kinases, and could increase synaptic NR1 and NR2B expression, prevent downregulation of dephosphorylated p38 MAPK and phosphorylated CREB in Aβ1-42 oligomer-treated neurons (Geng et al., 2013). EphB2 is considered a neuroprotective factor for hippocampal neurons with a potential therapeutic role in Alzheimer's disease (Miyamoto et al., 2016). LID is thought to share similar plasticity changes with drug addiction and memory/learning processes. The present study showed that the Ubash3b, Sik1, Map2k6, and Ephb2, all associated with mechanisms of synaptic plasticity, are regulated by L-DOPA as well as PD98059 treatment.
Hspa4, Hcrtr1, Fgfr2, and Cadm2 are also significantly regulated with L-DOPA and PD98059 treatment. The official full name of Hspa4 is heat shock protein family A member 4, a member of the heatshock protein (HSP) family 110 (Banduseela et al., 2013). Most HSPs act as chaperones and imperfections in their function can lead to an accumulation of misfolded proteins (Patterson, 2006). HSPA4 has interaction with PD-causing genes, including parkin, DJ-1 and PINK1, although the role of HSPA4 in PD progression and pathogenesis is unclear (van der Merwe et al., 2015). The official full name of Hcrtr1 is hypocretin receptor 1. HCRTR1 is one of G-protein coupled receptor, and takes part in positive regulation of ERK1/2 cascade and cytosolic calcium ion concentration. The official full name of Fgfr2 is fibroblast growth factor receptor 2. FGFR2 is an important neurotrophic factor, and may contribute to the development of mesencephalic dopaminergic neurons (Baron et al., 2012). The official full name of Cadm2 is cell adhesion molecule 2. CADM2 plays an important role for synapse organization, and provides regulated trans-synaptic adhesion (Tanabe et al., 2013). Clearly, a number of other genes could be listed, but these gene regulations need to be investigated in depth to establish their specific relation to behavioral changes and to gain insights into their functional roles.
In conclusion, our results reveal the interplay among p-ERK1/2, FosB, p-H3, TH, ARNT in the mechanisms of LID. Additionally, the series of genes identified in this study are opening new pathways for further insights into the molecular changes associated with LID. Most noticeably are the interrelated regulation of p-ERK1/2-NURR1-TH-ARNT genes and p-ERK1/2-FosB-UBASH3B genes. Further studies of these mechanisms may help identify targets for developing new LID treatments.

AUTHOR CONTRIBUTIONS
GC collected the data, clustered the literature data, and wrote the paper. SN collaborated in preparation of figures writing the paper. CH and KM collaborated in establishing the model. YX and XC collaborated in writing the paper, supervised the overall work, and revised the paper. ZZ collaborated in referencing the literature and revised the paper. SP advised for data presentation and revised the paper.

ACKNOWLEDGMENTS
This study was supported by the National Natural Science Foundation of China (NSFC Project No. 81171193, No.30700881, and No.81571249). SP was supported by NIH grants NS045962, NS073994, NCRR RR000165, and ORIP/OD OD011132. Dr. Hanqi Yin (Shanghai Biotechnology Corporation) provided substantial help in microarray data analysis.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins. 2017.00112/full#supplementary-material