Comparison of Perturbed Pathways in Two Different Cell Models for Parkinson's Disease with Structural Equation Model

Abstract

Increasing evidence indicates that different morphological types of cell death coexist in the brain of patients with Parkinson's disease (PD), but the molecular explanation for this is still under investigation. In this study, we identified perturbed pathways in two different cell models for PD through the following procedures: (1) enrichment pathway analysis with differentially expressed genes and the Reactome pathway database, and (2) construction of the shortest path model for the enriched pathway and detection of significant shortest path model with fitting time-course microarray data of each PD cell model to structural equation model. Two PD cell models constructed by the same neurotoxin showed different perturbed pathways. That is, one showed perturbation of three Reactome pathways, including cellular senescence, chromatin modifying enzymes, and chromatin organization, while six modules within metabolism pathway represented perturbation in the other. This suggests that the activation of common upstream cell death pathways in PD may result in various down-stream processes, which might be associated with different morphological types of cell death. In addition, our results might provide molecular clues for coexistence of different morphological types of cell death in PD patients.

1. Introduction

Parkinson's disease (PD) is a complex multifactorial disorder marked by extensive system-wide pathology, including a substantial loss of dopaminergic (DAergic) neurons in substantia nigra pars compacta (SNpc). This loss results in a concomitant depletion of dopamine (DA) in the striatum, impairing the nigrostriatal system such that the execution of coordinated movements is prevented (Ammal Kaidery et al., 2013). Although the pathogenesis of PD has not been elucidated, growing body of evidence suggests that oxidative stress and apoptosis may play a critical role in the accumulation of modified proteins and degeneration of DAergic neurons (Fernandez-Espejo, 2004). Increasing evidence indicates that different morphological types of cell death coexist in the brain of PD patients, all of which may result from the activation of common upstream programmed cell death pathways, which are mitochondrion dependent (Perier et al., 2012). This indicates that neuronal cell death in PD patients may not follow a unique mechanism. That is, various down-stream mechanisms may be responsible for neuronal cell death in PD. Therefore, an understanding of downstream molecular events in PD is essential for development of efficient PD drug development.

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 widely used as an animal model to replicate human PD symptoms in nonhuman primates and rodents. MPTP can rapidly cross the blood–brain barrier, and MPTP inside the brain is converted by MAO-B (monoamine oxidase) into a neurotoxin precursor, 1-methyl-4-phenylpyridinium (MPP⁺), which causes selective destruction of dopaminergic neurons in the SNpc (substantia nigra pars compacta) (Sharma and Taliyan, 2015). MPP⁺ induces neurotoxicity primarily by inhibiting complex I of the mitochondrial electron transport chain, resulting in ATP depletion and increased oxidative stress (Castro et al., 2012; Sharma and Deshmukh, 2014). Several cell culture studies of PD have reported other pathological effects of MPP⁺, which include synaptic dysfunction, α-synuclein aggregation, impairment of retrograde/anterograde transport, mitochondrial swelling and loss of cristae, upregulation of autophagy, lysosomal expansion, lysosomal dysfunction, and capase 3-independent cell death (Dagda et al., 2013). Thus, MPP⁺ has been commonly used in experimental PD both in vitro and in vivo models.

Increasing evidence indicates that different morphological types of cell death coexist in the brain of PD patients, all of which may result from the activation of common upstream programmed cell death pathways that are mitochondrion dependent (Perier et al., 2012). In this study, we identified pathways associated with MPP⁺-induced neuronal cell death in two types of human neuroblastoma cells, that is, SH-SY5Y and SH-EP cells, by using time-course microarray data for the estimation of cell type–specific effects by MPP⁺ treatment. Neuroblastoma cells have been categorized into broad groups that have distinct morphologies and biochemical properties: neuroblast-like N-type cells (e.g., SH-SY5Y), which are clonigenic and tumorigenic and grow rapidly in cultures; S-type cells (e.g., SH-EP), which are not; and I-type cells, which have characteristics of both N- and S-types (Do et al., 2011).

Biological pathway databases such as Kyoto Encyclopedia of Genes and Genomes (Kanehisa et al., 2012) and Reactome (Croft et al., 2014) can capture what is already known about the interplay of genes, proteins, and small molecules using a data model that is accessible to computation. Here, we focused on comparison of perturbed pathways between SH-SY5Y and SH-EP cells using the data model of Reactome, whose central class is a Reaction, and subclasses of Reaction depict these core biological events. Reactions are grouped into pathways, which in turn are assembled into a hierarchy of biological processes. To identify perturbed Reactome pathways from time-course microarray data, we employed structural equation model (SEM) with shortest path model (SPM) for each Reactome pathway as our previous approach (Pepe and Do, 2015). However, in this study we added time, an additional variable to the SEM for the consideration of time effect on the pathway. Perturbed pathways in SH-SY5Y and SH-EP cells showed difference. That is, chromatin/nucleosome modification-related pathways showed significant perturbation in SH-SY5Y cells, while metabolism disorder-related pathways showed significant perturbation in SH-EP cells. This suggests that neuronal cell death by MPP⁺ treatment may occur via various down-stream processes, which might be dependent on cell type or condition.

