Profiling transcriptomes of human SH-SY5Y neuroblastoma cells exposed to maleic acid

Reagents

All reagents were purchased from Sigma (St. Louis, MO, USA) unless otherwise stated.

Cell lines and culture conditions

The human neuroblastoma cell line, SH-SY5Y, was purchased from ATCC (Manassas, VA, USA). Cells were cultured in a 1:1 mixture of ATCC-formulated Eagle’s minimum essential medium and F12 medium supplemented with 2 mM l-glutamine, 1% of penicillin/streptomycin, and 10% heat-inactivated fetal bovine serum (GIBCO, Gaithersburg, MD, USA) at 37 °C in 5% CO₂. Cells were seeded at an initial density of 10⁴ cells/cm² and cultured overnight prior to the treatment with various concentrations of maleic acid (Sigma-Aldrich, St. Louis, MO).

Cell viability

Cell viability was evaluated by 3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-((4-sulfophenyl)-2H-tetrazolium) chloride (MTS) assay (The CellTiter 96® AQ_ueous One Solution Cell Proliferation Assay kit) according to the manufacturer’s instruction. Briefly, cells (4 × 10⁴) were seeded in a 96-well plate overnight and then exposed to maleic acid for 24 and 48 h. MTS solution was added to each well for 4 h at 37 °C and the absorbance was read at wavelength 490 nm by a microplate reader (Thermo Varioskan Flash, Thermo, Winooski, VT, USA).

Measurement of intracellular calcium

Fluo-4 AM, a membrane-permeable Ca²⁺-sensitive fluorescent dye, was used to measure the free cytosolic calcium. Cells were seeded at 2 × 10⁴ cells/well in 96-well clear-bottomed black microplates (Corning Costar, Tewksbury, MA, USA) for overnight. Cells were treated with serial concentrations of maleic acid and the supernatants of each well were removed. Next, the cells were washed with PBS twice and incubated with 5 μM of Fluo-4 AM (ThemoFisher) at 37 °C for 15 min. The change in fluorescence at excitation/emission of 494/506 nm was analyzed by a microplate reader (Thermo Varioskan Flash, Thermo, USA).

Measurement of intracellular thiols

The depletion of intracellular thiols in maleic acid-treated cells was determined by CMF-DA staining. SH-SY5Y cells were seeded at 2 × 10⁴ cells/well in 96-well clear-bottomed black microplates overnight. Then, the culture medium was replaced with new medium containing maleic acid (0.1–100 μM) for 1, 8, and 24 h. After treatment, cells were washed with PBS and then stained with 25 μM of CMF-DA for 30 min at 37 °C. Then, the cells were washed twice with PBS. The change in fluorescence at excitation/emission of 488/525 nm was analyzed by a microplate reader (Thermo Varioskan Flash, Thermo, USA).

Total RNA isolation

Total RNA was extracted from SH-SY5Y cells using the TRIzol (Qiagen, Valencia, CA, USA). The RNA samples were then treated with an RQ1 RNase-free DNase kit (Promega, Southampton, UK) according to the manufacturer’s instructions to remove contaminated DNA. The RNA concentration and purity were checked by OD260/OD280 (>1.8) and OD260/OD230 (>1.6), and the yield and quality were assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Human oligonucleotide DNA microarray

The Human Whole Genome OneArray® v6 (Phalanx Biotech Group, Hsinchu, Taiwan) contains 32,679 DNA oligonucleotide probes, and each probe is a 60-mer designed in the sense direction. Among the probes, 31,741 probes correspond to the annotated genes in RefSeq v51 and Ensembl v65 database. Besides, 938 control probes are also included. The detailed descriptions of the gene array are available from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL19137.

Microarray analysis

Fluorescent antisense RNA (aRNA) targets were prepared from 1 μg total RNA samples using OneArray® Amino Allyl aRNA Amplification Kit (Phalanx Biotech Group, Taiwan) and Cy5 dyes (Amersham Pharmacia, Piscataway, NJ, USA). Fluorescent targets were hybridized to the Human Whole Genome OneArray® with phalanx hybridization buffer using phalanx hybridization system. After hybridization at 50 °C for 16 h, non-specific binding targets were washed out three times. First, the slides were washed twice at 42 °C for 5 min followed by another wash at 25 °C for 5 min. The slides were then rinsed for 20 times and dried by centrifugation. Finally, the slides were scanned via an Agilent G2505C scanner (Agilent Technologies, Santa Clara, CA, USA). The Cy5 fluorescent intensities of each spot were analyzed by GenePix 4.1 software (Molecular Devices, Sunnyvale, CA, USA). The signal intensity of each spot was processed by Rosetta Resolver System® (Rosetta Biosoftware, Seattle, WA, USA). The error model of Rosetta Resolver System® could provide a judgment of the probe signal’s credibility by removing both systematic and random errors within the data. Unreliable spots with a flag less than zero were filtered out. Spots that passed the criteria were normalized by 50% median-scaling normalization method. The technical repeat data was tested by Pearson correlation coefficients to check the reproducibility (r > 0.975). Normalized spot intensities were transformed into gene expression log₂ ratios between the control and treatment groups. The data discussed in this study have been deposited in NCBI’s Gene Expression Omnibus (Edgar, Domrachev & Lash, 2002) and are accessible through GEO series accession number GSE86510.

