Parkinson’s disease (PD) is a progressive neurological disorder that results primarily from the death of dopaminergic (DAergic) neurons in the substantia nigra. The factors that trigger cell death in PD are currently unknown although neuromelanin accumulation, mitochondrial dysfunction, oxidative stress, exposure to iron and other metals, mutations in the α-synuclein gene, trauma, and dysfunction of the ubiquitin-proteasome system have all been implicated in the pathogenesis of PD (G?lvez-Jim?nez, 2007). Since the discovery that people who are intoxicated with 1-methyl-4-phenyl-1, 2, 3, 6-tetrahydropyridine (MPTP) develop a syndrome nearly identical to PD (Langston et al., 1983), MPTP has been used to generate PD models in non-human primates and mice. In the brain, MPTP is oxidized to 1-methyl-4-phenyl-2,3-dihydropyridinium (MPDP⁺) by monoamine oxidase B (MAOB) in glia and serotonergic neurons, which is then converted to 1-methyl-4-phenyl-pyridium (MPP⁺), an active metabolite of MPTP. MPP⁺ is taken up by DAergic neurons via dopamine and noradrenaline transporters, which results in inhibition of complex I of the mitochondrial electron transport chain and the formation of reactive oxygen species (ROS) (Lotharius and O’Malley, 2000; Nakamura et al., 2000), which in turn leads to cellular dysfunction and cell death (Nicotra and Parvez, 2002).

In vitro models such as pheochromocytoma (PC12) cells and human neuroblastoma (SH-SY5Y or SH-EP) cells have been widely used in studies to understand the pathogenesis of PD. When these cells were treated with neurotoxins such as MPP⁺, 6-hydroxydopamine (6-OHDA), or rotenone, many aspects of the DAergic neuron death are similar to those observed in PD. Conn et al. (2003) explored the transcriptional response of SH-SY5Y cells to MPP⁺ using cDNA microarray and identified transcription factors and cell cycle proteins as major classes of differentially expressed proteins. DNA microarray analysis is capable of profiling the expression levels of a large number of genes simultaneously and is therefore a promising technology for the elucidation of gene interactions. However, in many cases DNA microarray experiments yield a list of differentially expressed genes (DEGs) and it is not easy to link these DEGs to phenotypes that are the results of thousands of complex interactions occurring in various metabolic and signaling pathways. This indicates that the complex interactions in living cells should be taken into account to fully understand DNA microarray data. Recently, many efforts have been made to identify which biological pathways are perturbed using gene expression profiles (Pepe and Grassi, 2014; Tarca et al., 2009). Information for each pathway can be obtained from pathway databases such as KEGG (Kanehisa et al., 2012) and Reactome (Croft et al., 2014), which describe metabolic pathways and gene signaling networks and offer the potential for more complex and useful analysis. Identifying which pathways are perturbed in PD is important for improving our understanding of the molecular mechanisms underlying PD. In this study, I identified pathways that were significantly perturbed after treatment with MPP⁺ using the list of DEGs and their expression values with the commercial software Pathway Guide 3.0 (http://www.advaitabio.com).

Perturbed pathways at 1.5, 3, 9, 12, and 24 h after MPP⁺ treatment were examined using DEGs identified from whole genome expression data measured at each time point and their expression values. Two KEGG pathways, the mitogen-activated protein kinase (MAPK) signaling pathway and endoplasmic reticulum (ER) protein processing pathway, showed significant perturbation persistently and both resulted in positive perturbation of DNA-Damage-Inducible Transcript 3 (DDIT3). This provides a plausible explanation for the strong upregulation of DDIT3 observed in both microarray and real-time quantitative PCR (qPCR) analyses. Information on pathways that are perturbed by MPP⁺ treatment might provide insight into MPP⁺-induced neuronal cell death and the pathogenesis of PD.

Cell culture and MPP⁺ treatment

Human neuroblastoma SH-EP cells were kindly provided by Dr. Talia Hahn at Kaplan Medical Center (Rehovot, Israel). For treatment with MPP⁺ iodide (Sigma-Aldrich, USA), 1 × 10⁶ cells were plated in 100 mm² dishes (Corning, USA) in 10 ml Dulbecco’s modified Eagle’s Medium (DMEM; Sigma-Aldrich) with 10% fetal bovine serum (FBS), 100 units/ml penicillin, and 100 mg/ml streptomycin at 37°C with 5% CO₂ and cultured for 3 days. Freshly prepared MPP⁺ toxin was added to the cultures to a concentration of 1.25 mM and incubation continued at 37°C for 0 (control), 1.5, 3, 9, 12, and 24 h. The experiment was performed two times.

Assessment of cell viability

Cell viability was measured using the quantitative colorimetric MTT assay, which reveals the mitochondrial activity of living cells, as described previously (Nanjo et al., 1996). Briefly, MTT dissolved in phosphate-buffered saline was added at the end of incubation to a final concentration of 0.5 mg/ml. After 4 h incubation at 37°C and 5% CO₂, the supernatants were removed, and the formazan crystals that formed in the viable cells were measured at 550 nm using a microplate reader (Molecular Devices, USA).

RNA isolation and quality control

At each time point after MPP⁺ treatment, total RNA was extracted using TRIzol? (Invitrogen Life Technologies, USA) and purified using RNeasy columns (Qiagen, USA) according to the manufacturers’ protocol. Purified RNA samples with DNase treatment were quantified. All aliquots were stored at ?80°C until use. For quality control, RNA purity and integrity were evaluated by denaturing gel electrophoresis, optical density comparison of the 260/280 ratio, and analysis using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA).