2. Methods

2.1. Microarray data

Microarray data were obtained from our previous reports (Kim et al., 2013; Do, 2014) and included whole genome expression level at eight time points (0 [control] and 0.5, 1.5, 3, 6, 9, 12, and 24 h) after MPP⁺ treatment for human neuroblastoma cells SH-SY5Y and SH-EP. The number of replicated microarray data at each time point was three and two for SH-SY5Y and SH-EP cells, respectively. After log2 transformation, the expression data were normalized using quantile normalization. After filtering redundant probes, the number of common probes in two groups of experimental data, SH-SY5Y and SH-EP, was 30,483 and had one-to-one mapping of Entrez gene identifier.

2.2. Identification of differentially expressed genes from time-series microarray data

To find differentially expressed genes (DEGs) from time-series microarray data, a Bioconductor R package called maSigPro was employed (Conesa et al., 2006). This program uses a two-step regression approach in which experimental groups are defined by dummy variables. The first regression adjusts a global model and serves to select DEGs, while in the second step a variable selection strategy is applied to identify significant profile differences between experimental groups. In this study, there are no time-series expression data of control group, and we considered each set of time series data separately and found out the genes that show significant regression coefficient only for the variable time with third degree.

2.3. Identification of enriched pathways for significant genes and construction of their shortest path models

Enriched pathways for the set of DEGs during time-course were examined with a Bioconductor R package called ReactomePA, which allows enrichment analysis based on the Reactome pathways (Croft et al., 2014). After identification of enriched pathways in the set of DEGs, the connected structures of DEGs within them were surveyed by considering each pathway as a directed graph. All the shortest paths among differential genes for each enriched pathway were detected with the function get.shortest.paths() of the R package igraph and fused to construct the SPM for the corresponding pathway, as used in our previous work (Pepe and Do, 2015).

2.4. Measurement of modularity for shortest path model

The modularity for each SPM was estimated by the function modularity() of the R package igraph. When the division of the shortest pathway model produces a modularity value greater than 0.4, the community in the SPM was surveyed by using the walktrap.community() function of the R package igraph, and each community or module detected in the shortest model is tested separately.

2.5. Estimation of shortest path model/module showing perturbation during time-course

We evaluated SPM/module showing significant perturbation during time-course by using SEM. For this test, the first is to specify a model describing the causal relationships among variables. In general, all the variables can be represented as a system of structural linear equations as showed below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Y_i = \sum \nolimits_{j \in pa ( i ) } \beta_{ij} Y_j + U_i \quad i \in V\end{align*} \end{document}

and a covariance structure: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \rm cov} ( U_i; U_j ) = \begin{cases} \psi_{ij} \quad \ { \rm if} \ \ i = j \ \ { \rm or} \ \ j \in sib ( i ) \\ \ \ \ 0 \quad { \rm otherwise}\end{cases}\end{align*} \end{document}

The system of linear equations considers every ith node (Y_i), belonging to the set of variables V, as characterized by uni-directed relationships with its “parents” pa(i), defined by path coefficients (β_ij). The covariance structure describes the bi-directed relationships between the unmeasured component of the ith node (U_i) and its “siblings” sib(i), quantified by their covariance (ψ_ij). For each model we fitted two systems of linear equations, one in which we incorporated time and the other we did not. The system of linear equations with time is showed below: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Y_i = \sum \nolimits_{j \in pa ( i ) } \beta_{ij}Y_j + \alpha_iT + U_i\end{align*} \end{document}

where the only difference in the equation, compared to the first system, is that we add the influence of time for each gene in the model. α_i represent the influence of time on the node Y_i. For the estimation of the parameters θ = (β_ij; ψ_ij, α_i) the maximum likelihood estimation was used by the R package lavaan (Rosseel, 2012). To find the significant module/pathway models, the following omnibus test was performed: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}H0: \sum \nolimits_1 ( \theta ) = \sum \nolimits_2 ( \theta ) \ \ { \rm vs}. \ \ H1: \sum \nolimits_1 ( \theta ) \neq \sum \nolimits_2 ( \theta )\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sum \nolimits_1 ( \theta )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sum \nolimits_2 ( \theta )$$ \end{document} are the model-implied covariance matrices of samples with and without time effect, and θ represents the parameter of the model. If the number of genes in the set is greater than the number of samples, the shrinkage covariance (Schäfer and Strimmer, 2005) is applied to estimate the sample covariances given by R package parcor. The statistical significance is determined by comparing likelihood-ratio-test chi-square (χ²) values at a given degree of freedom. If there was a significant difference (p-value <0.05) in the chi-squared goodness-of-fit index, the module/pathway model is considered statistically significant.

2.6. Disease ontology enrichment analysis

To discover disease associations of significant genes, disease ontology (DO) enrichment analysis was carried out with the R package DOSE (Yu et al., 2015), which provides semantic similarity computations among DO terms and genes, and allows biologists to explore the similarities of diseases and of gene functions in disease perspective.

3. Results and Discussion