For the identification of differentially expressed genes, fold changes of probes are first calculated based on the ratio error model from Rosetta Resolver System® (Weng et al., 2006). Probes with a fold change larger than two are considered as differentially expressed genes. Compared to the conventional t-test, the error model helps to reduce false positives in comparison of gene expressions in two treatments with multiple replicates. In addition to the identification of differentially expressed genes based on the Rosetta error model, the average values of technical replicates were also calculated for creating a heatmap for visualization as shown in Fig. 1.

Enrichment analysis

Enrichment-analysis tools were implemented and integrated into ChemDIS (Tung, 2015) in order to analyze functions, pathways, and diseases from a given set of differentially expressed genes. For enriched functions, clusterProfiler (Yu et al., 2012) providing comprehensive enrichment analysis functionalities was applied to analyze enriched gene ontology (GO) terms for MF, biological process (BP), and cellular component (CC). The enriched KEGG (Kanehisa et al., 2014) and Reactome (Croft et al., 2014) pathways were analyzed using clusterProfiler and ReactomePA (Yu & He, 2016), respectively. For inferring diseases affected by a given chemical, enriched DO (Kibbe et al., 2015), and DOLite (Du et al., 2009) terms will be analyzed using DOSE package (Yu et al., 2015). DOLite is a simplified vocabulary list of DO, a standardized ontology connecting human proteins to diseases. All the enrichment analyses are based on hypergeometric tests with the Benjamini–Hochberg approach (Du et al., 2009) for multiple-testing corrections. Enriched terms with a corrected p-value <0.05 will be identified. The enrichment analysis function of ChemDIS is freely available at http://cwtung.kmu.edutw/chemdis and the analysis in this study is reproducible using the service.

Real-time reverse transcription-polymerase chain reactions

The concentration of RNA was determined by measuring its absorbance at 260 nm with a microplate reader (Thermo Varioskan Flash, Thermo, USA). Two micrograms of total RNA collected from each sample were reverse-transcribed through RevertAid RT Kit (Thermo) into cDNA products using oligo-dT as the primer. Real-time polymerase chain reaction (PCR) was performed in a 96-well optical tray by an ABI PRISM® 7900HT Sequence Detection System (Applied Biosystems, Cheshire, UK). During the real-time reverse transcription-polymerase chain reaction (RT-PCR) process, we used Luminaris Color HiGreen High ROX qPCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA) which offered a highly specific and sensitive method to quantify mRNA expressions. The HPRT gene was used as an endogenous control to normalize the expressions of target genes. The primers are: 5′-GTTCTGGACACCAACAAGGAC3′ and 5′-TCGTGGCAGTAGAGACA-GAGG-3′ for S100A3, 5′-AACTACTCTGGCAGGATCTGGA-3′ and 5′-AGGGGG-TTAGGGATGAGGCT-3′ for GGT5, 5′-GCTAAAGGTTTTGGAGCCCAC-3′ and 5′-AGAGTGGAGAGACATAGAGATGGT-3′ for GRIA1, 5′-GTTCTCCTTCATCCTATTTATGTGG-3′ and 5′-CCCTGTCCAAATCCTAACCCC-3′ for GNRHR, and 5′-TCAGTCAACGGGGGACATAAA-3′ and 5′-GGGGCTGTACTGCTTAACCAG-3′ for HPRT.

Statistical analysis

The data of cell viability, cellular levels of calcium and thiols, and RNA expression were presented as mean values and standard errors. One-way analysis of variance was adopted for multiple comparisons, and Dunnett’s two-tailed t-test was used to compare the results of the treatment groups with those of the control group using JMP statistical software (Cary, NC, USA). p < 0.05 was defined as statistically significant.

Transcriptome profiling and differentially expressed genes