Microarray expression profiling

Samples were collected at 0 h (control, before MPP⁺ treatment) and at 1.5, 3, 9, 12, and 24 h after MPP⁺ treatment. Total RNA was isolated from each sample, amplified, and purified using the Ambion Illumina RNA amplification kit (Ambion, USA) to generate biotinylated cRNA. Briefly, total RNA was reverse-transcribed to single stranded cDNA using a T7 oligo(dT) primer, converted into double-stranded cDNA, and purified. An in vitro transcription reaction was then carried out in the presence of biotinylated UTP and CTP to produce biotin-labeled cRNA from double stranded cDNA. After purification, the cRNA was quantified using a ND-1000 Spectrophotometer (NanoDrop, USA) and 750 ng of labeled cRNA was hybridized to each human HT-12 expression v.4 bead array for 16?18 h at 58°C, according to the manufacturer’s instructions (Illumina). Detection of the array signal was performed using fluorolink streptavidin-Cy3 (GE Healthcare Bio-Sciences, UK) according to the bead array manual. Arrays were scanned with an Illumina bead array reader confocal scanner according to the manufacturer’s instructions.

Normalization of microarray data

Probe signals exported from GenomeStudio software v2011.1 (Illumina) were checked for detection against negative controls with a GenomeStudio internal algorithm and missing values were introduced to replace signals under the detection limit. After the probe signals were log₂ transformed, quantile normalization was performed with the function lumiN from the lumi bioconductor package (www.bioconductor.org). As these normalized values were presented in terms of probe IDs, among which multiple IDs might correspond to the same gene, the probe IDs targeting the same gene were averaged and conversion of the probe ID to the gene name was performed using the lumi bioconductor package. After removing the probe IDs with no gene name, 48,803 probe IDs were finally mapped to 32,421 genes with a unique gene name (i.e., gene symbol).

Selection of differentially expressed genes

Genes that were differentially expressed between MPP⁺ treated and untreated (control) samples were identified by significance analysis of microarrays (SAM) (Tusher et al., 2001), which assigns a scored_i to each gene i based on the change in gene expression relative to the standard deviation of repeated measurements. SAM performs permutations of the repeated measurements to estimate the percentage of genes identified by chance, termed the false discovery rate (FDR). As the number of MPP⁺-treated and untreated samples is two for each time point in this study, the total number of distinct permutations is 6 (C²₄ = 6) for each time point. TheΔ value, which represents the distance from the expected d_i computed from permutation, was used for threshold for significance. A SAM analysis was carried out using samr in the R package (www.r-project.org). Genes with > 1.5-fold change andΔ = 0.4 were considered significantly regulated.

Identification of pathways that are perturbed by MPP⁺ treatment

