Differentially expressed genes and pathways induced by organophosphates in human neuroblastoma cells

Abstract

Organophosphates (OPs) are toxic chemicals commonly used as pesticides and herbicides. Some OPs are highly toxic to humans and have been used in warfare and terrorist attacks. In order to elucidate the molecular mechanisms of injury caused by OPs, the differentially expressed genes were analyzed in human SK-N-SH neuroblastoma cells induced by three OPs. The SK-N-SH cells were treated with one of the three OPs, chlorpyrifos, dichlorvos or methamidophos at LC20 (high-dose), the concentration causing 20% cell death, as well as 1/20 of LC20 (low-dose), a sub-lethal concentration with no detectable cell death, for 24 h. The genome-wide gene changes were identified by Agilent Microarray System, and analyzed by microarray analysis tools. The analysis revealed neuroblastoma cells treated with the high doses of all three OPs markedly activated cell apoptosis and inhibited cell growth and proliferation genes, which would most likely lead to the process of cell death. Interestingly, the analysis also revealed significant decrease in expressions of many genes in a specific spliceosome pathway in cells treated with the low doses of all three different OPs. The change of spliceosome pathway may represent an important mechanism of injury in neuronal cells exposed to low doses of various OPs. In addition to unraveling a potentially different form of OP pathogenesis, this finding could provide a new diagnostic marker in assessing OP-associated injury in cells or tissues. In addition, these results could also contribute to the development of new prevention and/or therapeutic regimens against OP toxicity.

Keywords

organophosphate neuroblastoma microarray pathway analysis gene regulation

Introduction

Organophosphates (OPs), toxic chemicals that exist in many pesticides, can cause severe chronic and acute illnesses or even death in humans.¹ In the USA, approximately half a million tons of pesticides are used in agriculture and households annually, and thousands of people were admitted to hospital because of OP poisoning, and children are the most common victims.^2,3 Therefore, OPs are a significant environmental threat to human health.

OPs are known to act through their inhibitory effects on the acetylcholineseterase enzyme (AChE) at nerve endings.⁴ The AChE enzyme is critical to normal control of nerve impulse transmission from nerve fibers to smooth and skeletal muscle cells, glandular cells, autonomic ganglia and in the central nervous system (CNS). The loss of AChE function allows accumulation of excess acetylcholine and results in over-stimulation of the effector organs.⁴ High acetylcholine concentrations accumulated in the CNS cause sensory and behavioral disturbances, incoordination, depressed normal motor function and respiratory depression. Increased pulmonary secretions coupled with bronchoconstriction are the usual causes of death from OP poisoning.⁵ On the other hand, it is evidently more difficult to document the damage of peripheral or central nervous system caused by chronic exposure to low doses of OPs.²

There have been many studies using biochemical, physiological, behavioral and histopathological approaches in animal models to search for cell injury mechanisms by organophosphorous agents.^6–9 Most of these studies have looked at specific aspects on OP poisoning, such as changes in phospholipase C activities in rat brain, cholin-esterase inhibition in whole blood and changes in the expression of synaptobrevin/vesicle-associated membrane protein in hippocampal neurons. Despite the general consensus of the disease progression of OP poisoning, few studies have been conducted using a systematic approach to research the detailed molecular mechanisms of injury by organophosphorous agents. Recently, with the development of microarray technology, it has become possible to study the global changes of gene expression profiles before and after exposures to toxic substances in biological systems.¹⁰ However, the few microarray studies that reported changes of the gene expression profiles associated with OP damage in cells^11–14 either focused on examining only a subset of gene expression profiles or examined the effect of one particular OP.

In addition to the common mechanism of inhibiting AChE, OPs could also be affecting many non-cholinergic neurochemical processes as well as non-neurochemical pro-cesses.¹ To develop effective diagnostic means, preventive measures and treatments against OP-associated toxicities and injury, it is critically important to better understand the molecular mechanisms of injury caused by various OPs. In the beginning, we examined genome-wide gene expression profiles in three different human cell lines exposed to two (high and low) doses of three different OPs. The alteration of gene expressions following exposure to all three different OPs for both doses were found to be much more prominent in the human neuronal cell line, SK-N-SH neuroblastoma cells, than in A549 lung alveolar epithelial cells and HepG2 hepatocellular cells. In the present study, we focused on analyzing the differentially expressed genes identified in human SK-N-SH neuroblastoma cells treated with three different OPs at high and low concentrations. The significant genes commonly expressed in the cells treated with all three different OPs may be used as potential candidates for biomarkers in the evaluation of the degree of OP injury.

Materials and methods

Cell culture

Human neuroblastoma SK-N-SH cell line, A549 lung alveolar epithelial cells and HepG2 hepatocellular cells were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). The cells were grown in usually Dulbecco's modified Eagle's medium medium with 10% fetal bovine serum (FBS), 2 mmol/L L-glutamine, 500 μmol/L sodium pyruvate, 50 U/mL penicillin and 50 μg/mL streptomycin. Cells were maintained at 37°C in a humidified atmosphere with 5% CO₂. The culture medium was changed every three days.

Measurement of cell viability and determination of LC20

SK-N-SH cells were seeded on 96-well plates at a density of 2 × 10⁴ cells/well and treated with three OPs: chlorpyrifos, dichlorvos or methamidophos for 24 h. The OPs were dissolved in 2% dimethyl sulfoxide (DMSO) (the final concentration of DMSO in medium is 2%, v/v) and 2% DMSO was used as the control. The LC80 (the concentration causing 80% cell death), starting concentrations of three organophosphorus agents were approximately determined as following: chlorpyrifos 800 μg/mL (2.28 mmol/L), dichlorvos 400 μg/mL (1.81 mmol/L) and methamidophos 4000 μg/mL (28.3 mmol/L). Serial two-fold dilutions were used to determine the LC20 (the concentration causing 20% cell death). Each condition was assessed in triplicate. Cell viability was determined by 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl-2H-tetrazolium bromide-based In Vitro Toxicology Assay Kit (Sigma, St. Louis, MO, USA). The samples were measured at the absorbance wavelength of 540 nm using a microplate reader (PerkinElmer, Waltham, MA, USA). Cell viability was expressed as a percentage of the control viability. The LC20 of each OP reagent was designed as the high-dose, and the 1/20 of the LC20 for each OP as low-dose for the cell treatment.