3.1. Identification of differentially expressed genes in two sets of time-course microarray data

The time-course microarray data for SH-SY5Y and SH-EP cells were obtained from our previous reports (Kim et al., 2013; Do, 2014). The microarray data consisted of eight time points, 0, 0.5, 1.5, 3, 6, 9, 12, and 24 h, after MPP⁺ treatment. The replicate number at each time point was three and two, respectively, for SH-SY5Y and SH-EP cells. To identify DEGs from time-course microarray data, an R package called maSigPro (Conesa et al., 2006) was employed. Before use of the maSigPro program, time-course microarray data were log2 transformed and normalized. After filtration of redundant probes, the common probes between two groups of time-course data (SH-SY5Y and SH-EP cells) were considered for further analysis. The total number of common probe was 30,483 and showed one-to-one mapping gene symbol and Entrez gene identifier. Since both data sets had no time series expression data for control groups, that is, SH-SY5Y and SH-EP cells without MPP⁺ treatment, MaSigPro was applied separately for each time-course data of SH-SY5Y and SH-EP cells. That is, we applied MaSigPro for single group of time-course data and found genes that show significant regression coefficient only for the variable time with third degree. If the regression coefficients of a gene are significant different from 0, the gene is considered as DEG.

The number of DEGs was 406 for SH-SY5Y, and 3,075 for SH-EP cells. The percentage of DEGs in the total genes was 1.3% and 10%, respectively, for SH-SY5Y and SH-EP cells. This suggests that SH-EP cells might be more reactive to MPP⁺ treatment than SH-SY5Y cells. The number of common DEGs between SH-EP and SH-SY5Y cells was 172 (Fig. 1). The common DEGs might be related to core mechanism of MPP⁺-induced neuronal cell death. Pearson correlation coefficients for common DEGs between gene expression profiles of SH-SY5Y and SH-EP cells were estimated and 66 genes with high coefficient values ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mid \rho \mid > 0.8$$ \end{document} ) were listed in Table 1. Except for 3 genes, C1orf21 (chromosome 1 open reading frame 21), ETS2 (v-ets erythroblastosis virus E26 oncogene homolog 2 [avian]), and RAB31 (RAB31, member RAS oncogene family[RAB31]), the remaining 63 genes showed positive correlation. This indicates that the regulation in most of common DEGs might be similar in both types of cells and might participate in the essential mechanism for MPP⁺-induced cell death.

FIG. 1.

DEGs in SH-SY5Y and SH-EP cells. The number within circle represents the number of DEGs and the number outside of two circles represents the number of non-DEGs. DEGs, differentially expressed genes.

Table 1.

Differentially Expressed Genes with High Correlation Coefficient ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mid \rho \mid > 0.8$$ \end{document} ) Between SH-SY5Y and SH-EP cells