To understand more about the effects of maleic acid on neuronal cells, the gene expression profiles of SH-SY5Y cells treated with three doses (10, 50, and 100 μM) of maleic acid were examined. The assessment of variability, which comes from biological and technical replicates, is an important concern for microarray experiments. However, a high intra-platform consistency of microarray experiments for both technical and biological replicates has been observed with correlation coefficients higher than 0.9 (Tan et al., 2003). Compound-induced changes have also been shown to outweigh the influence of the technical and biological variability in a microarray experiment (Bryant et al., 2011).

To minimize the selection bias caused by biological variability, this study adopted a scheme of three treatments with different doses in the experimental group versus one control, and only differentially expressed genes with monotonically increased or decreased expression values were considered for the following analysis. In addition, two technical replicates of each sample were utilized to confirm the reproducibility. The monotonically increased/decreased genes were identified by applying simple filtering rules. Given a candidate gene with four expression values, the following rule was applied to define a monotonically increased gene. IF (H > M) AND (M > L) AND (L > N) THEN the gene is a monotonically increased genes, where H, M, L, and N represent its expression values for the treatments of maleic acid in different doses (100, 50, 10, and 0 μM), respectively.

A total of 389 differentially expressed probes (191 upregulated and 198 downregulated) were identified with more than twofold changes over control after applying the highest dose (100 μM) of maleic acid (Table S1). To remove probes without corresponding gene annotations, 46 and 9 probes with keywords of “NA” and “uncharacterized,” respectively, in the corresponding description field were excluded from subsequent analysis. Finally, 243 differentially expressed genes were identified (113 upregulated and 130 downregulated). Detail information for the 243 differentially expressed genes is available in Table S1. Among the 243 differentially expressed genes identified by the Rosetta error model (Weng et al., 2006), 48 genes were also found to have a clear dose-dependent trend in response to maleic acid treatment based on the average expression values from technical replicates. Figure 1 shows the heatmap of the 48 genes, including 42 upregulated and six downregulated genes.

Enrichment analysis

The analysis of enriched GO terms was conducted for the 243 differential expressed genes to identify affected functions associated with maleic acid. As shown in Table 1, there were 12 significantly enriched GO terms (corrected p-value <0.05), including eight MFs, two CCs, and two BPs. The enriched MF GO terms were visually analyzed as shown in Fig. 2. Two main branches with hierarchy levels higher than binding (GO:005488) were identified. Metal ion binding (GO:0046872), cation binding (GO:0043169), and ion binding (GO:0043167) are located at the same branch, whereas DNA binding (GO:003677), nucleic acid binding (GO:0003676), heterocyclic compound binding (GO:1901363), and organic cyclic compound binding (GO:0097159) are located at the other branch. Both branches of GO terms of metal ion binding and DNA binding were consistent with our previous analysis from ChemDIS system (Lin, Wang & Tung, 2014; Tung, 2015). The calcium binding (GO:0005509), a direct descendant of the enriched term of metal ion binding (GO:0046872), was inferred to be associated with maleic acid. The enriched GO term of DNA binding (GO:0003677) is related to nucleotide binding (GO:0000166) that was inferred to be associated with maleic acid from ChemDIS (Lin, Wang & Tung, 2014; Tung, 2015). The enriched pathways are hematopoietic cell lineage (KEGG ID:hsa04640) and generic transcription pathway (Reactome ID: 212436), and there are no enriched DO and DOLite terms from the analysis. For validating the enrichment analysis in this study, the same gene set was applied to analyze enriched terms using the comparative toxicogenomics (CTD) database. As expected, except for the analysis of enriched disease terms, ChemDIS and CTD generated similar enriched GO and pathway terms as shown in Table S2. ChemDIS and CTD utilizing different disease ontologies could generate different results of enriched diseases (Tung, 2015).

Attenuation of cellular calcium and thiol levels by maleic acid