Treatment of SK-N-SH cells with OPs

Human neuroblastoma SK-N-SH cells were grown in DMEM medium with 10% FBS, and maintained at 37°C in a humidified atmosphere with 5% CO₂. The culture medium was changed every three days. When the cells were reaching 80% confluence, the media were changed to include various concentrations of OPs. For high-dose treatments, LC20 of each OP was determined as following: chlorpyrifos 100 μg/mL (0.29 mmol/L), dichlorvos 25 μg/mL (0.11 mmol/L) and methamidophos 500 μg/mL (3.54 mmol/L)). For low-dose treatments, the 1/20 of the LC20 for each OP was used to treat the cells. Control and treated SK-N-SH cells were incubated for 24 h before processing total RNA isolation.

Total RNA isolation and microarray experiments

Total RNAs from the treated and control cells were isolated by using Trizol (Invitrogen Valenica, CA, USA), and purified by RNeasy kit (Qiagen Valenica, CA, USA) following the manufacturers’ protocols. The quality of isolated RNA was evaluated by Agilent Bioanalyzer 2100 (Clara, CA, USA). RNA samples with the highest integrity and purity were further processed for microarray experiments.

Agilent human whole genomic 4 × 44K arrays were used in this study. Duplicates of the microarray experiments were performed with each dosage treatment of three OPs. The detailed information regarding target preparation, microarray hybridization and image quantification was described in the product protocol (Agilent Technologies). In brief, 200 ng of total RNA was used for complementary RNA (cRNA) amplification following Agilent two-color protocol. Cy3 dye was used to label the cRNA from a control sample, and Cy5 dye was used to label the cRNA from the samples treated with organophosphorus reagents. The final mixture of 825 ng amplified labeled cRNAs in the hybridization cocktail (Agilent) was placed onto the Agilent microarray slide, and hybridized at 42°C in Agilent hybridization oven for 14 h. After washing and spinning dry, all microarray slides were scanned by Agilent G2565CA Microarray Scanner with 2 μm resolution. Scanned raw images were analyzed and data files were generated with Agilent Feature Extraction Software (v10.1).

Data analysis and statistics

Data files were uploaded into the mAdb (a microarray database developed by National Cancer Institute¹⁵) and ArrayTrack system (a microarray analytic tool developed by Food and Drug Adminstration¹⁶). The advanced filters, minimum fluorescent intensity of 100 in both Cy3 and Cy5 channels, were applied before data analysis. All of the extracted data were normalized using a Lowess normalization method, analyzed by the software tools available through the database.

Statistic analyses were performed, and t-test and analysis of variance were used to compare the differences of gene expression levels between the cells treated with different organophosphorus reagents. Differentially expressed genes were selected based on P value and fold changes with the following criteria: P < 0.05 and fold change ≥2.

Significance Analysis of Microarray (SAM) provided by mAdb was applied to identify commonly differentially expressed genes induced by various OPs. Genes with missing values were removed from the SAM analysis. The number of permutations was set to be 100, and the maximum value of false discovery rate (FDR) was set to be < 0.025 as threshold criteria for counting the numbers of significant genes. Gene ontology (GO) and pathway analysis on differentially expressed genes identified by SAM were performed using the mAdb.

Differentially expressed genes induced by various OP reagents were exported into the Ingenuity Pathway Analysis (IPA) System (www.ingenuity.com). Biological functions and relevant pathways were constructed by the web software provided by IPA. The relationships between genes in the networks were examined with GO analysis to find the related functions. Functional analysis and canonical pathway analysis identified the functions and the pathways from the Ingenuity database that were most significant to the data sets.

Gene-specific reverse-transcription polymerase chain reaction and realtime quantitative PCR

Reverse-transcription polymerase chain reaction (RT-PCR) was performed using gene-specific primer pairs to verify the changes of expression in genes identified by microarray technology. The primers used in this study were shown inTable 1. Briefly, the first-strand cDNA was synthesized by using Superscript II Reverse Transcriptase (Invitrogen) in a 20 μL final volume reaction with 1 μg of total RNA and an oligo-dT primer. For PCR reactions, the resulting first-strand cDNA was amplified by using gene specific primers with AmpliTaq Gold (Applied Biosystems, Roche, Foster City, CA, USA). The RT-PCR reactions were carried out by initial activation at 94°C for five minutes; followed by 35 cycles of the three-step cycling at 94°C denaturation for 45 s, at 55°C annealing for 30 s and at 72°C extension for 45 s. The PCR products were detected by 2.0% agarose gels and visualized by ethidium bromide gel staining. A realtime quantitative PCR (qPCR) was used to confirm the expression level changes of the target significant genes by quantitatively assessing the amounts of the specific mRNA produced in the cells treated with OPs. The mRNA examined was similarly converted into cDNA templates and primer pairs used in RT-PCR were applied to conduct qPCR. The qPCR reactions were carried out by initial activation at 95°C for 10 min; followed by 40 cycles of two-step cycling at 95°C for 30 s and at 60°C for 60 s. The changes of cDNA level were expressed as log2 ratios in cells with high-dose OP treatments relative to the non-treated control.

Table 1

Primers used in reverse transcription polymerase chain reaction

Primer name	Primer sequence
MSH2-F	AGTGTTAACTGTCAGTGCCCAT
MSH2-R	GCAAACAGTCCTCAGTTACAGC
HPS3-F	AGAAGGATCAATGGGAAGGCA
HPS3-R	TGATATCACCGCTGTTAGGGT
CYP2W1-F	AACACTTCCTGCAGCTGTGC
CYP2W1-R	GCCAGTGCCTGTCTCTGAAG
GAPDH-F	AGAGAGACCCTCACTGCTGG
GAPDH-R	GTACATGACAAGGTGCGGCT
TNF-F	TCCGTGAAAACGGAGCTGAAC
TNF-R	AGGCTCAGCAATGAGTGACA