Gene	Description	Cor	Gene	Description	Cor
ACTB	Actin, beta (ACTB)	0.88	PPFIA1	Protein tyrosine phosphatase, receptor type, f polypeptide (PTPRF), interacting protein (liprin), alpha 1 (PPFIA1), transcript variant 2	0.81
AMT	Aminomethyltransferase (AMT)	0.89	QPRT	Quinolinate phosphoribosyltransferase (QPRT)	0.88
ANKRD6	Ankyrin repeat domain 6 (ANKRD6)	0.92	R3HCC1	PREDICTED: R3H domain and coiled-coil containing 1 (R3HCC1)	0.87
ARPM1	Actin-related protein M1 (ARPM1)	0.91	RBM12	RNA binding motif protein 12 (RBM12), transcript variant 1	0.96
ATP6V1E2	ATPase, H+ transporting, lysosomal 31kDa, V1 subunit E2 (ATP6V1E2)	0.88	RBM3	RNA binding motif (RNP1, RRM) protein 3 (RBM3), transcript variant 2	0.87
C12orf26	Chromosome 12 open reading frame 26 (C12orf26)	0.83	RFNG	RFNG O-fucosylpeptide 3-beta-N-acetylglucosaminyltransferase (RFNG)	0.89
C18orf56	Chromosome 18 open reading frame 56 (C18orf56)	0.90	RILPL1	Rab interacting lysosomal protein-like 1 (RILPL1)	0.93
C1orf97	Chromosome 1 open reading frame 97 (C1orf97)	0.91	RPL37	Ribosomal protein L37 (RPL37)	0.96
C20orf199	Chromosome 20 open reading frame 199 (C20orf199), transcript variant 3, noncoding RNA	0.82	RRAGB	Ras-related GTP binding B (RRAGB), transcript variant RAGBl	0.89
C7orf53	Chromosome 7 open reading frame 53 (C7orf53)	0.85	SERPINB1	Serpin peptidase inhibitor, clade B (ovalbumin), member 1 (SERPINB1)	0.85
COASY	Coenzyme A synthase (COASY), nuclear gene encoding mitochondrial protein, transcript variant 1	0.91	SNHG7	Small nucleolar RNA host gene 7 (nonprotein coding) (SNHG7), transcript variant 2, noncoding RNA	0.97
DNAJB2	DnaJ (Hsp40) homolog, subfamily B, member 2 (DNAJB2), transcript variant 2	0.93	SNORA45	Small nucleolar RNA, H/ACA box 45 (SNORA45), small nucleolar RNA	0.87
DUT	Deoxyuridine triphosphatase (DUT), nuclear gene encoding mitochondrial protein, transcript variant 1	0.97	SOD2	Superoxide dismutase 2, mitochondrial (SOD2), nuclear gene encoding mitochondrial protein, transcript variant 1	0.94
ELAC1	elaC homolog 1 (E. coli) (ELAC1)	0.80	SPATA20	Spermatogenesis associated 20 (SPATA20)	0.93
FAM24B	Family with sequence similarity 24, member B (FAM24B)	0.90	ST7L	Suppression of tumorigenicity 7 like (ST7L), transcript variant 4	0.91
GRK5	G protein-coupled receptor kinase 5 (GRK5)	0.91	SUOX	Sulfite oxidase (SUOX), nuclear gene encoding mitochondrial protein, transcript variant 1	0.82
GSDMB	Gasdermin B (GSDMB), transcript variant 1	0.85	TGOLN2	Trans-golgi network protein 2 (TGOLN2)	0.96
GYG2	Glycogenin 2 (GYG2)	0.86	TMEM108	Transmembrane protein 108 (TMEM108)	0.87
H3F3B	H3 histone, family 3B (H3.3B) (H3F3B)	0.94	TMEM70	Transmembrane protein 70 (TMEM70), transcript variant 1	0.96
HIST1H2BD	Histone cluster 1, H2bd (HIST1H2BD), transcript variant 2	0.93	TMEM79	Transmembrane protein 79 (TMEM79)	0.99
HIST2H2AA4	Histone cluster 2, H2aa4 (HIST2H2AA4)	0.96	TNFRSF21	Tumor necrosis factor receptor superfamily, member 21 (TNFRSF21)	0.96
HOXB6	Homeo box B6 (HOXB6), transcript variant 3	0.94	TNFSF4	Tumor necrosis factor (ligand) superfamily, member 4 (tax-transcriptionally activated glycoprotein 1, 34kDa) (TNFSF4)	0.92
ICA1	Islet cell autoantigen 1, 69kDa (ICA1), transcript variant 3	0.88	VWCE	Von Willebrand factor C and EGF domains (VWCE)	0.90
IFRD1	Interferon-related developmental regulator 1 (IFRD1), transcript variant 1	0.91	ZNF266	Zinc finger protein 266 (ZNF266)	0.80
LOC100134703	PREDICTED: misc_RNA (LOC100134703), partial miscRNA	0.93	ZNF274	Zinc finger protein 274 (ZNF274), transcript variant ZNF274c	0.94
LOC253039	Hypothetical LOC253039 (LOC253039), noncoding RNA	0.95	ZNF385D	Zinc finger protein 385D (ZNF385D)	0.83
LOC401431	Hypothetical gene LOC401431 (LOC401431)	0.94	ZNF397	Zinc finger protein 397 (ZNF397), transcript variant 2	0.91
LOC729559	PREDICTED: similar to ciliary rootlet coiled-coil, rootletin (LOC729559)	0.90	ZNF79	Zinc finger protein 79 (ZNF79)	0.86
MGC39900	PREDICTED: hypothetical protein MGC39900 (MGC39900)	0.96	AGAP1	ArfGAP with GTPase domain, ankyrin repeat and PH domain 1 (AGAP1), transcript variant 1	0.93
MIR1974	MicroRNA 1974 (MIR1974), microRNA	0.99	FLJ42627	Hypothetical LOC645644 (FLJ42627), noncoding RNA	0.92
MRP63	Mitochondrial ribosomal protein 63 (MRP63), nuclear gene encoding mitochondrial protein	0.86	C1orf21	Chromosome 1 open reading frame 21 (C1orf21)	−0.91
MRPS33	Mitochondrial ribosomal protein S33 (MRPS33), nuclear gene encoding mitochondrial protein, transcript variant 1	0.82	ETS2	v-ets erythroblastosis virus E26 oncogene homolog 2 (avian) (ETS2)	−0.82
PDCL3	Phosducin-like 3 (PDCL3)	0.92	RAB31	RAB31, member RAS oncogene family (RAB31)	−0.94

Cor, Pearson correlation coefficient.

3.2. Enrichment pathway analysis based on Reactome database

Based on the assumption that genes involved in the same biological processes, functions, or localizations present correlated behaviors in gene expression levels, signal intensities, allele occurrences, and so on, statistical tests can be applied to find perturbed or enriched pathways (Hung, 2013). In the case of microarray experiment that changes the expression pattern of thousands of genes, it may only affect the expression patterns of a small handful of biological pathways. Therefore, the identification of enriched pathways in DEGs would be helpful for understanding which molecular mechanisms were evolved in gene expression changes. Here, an R package called reactomePA was employed for the identification of enriched Reactome pathways in DEGs. To determine whether any Reactome pathways annotate a list of DEGs at frequency greater than that would be expected by chance, rectomePA calculates a p-value using the hypergeometric distribution: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}p = 1 - \mathop \sum \limits_ { i = 0 } ^ { k - 1 } \frac { \left( M \atop i \right) \left( N - M \atop n - i \right) } { \left( N \atop n \right) } \end{align*} \end{document}