Significantly perturbed pathways in MPP⁺ treated human neuroblastoma SH-EP cells were identified using the commercial software Pathway Guide 3.0 (http://www.advaitabio.com/). This program explores significantly affected pathways based on two types of evidence: (i) over-representation of DEGs in a given pathway, and (ii) abnormal perturbation of that pathway, measured by propagating measured expression changes across the pathway topology. These two types of evidence are captured by two independent probability values, P_NDE and P_PERT. This method provide the possibility to account for system-level dependencies and interactions, as well as identify perturbations and modifications at the pathway by overcoming the limitation of the most currently available approaches such as over-representation analysis (Dr?ghici et al., 2003) and functional class scoring (Goeman et al., 2004), which consider only the set of genes on any given pathway and ignore their position in those pathways.

The first probability, P_NDE= P(X ≥ N_DEH₀), captures the significance of the given pathway P_i as provided by an over-representation analysis (ORA) (e.g., Doniger et al., 2003; Pandey et al., 2004) of the number of DEGs (N_DE) observed on the pathway. H₀ stands for the null hypothesis: that the genes that appear as DEG on a given pathway are completely random. The P_NDE value represents the probability of obtaining a number of DEGs on the given pathway at least as large as the observed one, N_DE. The second probability, P_PERT, is estimated based on the amount of perturbation measured in each pathway. To calculate P_PERT for a given pathway, a gene perturbation factor (PF) is defined as described by Tarca et al. (2009): PF(gi)=ΔE(gi)+∑j=1nβijPF(gi)Nds(gj)

where ΔE(g_i) is the signed normalized expression change in the gene g_i (log₂ fold-change in this study) and the second term is the sum of the perturbation factors of the geneg_jdirectly upstream of the target gene g_i, normalized by the number of downstream genes of each such gene N_ds(g_j). The absolute value of β_ij quantifies the strength of the interaction between genes g_j and g_i. In this study, the value of β was +1 for induction/activation, and ?1 for repression/inhibition. A non-zero value of β was only assigned for genes that directly interacted with the gene g_i according to the pathway description. Given an edge directed from gene A to gene B in the pathway, A is considered upstream of B or B is considered downstream of A. The net perturbation accumulation at the level of each gene, ACC_g, is calculated as the difference in the perturbation factor PF of a gene and its log₂ fold-change: Acc(gi)=PF(gi)?ΔE(gi)

The total net accumulated perturbation of a given pathway is calculated as the sum of all perturbation accumulations for all genes in the pathway: tA=∑iAcc(gi). The probability, P_PERT, is defined as the probability of observing a total accumulated perturbation of the pathway, T_A, more extreme than t_A just by chance: PPERT=P(TA≥tA|H0)

Where H₀ is the null hypothesis for the distribution of T_A values. The null distribution of T_A values is empirically determined using a bootstrap approach. In this study, the bootstrapping number of iterations is set to 10,000.

To identify significantly perturbed pathways in human neuroblastoma SH-EP cells after MPP⁺ treatment, the two types of evidence, P_NDE and P_PERT, are combined into one global probability, P_G. The global probability P_G is calculated as in Tarca et al. (2009): PG=Ci?Ci?ln(Ci)

Where c_i= P_NDE(i)·P_PERT(i). Because P_G is a combined probability value, it can be used not only to rank the pathways but also to choose a desired type I error level. Since more than a hundred pathways (e.g., 146 pathways for KEGG database) are tested simultaneously, small P_G values can occur by chance. Therefore, the significantly perturbed pathways were identified by controlling the FDR with the popular FDR algorithm (Benjamini and Yekutieli, 2001).

cDNA synthesis and real-time qPCR

cDNA was produced using the Superscript™II RT-PCR System (Invitrogen, Germany) according to the manufacturer’s recommendations for oligo(dT)₂₀ primed cDNA-synthesis. cDNA synthesis was performed on 500 ng of total RNA at 42°C. Using this first-strand cDNA, the gene expression levels of the 21 test genes and one control gene were quantified using TaqMan technology on a ABI PRISM 7900HT Sequence Detection System (Applied Biosystems, USA) in 384-well microtiter plates using a final volume of 10 μl Gene-specific primers and probes were available as TaqMan Gene Expression Assays (Applied Biosystems, Table 1). Optimum reaction conditions were obtained with 5 μl of Universal Master Mix (Applied Biosystems, USA) containing dNTPs with UTP, MgCl₂, reaction buffer and AmpliTaq Gold? DNA polymerase, 90 nM of primer(s) and 250 nM fluorescence-labeled TaqMan probe. Finally, 2 μl template cDNA was added to the reaction mixture. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal control. Amplifications were performed with a 10 min template denaturation step at 95°C, followed by 40 cycles at 95°C for 15 s and 60°C for 1 min. All samples were amplified in triplicate and data were analyzed with Sequence Detector software (Applied Biosystems). The relative mRNA levels were estimated by the 2^ΔΔCT method (Livak and Schmittgen, 2001).

Effect of MPP⁺ concentration on cell viability

To investigate the neurotoxicity of MPP⁺ to SH-EP cells, cell viability was measured 48 h after MPP⁺ treatment. Greater than 50% cell death was observed at MPP⁺ concentrations of 1.25 mM or higher (Fig. 1A). This indicates that 1.25 mM MPP⁺ is sufficient to induce apoptosis in SH-EP cells. To capture the early cellular events triggered by MPP⁺ treatment, it is important to control the exposure time to MPP⁺. The time-dependent cell viability for 1.25 mM MPP⁺ is shown in Fig. 1B. Cell viability greater than approx imately 78% could be achieved for up to 24 h of treatment. Therefore, in this study I examined gene expression and pathway perturbation for up to 24 h after MPP⁺ treatment.

Selection of DEGs from a time series microarray experiment

Gene expression for MPP⁺-treated SH-EP cells was measured at six time points between 0 and 24 h using the human HT-12 expression v.4 bead array that included 47,231 probes with well-established or provisional annotation. A total of 12 bead arrays were used for examination of gene expression at the six time points (0, 1.5, 3, 9, 12, and 24 h; two replicates for each time point). The bead-summary data were log₂ transformed and normalized using quantile normalization (see “Materials and Methods”). After averaging of multiple probe IDs targeting the same gene, 48,803 probe IDs were finally mapped to 32,421 genes with a unique gene name, i.e., gene symbol. DEGs between MPP⁺-treated and untreated SH-EP cells (control) were identified for each time point by SAM analysis with R package samr (see “Materials and Methods”).

The proportion of DEGs that were downregulated increased with time of exposure to MPP⁺ (0, 20, 55, 57 and 53% respectively for 1.5, 3, 9, 12 and 24 h) (Fig. 2). The fact that no DEGs with down-regulation were detected at 1.5 h suggests that the effect of MPP⁺ treatment might be initialized by activation of MPP⁺ responsive genes. The numbers of up-regulated DEGs were 844, 803, 567, 797, and 1,679 respectively for 1.5, 3, 9, 12, and 24 h (Fig. 2). Of these, 55 genes were commonly up-regulated at all time points, as shown in Table 2. As the expression of these genes was robustly upregulated during the time course, they might play a core role in neuronal cell death induced by MPP+. The continued upregulation of genes such as brain-derived neurotrophic factor (BNDF) and fibroblast growth factor 2 (FGF2) is particularly interesting. Almeida et al. (2005) reported that BNDF protects neurons via transient activation of the Ras/MAPK pathway and the PI3-K/Akt pathway. In addition, it was reported that intrastriatal infusions of FGF2 induced recovery of striatal DAergic fibers and DAergic content in the mouse MPTP model, and increased protection in the neurotoxicity induced lesion of nigrostriatal DA system (Date et al., 1993). Continued upregulation of the BNDF and FGF2 genes might be a cell survival strategy against MPP⁺ toxicity. Extraordinary high expression was observed at two genes including MIR1974 and CLDN1. MIR1974 is considered non-canocical miRNA because it maps to the mitochondrial tRNA and rRNA genes (Bandiera et al., 2011), and its upregulation was also detected in MPP⁺ treated SH-SY5Y cells (Kimet et al., 2013). CLDN1 plays a major role in tight junction-specific obliteration of the intercellular space and it has been reported that its reexpression induces apoptosis in breat tumor spheroids (Hoevel et al., 2004). The percentage of DEGs for each time point compared with all genes on the array (i.e., 32,241 genes) was 2.6, 3.1, 3.88, 5.69, and 10.9% respectively for 1.5, 3, 9, 12, and 24 h.

Analysis of perturbed pathways in MPP⁺-treated SH-EP cells by Pathway Guide 3.0

Identifying the pathways that are significantly affected by MPP⁺ exposure is a crucial step in understanding the underlying molecular mechanisms of MPP⁺-induced cell death. To examine pathways that are significantly perturbed by MPP⁺ treatment, two independent probability values, P_NDE and P_PERT, were calculated for each pathway with incorporating parameters such as the normalized fold change of the DEGs, statistical significance of the set of pathway genes, and the topology of the signaling pathway using the commercial software Pathway Guide 3.0 (http://www.advaitabio.com; see “Materials and Methods” for details). These two types of evidence were combined into a global probability value, P_G, which is used to rank the pathways and test the hypothesis that the pathway is significantly perturbed by MPP⁺ treatment. The perturbation of pathways at 1.5, 3, 9, 12, and 24 h after MPP⁺ treatment was examined for 146 well-characterized human gene signaling pathways in KEGG database. In Fig. 3, two types of p-value are illustrated in two-dimensional plots for each time point. In this figure, the horizontal axis is used to plot the ?log P_NDE (over-representation evidence) and the vertical axis is used to plot ?log P_PERT (perturbation evidence). Since several pathways are tested simultaneously, the significance should be considered after adjustment of the p value for multiple comparisons. Pathways above the oblique red line are significant at 5% after Bonferroni correction (Hochberg, 1988) whereas those above the blue line are significant at 5% after FDR correction (Benjamini and Yekutieli, 2001). The number of significantly perturbed pathways dramatically increased at 24 h. This suggests that it might be difficult to restore the normal state of cells after 24 h of MPP⁺ treatment.

The top 10 pathways ranked by P_G at each time point are shown in Table 3. When the threshold for significance after FDR correction was 5%, all pathways listed in Table 3 satisfied this criterion. In particular, the MAPK signaling pathway showed significant perturbation at all time points. Closer inspection of genes associated with this pathway at 9 h after MPP⁺ treatment showed that several genes encoding ligands, including BDNF, FGF2, and transforming growth factor beta 3 (TGFB3), and, voltage-dependent L-type calcium channel subunit beta-3 (CACNB3) were differentially expressed (Fig. 4A). As these genes are located at very upstream in the pathway and affect entry points controlling this pathway, their perturbations are widely propagated throughout the pathway. MAPK genes including MAP2K5, MAP3K6, MAP4K2, MAPK12, MAPK8IP1 and MAPK7 were downregulated whereas genes near the end of pathways such as activating transcription factor 4 (ATF4), avian myelocytomatosis viral oncogene homolog (MYC), and DDIT3 showed upregulation. Based on these observed expression values, the total perturbation of each gene and its propagation in the MAPK signaling pathway was estimated (see “Materials and Methods”) and is shown in Fig. 4B. The red arrows represent propagation of perturbation. Through propagation of the perturbation of upstream genes, downstream genes were either upregulated (MAPK8/9/10 and MAP3K4) or downregulated (MAP3K5). The MAPK signaling pathway was more severely perturbed at 24 h after MPP⁺ treatment (Fig. 5A). At this time point, many ligand-encoding genes such as nerve growth factor (NGF), BNDF, FGF2, interleukin 1, alpha (IL1A), and TGFB2/3 were perturbed, as well as genes encoding receptors including platelet-derived growth factor alpha (PDGFRA), Fas cell surface death receptor (FAS), and monocyte differentiation antigen CD14 (CD14), with upregulation of MAP2K4, MAP3K1, and MAP3K7. The propagation of IL1A perturbation via interleukin 1 receptor type I (IL1R1) induced upregulation of caspase 3 (CASP3), whereas perturbation of FAS and TGFB1 resulted in upregulation of death-domain associated protein (DAXX) (Fig. 5B). Perturbation of the MAPK signaling pathway for 9 and 24 h ultimately had the same result of upregulation of DDIT3, ATF4, and v-rel avian reticuloendothe-liosisviral oncogene homolog B (RELB). Upregulation of DDIT3 has been also reported in MPP⁺-treated human neuroblastoma SH-SY5Y cells (Conn et al., 2002).

The other pathway that was frequently perturbed was the ERprotein processing pathway (Table 3). Perturbation of genes involved in this pathway at 24 h is shown in Fig. 6A. Altered expression of ubiquitin ligase complex genes such as ring finger protein 5 (RNF5), VCP-interacting membrane protein (VIMP), DnaJ homolog, subfamily b, member 2 (DNAJB2), ubiquitin-conjugating enzyme E2G 2(UBE2G2), and F-box protein 6 (FBXO6) might be a response to the ER stress caused by misfolded proteins. In this pathway, the most strongly upregulated gene was DDIT3, with a log₂ fold change of 4.4. The B-cell lymphoma 2 (BCL2) also showed upregulation (log₂ fold change = 0.95), but is inhibited by DDIT3, therefore overall BCL2 is strongly negatively perturbed by the propagation of DDIT3 perturbation (Fig. 6B). Moreover, apoptosis induced through a cascade starting with Bcl-2 homologous antagonist/killer (BAK1) and ending with Caspase 12 (CASP12) might be inhibited because BCL2 inhibits BAK1. In addition, downregulation of ER-to-nucleus signaling 1 (ERN1) might inhibit apoptosis through a cascade starting from ERN1 and ending with MAPK8. DDIT3 was markedly positively perturbed by upstream genes including ATF4, nuclear factor erythroid 2-related factor 2 (NFE2L2), and X-box binding protein 1 (XBP1). Together, these data suggest that DDIT3 might play a key role in MPP⁺-induced neuronal cell death.

Among the pathways listed in Table 3, MAPK signaling pathway and ER protein processing pathway showed larger values oft_A than other pathways at 12 and 24 h after MPP⁺ treatment. In addition, these two pathways were significantly perturbed at all observed time points, implicating an important role in neuronal death induced by MPP⁺. The number of significantly perturbed pathways at a level of 5% after FDR correction was 26, 19, 13, 24, and 51 for 1.5, 3, 9, 12, and 24 h respectively. The dramatic increase in the number of perturbed pathways at 24 h indicates that the cytotoxic effect by MPP⁺ might be propagated broadly into several cellular regulation networks. Most of the significantly perturbed pathways showed a tendency to have perturbed genes at or near the entry point of the pathway, thus differential expression of the ligand or receptor might have a large impact on the downstream pathway.

Validation of microarray data by real-time qPCR

To increase the fidelity of the microarray data, all gene expression values were measured twice with two independent samples for each time point. However, the quality of gene expression data obtained from microarrays can vary greatly according to the platform and procedures used because microarray experiments measure the expression of tens of thousands of genes simultaneously. For validation of the microarray data, real-time qPCR was carried out for 21 genes randomly selected from significantly impacted pathways shown in Table 3 at two time points, 3 and 24 h after MPP⁺ treatment, with TaqMan Gene Expression Assays (Table 1). All real-time qPCR experiments were performed in triplicate for each sample. Although there were some differences in the absolute value of the fold change between microarray and real-time qPCR results, the pattern of regulation showed good agreement. That is, five genes that were upregulated genes in microarray data at 3 h-BNDF,SMAD4, CREB5, FGF2 and HIF1A-also showed upregulation in real-time qPCR data at 3 h with a fold-change threshold of 1.5 and 2, respectively, for microarray and real-time qPCR data (Table 4). In addition, genes that were downregulated in microarray data, such as DDIT4 and EGR1, also showed downregulation in real-time qPCR experiment. Similarly, the regulation pattern at 24 h showed good agreement between microarray and real-time qPCR data, but the fold change values for DDIT3 and IFNB1 in real-time qPCR were very high compared with those for microarray data. This might be due to the high signal detection sensitivity and accuracy of real-time qPCR. The dynamic range for most commercial microarray platforms is usually 3?4 orders of magnitude, whereas TaqMan-based real-time PCR can achieve a dynamic range of 6?8 orders of magnitude (Wang et al., 2006). Therefore, the fold-changes obtained from real-time qPCR were considered more reliable than those from microarray in terms of detection sensitivity and accuracy. Nonetheless, the reliability of gene expression data obtained from microarray is comparable to that of real-time qPCR data in terms of the pattern of regulation.

Pathways perturbed in MPP⁺-treated human neuroblastoma SH-EP cells were identified using genome-wide gene expression data at five time points (1.5, 3, 9, 12, and 24 h) after MPP⁺ treatment. Two types of evidence, P_NDE and P_PERT, were used to identify pathways that were significantly perturbed by MPP⁺ treatment. Two pathways, MAPK signaling pathway and ER protein processing pathways, were highly ranked at each time point (Table 3). Perturbation of each these two pathways caused the common effect of upregulation of DDIT3 (Figs. 4B and 5B), which was confirmed by real-time qPCR. Tabuchi et al. (2014) demonstrated that treatment of the cells with the siRNA for DDIT3 significantly elevated the cell viability in sodium floride-elicited cell death accompanying ER stress with the induction of DDIT. It was also reported that the knockdown of DDIT3 in neuroblastoma cells significanlty reduced glucose deprivation-triggered cell death (K?gel et al., 2006). These evidences support that MPP⁺ induces the apoptosis with ER stress as described in previous reports (Ghribi et al., 2003; Selvaraj et al., 2012). Upregulation of DDIT3 was previously reported in MPP⁺-treated SH-SY5Y cells, which were generated by sub-cloning of the neuroblastoma SK-N-H cell line (Conn et al., 2002). The SH-EP cells used in this study were also originated from the SK-N-H cell line, but are morphologically distinct from SH-SY5Y cells; SH-SY5Y cells are N-type cell line expressing neurotransmitters and various neuronal cell surface markers whereas SH-EP cells are an epithelial substrate-adherent Schwannian (S-type) cell line that expresses proteins characteristic of Schwann cells and lacks neuronal markers (Do et al., 2011). The upregulation of DDIT3 by MPP⁺ treatment in both SH-SY5Y and SH-EP cells suggests that DDIT3 might play a key role in MPP⁺-induced neuronal cell death. Conn et al. (2002) suggested that upregulation of DDIT3 in MPP⁺-treated SH-SY5Y cells might be a cellular mechanism that is distinct from mitochondrial impairment or oxidative stress. Our results suggest that DDIT3 could be strongly upregulated by perturbation of two pathways-MAPK signaling pathway and ER protein processing pathway.

Perturbation of the MAPK signaling pathway was mainly driven by differential expression of ligand-encoding genes such as NFG, BDNF, IL1A, and TGFB2, and receptor genes such as FAS and CD14. These genes were upregulated at 24 h after MPP⁺ treatment and their perturbations were propagated to downstream genes in the MAPK signaling pathway (Fig. 5B). This suggests that perturbation of the pathway entry point might be more influential than perturbation of downstream genes in the pathway. Perturbation of the MAPK signaling pathway ultimately resulted in upregulation of DDIT3 and ATF4. Perturbation of the ER protein processing pathway involved ubiquitin ligase complex genes as well as ER associated degradation (ERAD)-related genes (Fig. 6A). Absence of the ERAD proteins that remove misfolded proteins from the ER might result in the accumulation of misfolded proteins and activation of the unfolded protein response (UPR) by ER stress. Indeed, UPR pathway genes such as eukaryotic translation initiation factor 2 subunit 1 (EIF2S1), NFE2L2, and ATF4 were upregulated after MPP⁺ treatment although there was no change in the expression of eukaryotic translation initiation factor 2-alpha kinase 1 (EIF2AK1) (Fig. 6A). DDIT3, which plays a role in apoptosis during ER stress (Marciniak et al., 2004; McCullough et al., 2001; Oyadomari and Mori, 2004), might be strongly upregulated as a result of perturbation of upstream genes including ATF4, NFE2L2, and XBP1.

In summary, the very high expression of DDIT3 in MPP⁺-treated human neuroblastoma SH-EP cells might result from the perturbation of the MAPK signaling pathway and ER protein processing pathway. Overexpression of DDIT3 results in a decrease in BCL2 protein and overexpression of BCL2 blocks DDIT3-induced apoptosis (Gotoh et al., 2004; Matsumoto et al., 1996). In this study, BCL2 was upregulated (fold-change = 1.9) at 24 h after MPP⁺ treatment, but this expression value was 11 times lower than that of DDIT3. Thus, the blocking effect of BCL2 on DDIT3-induced apoptosis might be weak. DDIT3 might also lead to glutathione depletion by overexpression of cation transport regulator homolog 1 (CHAC1), which is under the regulation of DDIT3 during the UPR of ER stress (Mungrue et al., 2009). The CHAC1 family of proteins function as γ-glutamylcyclotransferases acting specifically to degrade glutathione (Kumar et al., 2012), and overexpression of CHAC1 (fold change = 4.4) were observed at 24 h after MPP⁺ treatment. As overexpression of DDIT3 leads to translocation of BAX protein from the cytosol to the mitochondria (Gotoh et al., 2004), the DDIT-mediated death signal is finally transmitted to the mitochondria (Oyadomari and Mori, 2004), which function as integrators and amplifiers of the death pathway. This suggests that the toxicity signal of MPP⁺ is initiated by mitochondrial dysfunction caused by inhibition of complex I of the electron transport chain (ETC) and might feed back to the mitochondria via ER stress. This positive feedback mechanism could contribute to amplification of the death signal.