‘F’ represents forward primer and ‘R’ represents reverse primer

Terminal deoxynucleotidyl transferase-mediated dUTP nick-end labeling assay

The terminal deoxynucleotidyl transferase-mediated dUTP nick-end labeling (TUNEL) method was performed to label 3’-end of fragmented DNA of the apoptotic SK-N-SH cells in cultures treated with three different OPs at the high-dose for 24 h. The assay was conducted using Deadend Fluorometric TUNEL system (Promega, Madison, WI, USA). The cell staining procedure was performed following the product instruction. The images of stained cells were captured on a Zeiss LSM-5 (Cambridge, UK) confocal microscope. Furthermore, 100 cells for each sample were examined to determine the percentage of apoptotic cells, where green fluorescent cells were counted as the cells that were undergoing apoptosis.

Results

Cell viability and determination of LC20 and sublethal concentration

The pilot experiments were performed to determine the approximate concentrations causing 80% cell death (LC80) for three OPs: chlorpyrifos, dichlorvos and methamidophos. The LC80 were used as the starting concentrations for determining the LC20 (the concentration causing 20% cell death) for three OPs. Two-fold serial dilutions of three OPs were made from the starting concentrations (LC80), and were used to treat the cells. The survival curves for SK-N-SH neuroblastoma cells treated with three different OPs in serial dilutions were shown in Figure 1. The survival curves for the three OPs not only indicate the toxicity differences among three OPs, but also reflect a dose-dependent variation of the toxicity for each OP. The LC20 for the three OPs was determined as: chlorpyrifos 100 μg/mL (0.29 mmol/L), dichlorvos 25 μg/mL (0.11 mmol/L) and methamidophos 500 μg/mL (3.54 mmol/L), respectively. Subsequently, 1/20 of LC20 was designed as the sublethal concentration at which no obvious cell death was detected in the culture (data not shown). The LC20 so determined was used in the high-dose treatments, and the sublethal concentration was used in low-dose treatments in the present OP injury study.

Figure 1

The survival curve for neuronal cells treated with three different OPs. SK-N-SH cells were treated with three OPs: chlorpyrifos, dichlorvos and methamidophos, respectively. The starting concentrations were the approximate LC80 and two-fold serial dilutions were used to determine the LC20. The average and standard deviation for triplicate sample treatments are shown on each spot. OP, organophosphate. (A color version of this figure is available in the online journal)

Identification of common significant genes and pathways induced by high-dose OPs in neuroblastoma cells

The gene expression profiles of SK-N-SH neuroblastoma cells treated with high-dose (LC20) OPs were analyzed by microarray technology. Duplicate microarray experiments were performed for each OP treatment. In the beginning, we analyzed the gene expression profiles with ArrayTrack. The significant genes were identified by the cut-offs of P < 0.05 and fold change >1.5 (the mean intensity is greater than 100). The identification of significant genes was based on Welch's t-test. As shown in Figure 2, 1342 genes were identified significantly expressed in SK-N-SH cells treated with chlorpyrifos; 1768 genes for dichlorvos and 1033 genes for methamidophos. However, only 72 genes (less than 10%) were found common to all the three high-dose OP treatments. The small number of common significant genes limited the pathway analysis and function study. As our focus was on the identification of the common significant genes induced by high-dose treatments with the three different OPs, we put all six high-dose data-sets together and analyzed with mAdb using SAM with FDR <0.025. A total of 496 common significant genes were differentially expressed in response to treatments of three OPs: 175 upregulated genes and 321 downregulated genes (Figure 3). The selected gene list (Table 2) shows all the common significant genes induced by high-dose treatments which were discussed in this article. Many of the upregulated genes play important roles in cell apoptotic process and cell cycle arrest, such as ERN1 (endoplasmic reticulum to nucleus signaling 1), PRNP (prion protein), NET1 (neuro-epithelial cell transforming 1), CLN3 (ceroid-lipofuscinosis, neuronal 3), and MT1A and MT3 (metallothionein 1A and 3).

Table 2

Selected common significant genes in SK-N-SH cells treated with three high-dose organophosphates identified by ArrayTrack and mAdb analysis

RefSeq	Gene symbol	Annotations	Fold^*
Upregulated genes
NM_001433	ERN1	Endoplasmic reticulum to nucleus signaling 1	4
NM_000311	PRNP	Prion protein (PRNP)	2.3
NM_017781	CYP2W1	Cytochrome P450, family 2, subfamily W, polypeptide 1	2.3
NM_032383	HPS3	Hermansky–Pudlak syndrome 3	2.2
NM_000594	TNF	Tumor necrosis factor (TNF superfamily, member 2)	2.0
NM_005863	NET1	Neuroepithelial cell transforming 1 (NET1),	1.7
NM_005946	MT1A	Metallothionein 1A (MT1A)	1.6
NM_000086	CLN3	Ceroid-lipofuscinosis, neuronal 3 (CLN3),	1.4
NM_005954	MT3	Metallothionein 3 (MT3)	1.3
Downregulated genes
NM_002467	MYC	v-myc myelocytomatosis viral oncogene homolog	-1.0
NM_001256	CDC27	Cell division cycle 27 homolog	-1.2
NM_000251	MSH2	mutS homolog 2, colon cancer, nonpolyposis type 1	-1.3
NM_001786	CDK1	Cyclin-dependent kinase 1 (CDK1),	-1.4
NM_002643	PIGF	Phosphatidylinositol glycan anchor biosynthesis, class F	-1.4
NM_001880	ATF2	Activating transcription factor 2	-1.4
NM_024570	RNASEH2B	Ribonuclease H2, subunit B (RNASEH2B),	-1.4
NM_020686	ABAT	4-Aminobutyrate aminotransferase	-1.5
NM_003659	AGPS	Alkylglycerone phosphate synthase (AGPS)	-1.6
NM_001789	CDC25A	Cell division cycle 25 homolog A	-1.7
NM_006030	CACNA2D2	Ccalcium channel, voltage-dependent, alpha 2/delta subunit 2	-1.7
NM_020791	TAOK1	TAO kinase 1 (TAOK1)	-1.7
NM_000946	PRIM1	Primase, DNA, polypeptide 1 (49 kDa)	-2.0