In the above equation, N is the total number of genes in the background distribution (all Reactome pathways), M is the number of genes within that distribution annotated to the pathway of interest, n is the size of the list of DEGs, and k is the number of genes within that list that are annotated to the pathway. The p-values were adjusted for multiple comparisons. Table 2 shows enriched Reactome pathways of SH-SY5Y and SH-EP cells when the cutoff of adjusted p-value is 0.05. Except for metabolism, no common enriched pathways were observed between SH-SY5Y and SH-EP cells. That is, SH-SY5Y cells mainly showed enrichment on the pathways related to oxidative stress-induced senescence and chromatin organization, while pathways such as mitotic G1-G1/S phase and protein folding were enriched for SH-EP cells. This indicates that MPP⁺-induced neuronal cell death mechanism might not be unique. That is, the inhibition of mitochondrial complex I inhibition by MPP⁺ might cause various detrimental down-stream processes depending on cell type.

Table 2.

Enriched Reactome Pathways in SH-SY5Y and SH-EP cells

Cell type	Reactome ID	Description	GeneRatio	BgRatio	p.adjust
SH-SY5Y	3214847	HATs acetylate histones	8/128	119/6750	0.0308
	1430728	Metabolism	42/128	1465/6750	0.0308
	2559583	Cellular senescence	8/128	143/6750	0.0360
	2559580	Oxidative stress-induced senescence	6/128	88/6750	0.0360
	3247509	Chromatin modifying enzymes	10/128	216/6750	0.0360
	4839726	Chromatin organization	10/128	216/6750	0.0360
	71291	Metabolism of amino acids and derivatives	9/128	188/6750	0.0368
SH-EP	453279	Mitotic G1-G1/S phases	37/961	126/6750	0.0034
	68962	Activation of the prereplicative complex	14/961	30/6750	0.0046
	389977	Postchaperonin tubulin folding pathway	10/961	18/6750	0.0073
	390466	Chaperonin-mediated protein folding	13/961	29/6750	0.0077
	379724	tRNA aminoacylation	16/961	42/6750	0.0084
	391251	Protein folding	14/961	34/6750	0.0084
	379716	Cytosolic tRNA aminoacylation	11/961	24/6750	0.0107
	389958	Cooperation of prefoldin and TriC/CCT in actin and tubulin folding	12/961	28/6750	0.0107
	69190	DNA strand elongation	13/961	32/6750	0.0107
	1430728	Metabolism	251/961	1465/6750	0.0107
	389960	Formation of tubulin folding intermediates by CCT/TriC	10/961	21/6750	0.0107
	176187	Activation of ATR in response to replication stress	14/961	37/6750	0.0121
	69231	Cyclin D associated events in G1	12/961	30/6750	0.0153
	69236	G1 phase	12/961	30/6750	0.0153
	352230	Amino acid transport across the plasma membrane	12/961	31/6750	0.0191
	389957	Prefoldin mediated transfer of substrate to CCT/TriC	11/961	27/6750	0.0191

BgRatio, ratio of the corresponding pathway gene number to the all reactome pathways' gene number; GeneRatio, ratio of the corresponding pathway gene number within the list of differentially expressed genes (DEGs) to the all reactome pathways' gene number in the list of DEGs.

3.3. Construction of shortest path model for enriched pathways and identification of perturbed regions during time-course with SEM

Although enriched pathways identified in the above section might play a role in MPP⁺-induced neuronal cell death, there is no information on which areas of pathways are actually perturbed over time. To detect perturbed regions of pathways over the time, we employed the model-implied covariance approach by SEM. This analysis consisted of three steps: the first is to construct an SPM for each enriched pathway by using shortest path, which is defined as minimum step from one DEG to the other DEG on the pathway. Thus, an SPM for each pathway was built by the union of all shortest paths connecting the DEGs, and may include DEGs or non-DEGs (Pepe and Do, 2015); the second step is to fit each SPM to two SEMs, one describing the causal relationships among variables (genes) and the other describing the causal relationships among variables (genes) affected by time. When the modularity of SPM was higher than 0.5, communities/modules were detected by the walktrap algorithm of the R package igraph. In this case, the fitting to SEM was carried out at level of each module within the SPM rather than the SPM itself (see Materials and Methods for detail); the third step is to compare the likelihood-ratio-test chi-square (χ²) values of two fitted SEMs. If there was a significant difference (p-value <0.05) in the chi-squared goodness-of-fit index, the module/pathway model was considered statistically significant.