According to our previous transcriptome and enrichment analysis, calcium binding was inferred to be affected by maleic acid. Among the differentially expressed genes, eight genes are associated with the GO term of calcium binding (GO:0005509) including S100A3, GNPTAB, CDH23, PCDHAC1, PADI1, CRACR2B, MAN1A2, and CDH8. In addition, GNRHR and GRIA1 among the differentially expressed genes are involved in the regulation of cellular calcium levels. Altogether, cellular calcium is potentially affected by maleic acid. To validate this finding, the kinetic change of cellular calcium and thiol levels by maleic acid was analyzed by specific fluorescence Fluo-4 and CMF-DA probes, respectively. To avoid measuring the direct cytotoxic effects of maleic acid on SH-SY5Y cells, we first evaluated the cell viability of maleic acid-treated SH-SY5Y cells by an MTS assay. Figure 3A showed that there was no significant change in cell viability after 48 h of maleic acid treatment. By contrast, maleic acid (0.1–100 μM) significantly decreased the cellular calcium level by 12–41%, determined by the change of fluorescence intensity of Fluo-4, in a concentration- and time-dependent manner (Fig. 3B). Maleic acid (1–100 μM) decreased the cellular thiol levels by 16–31% after 24 h treatment (Fig. 3C). The exact p-values of Fig. 3 were shown in Supplemental Information 2. These results indicated that maleic acid resulted in a loss of cellular calcium level earlier than glutathione depletion. These observations could give the evidence that maleic acid at the concentrations higher than 0.1 μM interferes cellular calcium homeostasis in human neuronal cells. The homeostasis of intracellular Ca²⁺ plays important roles in the regulation of physiological responses or controlling cell death (Demuro et al., 2005; Bading, 2013; Brini et al., 2014). Depletion of intracellular Ca²⁺ has been reported to induce SH-SY5Y cell death by activating caspase activation (McGinnis, Wang & Gnegy, 1999). Our data showed that maleic acid significantly attenuated intracellular Ca²⁺ levels without induction of cell death, indicating that other pathways might be involved and need to be further studied.

qPCR validation

According to our transcriptome analysis, there were 243 genes whose expression were significantly affected by maleic acid in a dose-dependent manner. To confirm the association between the gene regulation and the attenuated cellular calcium and thiol levels by maleic acid, we evaluated the expression of GGT5, GNRHR, GRIA1, and S100A3 by real-time qPCR using three independent biological replicates to further prove the reliability of the microarray data.

Maleic acid upregulated the expression of GGT5 in a dose-dependent manner (Fig. 4A). Gamma-glutamyl transferase (GGT) family genes encode key enzymes widely distributed in various tissues in response to different stimuli (Lieberman et al., 1995). They might play important roles in glutathione salvage, detoxification of xenobiotics, and metabolism of endogenous mediators such as prostaglandins and leukotrienes (Lieberman et al., 1995). In neurons, GGT enzymes liberate glutamate from glutathione. Glutathione has been suggested as a reservoir of neuronal glutamate (Koga et al., 2011), which is the main excitatory transmitter in the central nervous system to participate in many important physiological processes including long-term potentiation and developmental plasticity (McEntee & Crook, 1993). Deficiency of glutathione has been reported in several neuropsychiatric disorders (Gu, Chauhan & Chauhan, 2015). Dysregulation of GGT may interfere with the glutathione synthesis and cellular glutamate levels.

Maleic acid also upregulated GRIA1, GNRHR, and S100A3 in a dose-dependent manner (Figs. 4B–4D). The exact p-values of Fig. 4 were shown in Supplemental Information 2. GRIA1, GNRHR, and S100A3 are associated with calcium influx or binding. GRIA1 encodes glutamate receptor 1, which is a subunit of calcium-permeable α-amino-3-hydroxy-5-methyl-4-isoxazolepropionate (AMPA) glutamate receptors. Glutamate receptors are activated in a variety of normal neurophysiologic processes and are responsible for the glutamate-mediated postsynaptic excitation of neuronal cells. GNRHR gene encodes gonadotropin-releasing hormone receptor, which belongs to seven-transmembrane G-protein coupled receptor (Millar, 2005). Activation of GNRHR will further activate a phosphatidylinositol-calcium second messenger system and further regulate the release of follicle-stimulating hormone and luteinizing hormone (Romanelli et al., 2004). S100A3 encodes S100 calcium-binding protein A3 which contains two EF-hand calcium-binding motifs and plays roles in cell cycle progression and differentiation (Schäfer & Heizmann, 1996). The increased expression of S100A3 has been associated with methamphetamine-mediated neuronal toxicity (Cadet et al., 2011). As GNRHR, GRIA1, and S100A3 were all significantly augmented by maleic acid, the finding suggested that the dysregulation of these genes may participate in maleic acid-mediated toxicity.

Collectively, the present data demonstrated that maleic acid interfered with the expression level of genes associated with calcium influx or binding and further gave the mechanical evidence on how maleic acid affects human neuronal cells. However, the alteration of functions of the neuronal cells by these proteins in response to maleic acid treatment was not evaluated in this study and warrant further study. Nonetheless, the study has set an example of utilizing ChemDIS system to prioritize the potential effects of maleic acid on the human health and to generate testable experimental approaches which could validate the inferred association. This approach might also be useful for understanding other poorly characterized chemicals.