Significant genes were selected based on P value and fold changes with the following criteria: P < 0.05 and fold change ≥2

Log2 fold ratio: 1 equals two-fold difference; 2 equals four-fold difference, and so on

Figure 2

Common significant genes identified in SK-N-SH cells using ArrayTrack. Cells were treated with high-dose of three OPs (Chlorpyrifos, Dichlorvos and Methamidophos) and the significant genes were extracted by ArrayTrack. The common significant genes for high-dose (72) were the overlap among all the three different OP treatments. (A color version of this figure is available in the online journal)

Figure 3

Common significant genes identified in SK-N-SH cells using mAdb. Cells were treated with high or low dose of three different OPs (duplicate for each OP) and the significant genes were extracted by mAdb. Six arrays for high-dose treatments and six arrays for low-dose treatments were used for the analysis. The common significant genes for both high-dose (496) and low-dose (668) treatments were further analyzed to find the overlapping common significant genes (153). The upregulated (Up) and downregulated (Down) genes were also shown for high- and low-dose treatments. OP, organophosphate. (A color version of this figure is available in the online journal)

In addition to upregulated genes, many genes implicated in cell cycle pathway, DNA replication and metabolism were downregulated.Table 3 shows the major pathways obtained by KEGG pathway analysis performed on the common differentially expressed genes induced by high-dose OPs. KEGG pathway analysis indicated inhibitory effects of OP treatments on many pathways, including metabolic pathway (AGPS, ABAT, PIGF, etc), cell cycle pathway (CDC25A, CDK1, CDC27, PRIM1, RNASEH2B, MYC, etc) and MAPK (mitogen-activated protein kinase) signaling pathway (TAOK1, CACNA2D2, ATF2, etc), as shown in Tables 3 –5.

Table 3

Major pathways obtained from KEGG analysis of significant genes induced by high- or low-dose OP treatments

Pathway name	Genes involved
High-dose treatment-induced pathway
Metabolic pathway	39
Cell cycle	10
MAPK signaling pathway	9
p53 signal pathway	9
Low-dose treatment-induced pathway
Metabolic pathway	42
Spliceosome pathway	18
MAPK signaling pathway	14
Cell cycle	12

OP, organophosphate; MAPK, mitogen-activated protein kinase

Table 4

Significant genes in cell cycle pathway induced by high- or low-dose organophosphate treatments obtained from KEGG analysis

RefSeq	Gene symbol	Annotations	Fold^*
High-dose cell cycle pathway
NM_001789	CDC25A	Cell division cycle 25 homolog A	-1.7
NM_006142	SFN	Stratifin	-1.7
NM_001786	CDK1	Cyclin-dependent kinase 1	-1.4
NM_001256	CDC27	Cell division cycle 27 homolog, transcript variant 2	-1.2
NM_078626	CDKN2C	Cyclin-dependent kinase inhibitor 2C (p18, inhibits CDK4)	-1.2
NM_006716	DBF4	DBF4 homolog	-1.1
NM_006607	PTTG2	Pituitary tumor-transforming 2	-1.0
NM_031966	CCNB1	Cyclin B1	-1.0
NM_002467	MYC	v-myc myelocytomatosis viral oncogene homolog	-1.0
NM_001237	CCNA2	Cyclin A2	-1.0
Low-dose cell cycle pathway
NM_005445	SMC3	Structural maintenance of chromosomes 3	-1.4
NM_001007793	BUB3	Budding uninhibited by benzimidazoles 3 homolog	-1.2
NM_078626	CDKN2C	Cyclin-dependent kinase inhibitor 2C (p18, inhibits CDK4)	-1.2
NM_022662	ANAPC1	Anaphase promoting complex subunit 1	-1.1
NM_012381	ORC3	Origin recognition complex, subunit 3	-1.1
	YWHAZ	Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein	-1.1
NM_002467	MYC	v-myc myelocytomatosis viral oncogene homolog	-1.1
NM_001527	HDAC2	Histone deacetylase 2	-1.0
NM_004964	HDAC1	Histone deacetylase 1	-1.0
NM_001146156	GSK3B	Glycogen synthase kinase 3 beta	-1.0
NM_031966	CCNB1	Cyclin B1	-1.0
NM_053056	CCND1	Cyclin D1	-1.0

Significant genes were selected based on P value and fold changes with the following criteria: P < 0.05 and fold change ≥2

Log2 fold ratio: 1 equals two-fold difference; 2 equals four-fold difference, and so on

Table 5

Significant genes in MAPK pathway induced by high- or low-dose OP treatments obtained from KEGG analysis

RefSeq	Gene symbol	Annotations	Fold^*
Low-dose MAPK pathway
NM_020791	TAOK1	TAO kinase 1	-1.7
NM_006030	CACNA2D2	Calcium channel, voltage-dependent, alpha 2/delta subunit 2	-1.7
	ATF2	Activating transcription factor 2	-1.4