Three SPMs for enrich pathways of SH-SY5Y cells were detected with p-values less than 0.00001. These are relative to the Reactome pathways: (1) cellular senescence; (2) chromatin modifying enzymes; (3) chromatin organization (Fig. 2). In the case of SH-EP cells, there were no significant pathways at SPM level but six significant modules (five with p-values <0.00001 and one with p-value = 0.009) were detected within SPM of the metabolism pathway (Fig. 3). These results show that the perturbed pathway regions in SH-SY5Y and SH-EP cells are different. This suggests that the cell death mechanism by MPP⁺ might depend on cell type. All the genes on the significant pathway/module might play a key role in the process of MPP⁺-induced neuronal cell death. Thus, we defined these genes as significant. The number of significant genes was 27 and 109, respectively, for SH-SY5Y and SH-EP cells.

FIG. 2.

Significant shortest path models for enriched reactome pathways in SH-SY5Y cells. (A) Cellular senescence, (B) chromatin modifying enzymes, and (C) chromatin organization.

FIG. 3.

Six significant modules for the reactome metabolism pathway in SH-EP cells.

To survey disease associations of these significant genes, we carried out DO enrichment analysis with the R package DOSE (Yu et al., 2015). Table 3 shows significantly enriched DO terms for SH-SY5Y and SH-EP cells. DO terms, including cancer and disease of cellular proliferation were zenriched for significant genes for SH-SY5Y cells. These DO terms included five genes related to chromatin modification. Entrez genes, including 10919 (EHMT2), 2648 (KAT2A), 8295 (TRRAP), 5931 (RBBP7), and 11091 (WDR5), are involved in histone modification process within chromatin. The modification of histones is often dynamic and reversible and mediated by two antagonistic sets of enzymatic complexes that attach or remove the corresponding chemical groups in a site-specific manner (Sharma and Taliyan, 2015). The purpose of the chromatin remodeling proteins is to alter the nucleosome architecture such that genes are exposed to or hidden from the transcriptional machinery. Recent investigations suggest that gene expression modulated by histone acetylation might be associated with neurodegenerative processes (Sadri-Vakili and Cha, 2006; Harrison and Dexter, 2013). However, in SH-SY5Y cells the expression values of those genes involved in histone modification process were not changed during time-course, while several histone genes such as HIST1H2BD, HIST1H2BK, HIST1H2AC, HIST2H2BE, and HIST2H2AA4 showed downregulation (Table 4).

Table 3.

Enriched Disease Ontology Terms for the Significant Genes in SH-SY5Y and SH-EP Cells

Cell type	Disease ontology term	p.adjust	Entrez gene ID
SH-SY5Y	Cancer	0.009	6233/7157/7316/1432/10919/5928/11091/2648/8295/5931/60/8473.
	Disease of cellular proliferation	0.009	6233/7157/7316/1432/10919/5928/11091/2648/8295/5931/60/8473.
SH-EP	Inherited metabolic disorder	1.54E-09	1738/3030/1737/1593/594/3425/6448/4669/3423/3074/335/5468/336/3931/5465/2646/2799/3073/2717.
	Lysosomal storage disease	9.42E-09	3030/1593/3425/6448/4669/3423/3074/5468/5465/2646/2799/3073/2717.
	Disease of metabolism	9.42E-09	1738/3030/1737/8802/2746/5019/5106/5009/5160/5091/1593/3033/23479/594/3425/6448/46.69/3423/3074/6256/335/5468/336/3931/5465/2645/2646/2799/3073/2717/2026.
	Lipid storage disease	0.0002	3030/1593/3074/5468/5465/2646/3073/2717.
	Syndrome	0.0039	8801/3030/3418/1376/2746/54675/4967/2806/6901/5009/10166/6390/2271/4833/6448/4669/3423/1387/335/5468/336/3931/406/5465/2645/2646/4928/2799.
	Acquired metabolic disease	0.0073	8802/2746/5019/5106/5009/5160/5091/3033/23479/6256/335/5468/336/3931/5465/2645/2646/2026.
	Carbohydrate metabolism disease	0.0197	2746/5019/5106/5091/3033/6256/335/5468/336/3931/5465/2645/2646/2026.
	Glucose metabolism disease	0.0197	2746/5019/5106/5091/3033/6256/335/5468/336/3931/5465/2645/2646/2026.

Table 4.

Fold Change of Selected Significant Genes at 6, 12, and 24 h After MPP⁺ Treatment in SH-SY5Y and SH-EP Cells

	SH-SY5Y			SH-EP