NR_038322	MAP3K13	Mitogen-activated protein kinase kinase kinase 13	-1.2
NM_002467	MYC	v-myc myelocytomatosis viral oncogene homolog	-1.0
NM_203401	STMN1	Stathmin 1	-1.0
NM_031896	CACNG7	Calcium channel, voltage-dependent, gamma subunit 7	1.1
NM_005247	FGF3	Fibroblast growth factor 3	1.2
NM_001205294	CACNA1E	Calcium channel, voltage-dependent, R type, alpha 1E subunit	1.4
Low-dose MAPK pathway
NM_002577	PAK2	p21 protein (Cdc42/Rac)-activated kinase 2	-1.5
NM_020791	TAOK1	TAO kinase 1	-1.4
NM_033360	KRAS	v-Ki-ras2 Kirsten rat sarcoma viral oncogene homolog	-1.3
NM_006030	CACNA2D2	Calcium channel, voltage-dependent, alpha 2/delta subunit 2	-1.3
NM_030662	MAP2K2	Mitogen-activated protein kinase kinase 2	-1.3
NR_038322	MAP3K13	Mitogen-activated protein kinase kinase kinase 13	-1.2
NM_199054	MKNK2	MAP kinase interacting serine/threonine kinase 2	-1.2
NM_030662	MAP2K2	Mitogen-activated protein kinase kinase 2	-1.1
NM_199004	ARRB2	arrestin, beta 2	-1.1
NM_002467	MYC	v-myc myelocytomatosis viral oncogene homolog	-1.1
NM_005343	HRAS	v-Ha-ras Harvey rat sarcoma viral oncogene homolog	-1.0
NM_182810	ATF4	Activating transcription factor 4	-1.0
NM_001946	DUSP6	Dual specificity phosphatase 6	-1.0
	MAPK8	Mitogen-activated protein kinase 8	-1.0

OP, organophosphate; MAPK, mitogen-activated protein kinase

Significant genes were selected based on P value and fold changes with the following criteria: P < 0.05 and fold change ≥2

Log2 fold ratio: 1 equals two-fold difference; 2 equals four-fold difference, and so on

Identification of common significant genes and pathways induced by low-dose OPs in neuroblastoma cells

We also identified the common significant genes induced by low-dose treatments with three different OPs, using SAM with FDR <0.025. A total of 668 common significant genes were differentially expressed in response to low-dose treatments: 58 upregulated genes and 610 downregulated genes (Figure 3). As shown in Figure 3, there were 153 common significant genes between high- and low-dose treatments.Table 3 shows the major pathways induced by the low-dose treatments with three different OPs. Three of the four major pathways were common between high- and low-dose treatments: metabolic pathway; cell cycle pathway and MAPK signaling pathway. Interestingly, in addition to the common pathways induced by both high- and low-dose treatments of OPs, the low-dose treatments of OPs induced strong inhibitory effects on the genes involved in spliceosome pathway. A total of 18 genes involved in the spliceosome pathway was downregulated by low-dose OPs (Table 6). Figure 4 shows the standard spliceosome pathway, which is made up of five small nuclear ribonucleoproteins (snRNPs), U1,U2, U4, U5 and U6 snRNPs, and several spliceosome-associated proteins. SNRNP70 encoded U1 snRNP, SF3A3 encoded U2 snRNP, LSM encoded U4/U6 snRNP and PRPF8 encoded U5 snRNP were down-regulated by low-dose treatments. Other snRNP-associated proteins, such as HNRNPK, HNRNPA1, NCBP2, CWC15, SRSF1, RBM8A, SART1, etc, were also downregulated (Table 6).

Table 6

Genes involved in spliceosome pathway obtained from KEGG analysis of significant genes induced by low-dose treatments with three different organophosphates

RefSeq	Gene symbol	Annotation	Fold^*
NM_007362	NCBP2	Nuclear cap-binding protein subunit 2	-1.5
NM_002136	HNRNPA1	Heterogeneous nuclear ribonucleoprotein A1	-1.4
NM_005105	RBM8A	RNA-binding motif protein 8A	-1.3
NM_005146	SART1	Squamous cell carcinoma antigen recognized by T	-1.1
NM_006802	SF3A3	Splicing factor 3a, subunit 3, 60 kDa	-1.1
NM_021239	RBM25	RNA-binding motif protein 25	-1.1
NM_002139	RBMX	RNA-binding motif protein, X-linked	-1.1
NM_032864	PRPF38A	PRP38 pre-mRNA processing factor 38	-1.1
NM_016403	CWC15	CWC15 spliceosome-associated protein	-1.1
NM_006445	PRPF8	PRP8 pre-mRNA processing factor 8	-1.1
NM_002370	MAGOH	Mago-nashi homolog, proliferation-associated	-1.1
NM_014463	LSM3	LSM3 homolog, U6 small nuclear RNA associated	-1.1
NM_006924	SRSF1	Serine/arginine-rich splicing factor 1	-1.0
NM_006590	USP39	Ubiquitin specific peptidase 39	-1.0
NM_002140	HNRNPK	Heterogeneous nuclear ribonucleoprotein K	-1.0
NM_003089	SNRNP70	Small nuclear ribonucleoprotein 70kDa (U1)	-1.0
NM_019610	RBMXL1	RNA-binding motif protein, X-linked-like 1	-1.0
NM_018048	MAGOHB	Mago-nashi homolog	-1.0

Significant genes were selected based on P value and fold changes with the following criteria: P< 0.05 and fold change ≥2

Log2 fold ratio: 1 equals two-fold difference; 2 equals four-fold difference, and so on

Figure 4

Representative spliceosome pathway. One of the most significant pathways identified by KEGG analysis of significant genes induced by low-dose treatments was shown here. (A color version of this figure is available in the online journal)

Verification of gene changes by RT-PCR and qPCR

In order to verify the microarray results, we selected three upregulated and three downregulated genes induced by high-dose treatments for RT-PCR assay. The results (Figure 5a) showed that all the three upregulated genes (CYP2W1, HPS3 and TNF-α) and one downregulated gene (MSH2) identified by our microarray analysis could be confirmed by semi-quantitative assessment of the amounts of mRNA/cDNA for these genes using RT-PCR. There were no obvious RT-PCR bands detected with the other two downregulated genes (data not shown). Using a more quantitative measurement of the changes in expression of these genes, we compared mRNA/cDNA levels for these genes in cells treated with high-dose OPs and in cells of the non-treated control culture by qPCR. The changes of expression of these genes in OP-treated cells, both upregulation and downregulation identified by the microarray analysis, could similarly be verified by examining the amounts of their mRNAs using qPCR (Figure 5b). However, the degree of changes in expression of these genes varied somewhat in cells treated with the three different OPs.

Figure 5

Verification of gene changes by RT-PCR and qPCR. The first-strand cDNA was synthesized using Superscript II Reverse Transcriptase from 1 μg of total RNA with an oligo-dT primer and used as a template for PCR. OP-C, OP-D and OP-M represent neuroblastoma cells treated with high-dose chlorpyrifos, dichlorvos and methamidophos, respectively. (a) The RT-PCR results are shown: Control represents neuroblastoma cells treated with 2% DMSO. Upregulated genes are marked with ↑; downregulated with ↓, and human GAPDH gene is used as an internal standard marked with —. (b) The log2 ratios for microarray assay, three high-dose OP treatments are shown. The gene changes were normalized to GAPDH. RT-PCR, reverse transcription polymerase chain reaction; qPCR, quantitative PCR; OP, organopho-sphate; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; DMSO, dimethyl sulfoxide; TNF-α, tumor necrosis factor. (A color version of this figure is available in the online journal)

Detection of cell apoptosis by TUNEL assay

Our microarray analysis revealed that many pathways associated with cell apoptosis and/or cell death had been activated in the neuroblastoma cells treated with high-dose OPs. The neuroblastoma cells in cultures treated with the OPs were examined using TUNEL assay for evidence of an increase of apoptosis. All three high-dose OP treatments evidently induced many of the neuroblastoma cells to undergo apoptosis, seen with positive green fluorescent staining in the assay, less found in the non-treated controls (Figure 6a). All three different OP treatments induced many folds increase of cell apoptosis in the cultures in comparison with those in the control cultures (Figure 6b). OP-D and OP-M appeared to be more toxic than OP-C in terms of the effectiveness of inducing apoptosis in neuroblastoma cells.

Figure 6

The TUNEL assay was carried out by DeadEnd Fluorometric TUNEL System Kit. The high-dose chlorpyrifos (OP-C), dichlorvos (OP-D) and methamidophos (OP-M) were used to treat SK-N-SH cells for 24 h. (A) The images of TUNEL positive cells were captured by a confocal microscope (× 400). (B) Quantitative result of TUNEL assay was analyzed: 100 cells were counted for each sample. The percentage of TUNEL-positive cell (green fluorescent) number to the total cell number is shown. TUNEL, terminal deoxynucleotidyl transferase-mediated dUTP nick-end labeling; OP, organophosphate; (A color version of this figure is available in the online journal)

Discussion

The main goals for our study were to elucidate the molecular mechanism of OP injury and find potential biomarkers for evaluating the degree of OP poisoning to human cells. Three representative human cell lines (alveolar epithelial cells A549, neuroblastoma cells SK-N-SH and hepatocellular cells HepG2) were initially studied and treated, respectively, with one of the three OPs (chlorpyrifos, dichlorvos or methamidophos) at high-dose (LC20) or low-dose (1/20 of LC20) for 24 h. Changes in global gene expression profiles in three different cell lines treated with each of the three OPs were analyzed by ArrayTrack. Hundreds of genes with significant changes in expression (P < 0.05 and fold changes >1.5) were identified as ‘significant genes’ for each cell line treated with the three different Ops, respectively. We found few common significant genes with similar changes in expression among the three different cell lines treated by either high- or low-dose of any OP studied (data not shown). OP-induced injury appeared to be rather tissue-specific and dose-related.

Among the three different cell lines examined in this study, human neuroblastoma cells exposed to all the three different OPs at either high or low dose clearly had the highest numbers of significant genes that were further analyzed using mAdb. We examined the significant genes in neuroblastoma cells treated with three different OPs in order to identify the common mechanisms of OP injury. In addition, we compared the significant genes and the pathways involved in the neuroblastoma cells exposed to the high and the low dose of OPs. Human neuroblastoma cells were previously used as a model test system for OP-associated injury.¹⁷ It was reported that mipafox (a well-known organophosphorus cholinesterase inhibitor) caused organophosphorus-induced delayed neuropathy not only by their effects on neurite length, but also by their effects on Ca²⁺ homeostasis in differentiated neuroblastoma cells.¹⁸

In our genome-wide analysis of the common significant genes induced by all three different OPs, we found that the high- and low-dose treatments in neuroblastoma cells were apparently affecting different genes that showed significant changes in expression. In other words, not too many common significant genes were identified in cells with high- and low-dose treatments of all three OPs (Figure 3). However, we noticed that significant genes identified in both high-dose and low-dose treatments were often involved in similar pathways. The KEGG pathway analysis revealed three out of four top significant pathways involved in the high-dose and low-dose treatments were actually the same (Table 3). Specifically, metabolic, cell cycle and MAPK signaling pathways were the three pathways found to be downregulated by both high- and low-dose treatments with all three different OPs (Tables 3, 4 and 5). MAPK signaling pathway plays an important role in several key cellular functions, including cell proliferation, differentiation and migration.¹⁹ Downregulated effects on metabolic, cell cycle and MAPK signaling pathways were consistent with the apparent biological conditions that cell growth and proliferation were markedly inhibited in the cells with both high-and low-dose treatments.

Among the downregulated significant genes (Table 2), PRIM1 (primase, DNA, polypeptide 1) and CDC25A (cell division cycle A) are two important genes known to be involved in cell growth and cell cycle processes. PRIM1 is a small unit of a heterodimer in DNA primase, which plays key roles in DNA replication and synthesis of RNA primer.²⁰ CDC25A is a member of the CDC25 family of phosphatases, and plays important roles in cell cycle regulation, particularly at G1/S and G2/M check points.^21,22 Downregulation of the PRIM1 and CDC25A would directly inhibit DNA replication, cause cell cycle arrest and result in cell death. It is important to note that changes of expression (down- or up-regulation) of the ‘significant genes’ identified by the microarray analysis could be verified by examining specifically the amounts of mRNA in neuroblastoma cells treated with various OPs using both RT-PCR and qPCR (Figures 5a and b).