Gene	6 h	12 h	24 h	6h	12 h	24 h
Chromatin modification
KAT2A	0.89	0.90	0.94	0.91	0.82	0.60
EHMT2	0.90	0.94	0.95	0.77	0.74	0.66
TRRAP	0.94	1.00	1.02	0.82	0.85	0.66
RBBP7	1.03	1.08	1.11	1.26	1.09	1.07
WDR5	1.16	1.08	1.00	1.06	1.13	1.22
Histone
HIST1H2BD	0.57	0.52	0.58	0.70	0.57	0.51
HIST2H3D	1.08	1.07	1.02	1.00	1.02	0.96
HIST1H3H	0.82	0.79	0.83	0.72	0.90	0.64
HIST2H2AA4	0.87	0.63	0.49	0.72	0.56	0.30
HIST1H2BK	0.90	0.63	0.30	0.80	0.52	0.43
HIST2H4B	0.86	0.84	0.89	0.92	1.64	1.12
HIST2H2BE	0.82	0.53	0.46	0.63	0.72	0.31
HIST1H2BG	0.79	0.72	0.73	0.65	0.73	0.44
HIST1H2AC	0.76	0.63	0.68	0.81	1.01	0.50
HIST1H4C	0.92	0.86	0.90	1.58	1.81	1.30
Glucose metabolism
GLUD1	1.05	1.01	1.05	0.90	0.89	1.11
OXCT1	1.23	1.22	1.13	1.13	1.06	1.22
PCK2	0.75	0.93	1.18	1.13	3.22	2.66
PC	1.02	1.10	1.04	0.67	0.63	0.61
HADH	0.95	1.01	1.05	0.93	0.83	0.67
RXRA	0.93	0.92	1.08	0.65	0.40	0.62
APOA1	1.13	1.05	1.12	1.02	1.02	0.95
PPARG	1.09	1.07	1.08	0.85	0.85	1.54
APOA2	1.01	0.94	0.92	1.05	0.95	0.94
LCAT	1.07	1.08	1.29	1.41	1.68	2.71
PPARA	0.98	0.99	1.04	0.96	0.96	1.13
GCK	1.03	1.08	1.13	1.03	1.02	1.05
GCKR	1.03	1.01	1.06	0.98	1.01	1.04
ENO2	0.87	0.96	1.10	0.60	0.54	0.64

In the case of SH-EP cells, downregulation was observed not only in histone genes, including HIST1H3H, HIST2H2AA4, HIST1H2BK, HIST2H2BE, HIST1H2BG, and HIST1H2AC, but also histone modification process-related genes such as EHMT2, KAT2A, and TRRAP. This suggests that MPP⁺ might induce chromatin structural defect rather than modification in SH-SY5Y cells. This might lead to chromatin modification as well as chromatin structural defect. Enriched DO terms in SH-EP cells included many diseases associated to metabolism such as inherited metabolic disorder, acquired metabolic disease, carbohydrate metabolism disease, and glucose metabolism disease. Entrez genes such as 5091 (PC), 3033 (HADH), and 2026 (ENO2), which are involved in carbohydrate/glucose metabolism disease, showed downregulation in SH-EP cells, while expression change of these genes was not observed in SH-SY5Y cells (Table 4). It has been reported that lysosomal storage disorder and dysregulation of glucose metabolism are associated with PD (Shachar et al., 2011). Especially, as dysregulation of glucose metabolism is an early event in sporadic PD (Dunn et al., 2014), the understanding of metabolism changes occurring in SH-EP cells might provide a clue for PD development. Our result suggests that MPP⁺ in SH-EP cells might lead to metabolic disorder or stress as well as changes in chromatin structure.

4. Conclusions

The effects of a parkinsonian neurotoxin, MPP⁺, on two types of human neuroblastoma cells, that is, N-type SH-SY5Y and S-type SH-EP cells, were investigated using time-course microarray data. The number of DEGs in SH-EP cells was ∼7.6 times of those in SH-SY5Y cells. This indicates that SH-SY5Y and SH-EP cells may lead to different downstream responsiveness in MPP⁺ treatment although they share a common upstream molecular event, such as the inhibition of mitochondrial complex I by MPP⁺. Pathway analysis with Reactome clearly showed the difference of perturbed pathways between SH-SY5Y and SH-EP cells. Pathways related to cellular senescence and chromatin modification/organization showed perturbation in SH-SY5Y cells, while pathways related to many metabolisms represented perturbation in SH-EP cells. In SH-SY5Y cells, expression values of genes such as KAT2A, EHMT2, TRRAP, RBBP7, and WDR5, which are involved in chromatin modification, were not changed over time, but histone genes, including HIST1H2BD, HIST1H2BK, HIST1H2AC, HIST2H2BE, and HIST2H2AA4, showed downregulation. This suggests that MPP⁺ in SH-SY5Y cells might induce structural defect rather than modification in nucleosome.

On the other hand, MPP⁺ in SH-EP cells led to perturbation of metabolisms that include Entrez genes such as 5091 (PC), 3033 (HADH), 2026 (ENO2), and 3030 (HADHA). Pyruvate carboxylase (PC) is involved in gluconeogenesis, lipogenesis, insulin secretion, and synthesis of the neurotransmitter glutamate, while hydroxyacyl-CoA dehydrogenase (HADH) and long-chain-3-hydroxyacyl-CoA dehydrogenase (HADAH) are mitochondrial functional proteins. Neural enolase (ENO2) has neurotrophic and neuroprotective properties on a broad spectrum of central nervous system (CNS) neurons. These four metabolic genes showed downregulation. This indicates that metabolic stress or disorder may play a key role in MPP⁺-induced cell death for SH-EP cells. In addition, MPP⁺ might lead to chromatin modification as well as chromatin structural defect in SH-EP cells via expression changes of corresponding genes. This suggests that effect of MPP⁺ on chromatin might play an important role in both cell types: N-type SH-SY5Y and S-type SH-EP cells. Conclusively, the mitochondrial complex I inhibition by MPP⁺ might lead to neuronal cell death via cell type-dependent downstream processes. This might provide a hint of explanation for coexistence of different morphological types of cell death in PD patients. In addition, our results might be useful in the development of multitarget PD drugs because MPP⁺-induced cell death mechanisms in both SH-SY5Y and SH-EP cells show association with pathogenesis of PD.