In addition to the inhibitory effects found on the three pathways involved in cell metabolism, growth and proliferation, our results showed that high-dose treatments with all three different OPs (LC20) induced increased expression of genes such as ERN1, PRNP, NET1, CLN3, MT1A and MT3 in the neuroblastoma cells. ERN1 encodes a small nucleus protein which activates cell apoptotic process and induces cell cycle arrest.^23,24 NET1 gene is involved in the induction of apoptosis and regulation of cell growth.²⁵ CLN3 is a multifunctional gene which is particularly involved in the regulation of cell death and neuron apoptosis.^26,27 MT1A and MT3 play a negative role in the regulation of cell growth and neurogenesis.^28,29 Activation of these genes would facilitate the cellular apoptotic process and further inhibit cell growth, thus results in rapid cell death. The marked increase of cell apoptosis in the OP-treated cultures could also be documented by TUNEL assay (Figures 6a and b). In the present study, OP-D and OP-M appeared to be more toxic than OP-C in the effect of inducing neuroblastoma cells to undergo apoptosis.

In comparison, our studies of significant genes in the neuroblastoma cells treated with low dose of the three OPs did not identify any of those apoptotic genes activated significantly by high-dose treatments of the three OPs. In contrast, the KEGG analysis showed that low-dose treatments of all three OP induced significant decreases in the expression of the genes involved in a unique pathway with spliceosome process (Table 6 and Figure 4). Most eukaryotic genes are expressed as precursor mRNAs that are converted to mRNA by a macromolecular complex ‘spliceosome’ through splicing, an essential step of gene expression in which non-coding sequences (introns) are removed and coding sequences (exons) are ligated together.³⁰ Spliceosomes are a dynamic family of complex particles that assemble on the mRNA precursor and help fold it into a conformation that will then allow transesterification to proceed. The proteins coded by the downregulated spliceosom-associated genes, HNRNPK (heterogeneous nuclear ribonucleoprotein K³¹) and HNRNPA1 (heterogeneous nuclear ribonucleoprotein A1³²) genes, belong to the subfamily of ubiquitously expressed heterogeneous nuclear ribonucleoproteins (hnRNPs). The hnRNPs are RNA-binding proteins that complex with heterogeneous nuclear RNA (hnRNA) and pre-mRNAs in the nucleus. They participate in pre-mRNA processing, mRNA metabolism, nuclear/cytoplasm transport and modulation of splicing site selection. Our results suggest that the low-dose treatments, unlike the high-dose treatments, might not lead to rapid cell death. Instead, exposure to low-dose of OPs could significantly decrease expression in genes of the spliceosome pathway and affect the normal processing of mRNA representing a new form of OP-associated patho-genesis. In this context, a recent study in patients with mye-lodysplasia syndromes showed frequent pathway mutations of splicing machinery spliceosome.³³

Our results also showed that the expressions of some genes were commonly affected in the neuroblastoma cells exposed to all three different OPs at two different doses, respectively (Figure 3). This finding suggests that there could be a common mechanism in OP-associated injury. Identification of these common significant genes might be important for understanding the molecular mechanisms in cell injury elicited by OPs. They could also serve as potential candidates for biomarkers to evaluate the degree of OP-associated injury in neuronal cells as well as potential targets for therapeutic regimens against OP-associated toxicities.

Author contributions: All authors contributed to the design, conduct, collection and analysis of data and the preparation of the manuscript. Specially, TL, HZ and JZ conducted the experiments and collected the data; G-CH, JH, ST, BL, RKP and S-CL analyzed and wrote the manuscript.

Footnotes

Acknowledgements

This research was supported by a grant from Defense Threat Reduction Agency (DTRA: #CBS.MEDCHEM.01.10.FDA. PP008). We thank Dr Amy Yang for her help on microarray experiments and data analysis, Dr Ephrem Tekle for his help on confocal microscope.

References

Pope

CN.

Organophosphorus pesticides: do they all have the same mechanism of toxicity?

J Toxicol Environ Health B Crit Rev 1999; 2: 161–81.

Clegg

, van Gemert

Expert panel report of human studies on chlorpyrifos and/or other organophosphate exposures. J Toxicol Environ Health B Crit Rev 1999; 2: 257–79.

Eskenazi

, Bradman

, Castorina

Exposures of children to organophosphate pesticides and their potential adverse health effects. Environ. Health Perspect 1999; 107 (Suppl. 3):409–19

Vidair

CA.

Age dependence of organophosphate and carbamate neurotoxicity in the postnatal rat: extrapolation to the human. Toxicol Appl Pharmacol 2004; 196: 287–302.

Brown

, Brix

KA.

Review of health consequences from high-, intermediate- and low-level exposure to organophosphorus nerve agents. J Appl Toxicol 1998; 18: 393–408.

Bloch-Shilderman

, Rabinovitz

, Egoz

, Raveh

, Allon

, Grauer

, Gilat

, Weissman

BA.

Subchronic exposure to low-doses of the nerve agent VX: physiological, behavioral, histopathological and neurochemical studies. Toxicol Appl Pharmacol 2008; 231: 17–23.

Prendergast

, Self

, Smith

, Ghayoumi

, Mullins

, Butler

, Buccafusco

, Gearhart

, Terry

AV.

Microtubule-associated targets in chlorpyrifos oxon hippocampal neurotoxicity. Neuroscience 2007; 146: 330–9.

Niijima

, Nagao

, Nakajima

, Takatori

, Matsuda

, Iwase

, Kobayashi

Sarin-like and soman-like organophosphorous agents activate PLC gamma in rat brains. Toxicol Appl Pharmacol 1999; 156: 64–9.

Ivanović-Matić

, Poznanović

, Grigorov

, Dinić

, Mihailović

, Grdović

, Uskoković

, Martinović

, Arambasić

, Petrović

, Bogojević

The organophosphate-induced acute-phase response is characterized by synthesis of alpha 1-acid glycoprotein that exhibits an immunomodulatory effect. J Appl Toxicol 2008; 28: 63–71.

10.

Tang

, François

, Glatigny

, Agier

, Mucchielli

, Aggerbeck

, Delacroix

Expression ratio evaluation in two-colour microarray experiments is significantly improved by correcting image misalignment. Bioinformatics 2007; 23: 2686–91.

11.

Mense

, Sengupta

, Lan

, Zhou

, Bentsman

, Volsky

, Whyatt

, Perera

, Zhang

The common insecticides cyfluthrin and chlorpyrifos alter the expression of a subset of genes with diverse functions in primary human astrocytes. Toxicol Sci 2006; 93: 125–35.

12.

Slotkin

, Seidler

FJ.

Protein kinase C is a target for diverse developmental neurotoxicants: transcriptional responses to chlorpyrifos, diazinon, dieldrin and divalent nickel in PC12 cells. Brain Res 2009; 1263: 23–32.

13.

Gultkin

, Patat

, Akca

, Akdogan

, Altuntas

Melatonin can suppress the cytotoxic effects of chlorpyrifos on human HepG2 cell line. Human Exp Toxicol 2006; 25: 47–55.

14.

Slotkin

, Seidler

FJ.

Transcriptional profiles reveal similarities and differences in the effects of developmental neurotoxicants on differentiation into neurotransmitter phenotypes in PC12 cells. Brain Res Bull 2009; 78: 211–25.

15.

Greene

, Asaki

, Bian

, Bock

, Castillo

, Chandramouli

, Martell

, Meyer

, Ruppert

, Sundaram

, Tomlin

, Yang

, Powell

The NCI/CIT microArray database (mAdb) system - bioinformatics for the management and analysis of Affymetrix and spotted gene expression microarrays. AMIA Annu Symp Proc 2003; 1066

16.

Tong

, Cao

, Harris

, Sun

, Fang

, Fuscoe

, Harris

, Hong

, Xie

, Perkins

, Shi

, Casciano

ArrayTrack-supporting toxicogenomic research at the U.S. Food and Drug Administration National Center for Toxicological Research. Environ Health Perspect 2003; 111: 1819–26.

17.

Cho

, Wild

, Donnelly

, Tiffany-Castiglioni

Degradation of organophosphorus neurotoxicity in SY5Y neuroblastoma cells by organophosphorus hydrolase (OPH). J Toxicol Environ Health A 2006; 69: 1413–29.

18.

Cho

, Tiffany-Castiglioni

Neurofilament 200 as an indicator of differences between mipafox and paraoxon sensitivity in Sy5Y neuroblastoma cells. J Toxicol Environ Health A 2004; 67: 987–1000.

19.

Wancket

, Frazier

, Liu

Mitogen-activated protein kinase (MKP)-1 in immunology, physiology, and disease. Life Sci 2011; 90: 237–48.

20.

Urban

, Joubert

, Purse

, Hocek

, Kuchta

RD.

Mechanisms by which human DNA primase chooses to polymerize a nucleoside triphosphate. Biochemistry 2010; 49: 727–35.

21.

Galaktionov

, Jessus

, Beach

Raf1 interaction with Cdc25 phosphatase ties mitogenic signal transduction to cell cycle activation. Genes Dev 1995; 9: 1046–58.

22.

Vazquez-Novelle

, Mailand

, Ovejero

, Bueno

, Sacristan

MP.

Human Cdc14A phosphatase modulates the G2/M transition through Cdc25A and Cdc25B. J Biol Chem 2010; 285: 40544–53.

23.

Iwawaki

, Hosoda

, Okuda

, Kamigori

, Nomura-Furuwatari

, Kimata

, Tsuru

, Kohno

Translational control by the ER transmembrane kinase/ribonuclease IRE1 under ER stress. Nat Cell Biol 2001; 3: 158–64.

24.

Bouchecareilh

, Higa

, Fribourg

, Moenner

, Chevet

Peptides derived from the bifunctional kinase/RNase enzyme IRE1alpha modulate IRE1alpha activity and protect cells from endoplasmic reticulum stress. FASEB J 2011; 25: 3115–29.

25.

Srougi

, Burridge

The nuclear guanine nucleotide exchange factors Ect2 and Net1 regulate RhoB-mediated cell death after DNA damage. PLoS One 2011; 6: e17108

26.

Janes

, Munroe

, Mitchison

, Gardiner

, Mole

, Wallace

BA.

A model for Batten disease protein CLN3: functional implications from homology and mutations. FEBS Lett 1996; 399: 75–7.

27.

Cao

, Staropoli

, Biswas

, Espinola

, Macdonald

, Lee

, Cotman

SL.

Distinct early molecular responses to mutations causing vLINCL and JNCL presage ATP synthase subunit C accumulation in cerebellar cells. PLoS One 2011; 6: e17118

28.

Toh

, Li

, Yip

, Lo

, Guo

, Phan

, Bay

BH.

Modulation of metallothionein isoforms is associated with collagen deposition in proliferating keloid fibroblasts in vitro. Exp Dermatol 2010; 19: 987–93.

29.

Pountney

, Dickson

, Power

, Vickers

, West

, Gai

WP.

Association of metallothionein-III with oligodendroglial cytoplasmic inclusions in multiple system atrophy. Neurotox Res 2011; 19: 115–22.

30.

Will

, Luhrmann

Spliceosome structure and function. Cold Spring Harb Perspect Biol 2011; 3: a003707

31.

Mikula

, Bomsztyk

Direct recruitment of ERK cascade components to inducible genes is regulated by heterogeneous nuclear ribonucleoprotein (hnRNP) K. J Biol Chem 2011; 286: 9763–75.

32.

Michlewski

, Guil

, Caceres

JF.

Stimulation of pri-miR-18a processing by hnRNP A1. Adv Exp Med Biol 2010; 700: 28–35.

33.

Yoshida

, Sabada

, Shiraishi

, Nowak

, Nagata

, Yamamoto

, Sato

, Sato-Otsubo

, Kon

, Nagasaki

, Chalkidis

, Shiosaka

, Kawahata

, Yamaguchi

, Otsu

, Obara

, Sakata-Yanagimoto

, Ishiyama

, Mori

, Nolte

, Hofmann

, Miyawaki

, Sugano

, Haferlach

, Koeffler

, Shih

, Haferlach

, Chiba

, Nakauchi

, Miyano

, Ogawa

Frequent pathway mutations of splicing machinery in myelodysplasis. Nature 2011; 478: 64–9.