Footnotes

Acknowledgment

This study was supported by the Basic Science Research Program through the National Research Foundation (NRF) of Korea funded by the Ministry of Education, Science, and Technology (MEST) (No. 2012R1A1A2005622).

Author Disclosure Statement

The authors declare that no financial conflicts of interest exist.

References

Ammal Kaidery

, Tarannum

, and Thomas

2013. Epigenetic landscape of Parkinson's disease: Emerging role in disease mechanisms and therapeutic modalities. Neurotherapeutics. 10, 698–708.

Castro

A.A.

, Ghisoni

, Latini

, et al. 2012. Lithium and valproate prevent olfactory discrimination and short-term memory impairments in the intranasal 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) rat model of Parkinson's disease. Behav. Brain Res. 229, 208–215.

Conesa

, Nueda

M.J.

, Ferrer

, et al. 2006. maSigPro: A method to identify significantly differential expression profiles in time-course microarray experiments. Bioinformatics. 22, 1096–1102.

Croft

, Mundo

A.F.

, Haw

, et al. 2014. The Reactome pathway knowledgebase. Nucleic Acids Res. 42, D472–D477.

Dagda

R.K.

, Das Banerjee

, and Janda

2013. How Parkinsonian toxins dysregulate the authophasy machinery. Int. J. Mol. Sci. 14, 22163–22189.

J.H.

, 2014. Neurotoxin-induced pathway perturbation in human neuroblastoma SH-EP cells. Mol. Cells. 37, 672–684.

J.H.

, Kim

I.S.

, Lee

J.D.

, et al. 2011. Comparison of genomic profiles in human neuroblastic SH-SY5Y and substrate-adherent SH-EP cells using array comparative genomic hybridization. BioChip J. 5, 165–174.

Dunn

, Allen

G.F.

, Mamais

, et al. 2014. Dysregulation of glucose metabolism is an early event in sporadic Parkinson's disease. Neurobiol. Aging., 35, 1111–1115.

Fernandez-Espejo

2004. Pathogenesis of Parkinson's disease: Prospects of neuroprotective and restorative therapies. Mol. Neubiol. 29, 15–30.

10.

Harrison

I.F.

, and Dexter

D.T.

2013. Epigenetic targeting of histone deacetylase: Therapeutic potential in Parkinson's disease?. Pharmacol. Ther. 140, 34–52.

11.

Hung

J.H.

2013. Gene Set/Pathway enrichment analysis. Methods Mol. Biol. 939, 201–213.

12.

Kanehisa

, Goto

, Sato

, et al., 2012. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 40, D109–D114.

13.

Kim

I.S.

, Choi

D.-K.

, and Do

J.H.

2013. Genome-wide temporal responses of human neuroblastoma SH-SY5Y cells MPP+ neurotoxicity. BioChip J. 7, 247–257.

14.

Langston

J.W.

, Ballard

, Tetrud

J.W.

, et al. 1983. Parkinsonism in humans due to a product of meperidine-analog synthesis. Science. 219, 979–980.

15.

Pepe

, and Do

J.H.

2015. Estimation of dysregulated pathway regions in MPP⁺ treated human neuroblastoma SH-EP cells with structural equation model. BioChip J. 9, 131–138.

16.

Perier

, Bové

, and Vila

2012. Mitochondria and programmed cell death in Parkinson's disease: Apoptosis and beyond. Antioxid. Redox Signal. 16, 883–895.

17.

Rosseel

2012. lavaan: An R package for structural equation modeling. J. Stat. Softw. 48, 1–36.

18.

Sadri-Vakili

, and Cha

J.H.

2006. Histone deacetylase inhibitors: A novel therapeutic approach to Huntington's disease (complex mechanism of neuronal death). Curr. Alzheimer Res. 3, 403–408.

19.

Schäfer

, and Strimmer

2005. A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics. Stat. Appl. Genet. Mol. Biol. 4, 32.

20.

Shachar

, Lo Bianco

, Recchia

, et al. 2011. Lysosomal storage disorders and Parkinson's disease: Gaucher disease and beyond. Mov. Disord. 26, 1593–1604.

21.

Sharma

, and Deshmukh

2014. Vinpocetin attenuates MPTP-induced motor deficit and biochemical abnormalities in Wistar rats. Neuroscience. 286, 393–403.

22.

Sharma

, and Taliyan

2015. Targeting histone deacethylases: A novel approach in Parkinson's disease. Parkionsons Dis. 2015, 303294.

23.

, Wang

L.G.

, Yan

G.R.

, et al. 2015. DOSE: An R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. 31, 608–609.