Identification of molecular signatures predicting the carcinogenicity of polycyclic aromatic hydrocarbons (PAHs)

Toxicology Letters 212 (2012) 18–28

Contents lists available at SciVerse ScienceDirect

Toxicology Letters journal homepage: www.elsevier.com/locate/toxlet

Identification of molecular signatures predicting the carcinogenicity of polycyclic aromatic hydrocarbons (PAHs) Mi-Kyung Song a,b , Mee Song a , Han-Seam Choi a , Youn-Jung Kim c , Yong-Keun Park b , Jae-Chun Ryu a,∗ a b c

Cellular and Molecular Toxicology Laboratory, Korea Institute of Science & Technology, P.O. Box 131, Cheongryang, Seoul 130-650, Republic of Korea School of Life Sciences and Biotechnology, Korea University, Anam-Dong, Seoungbuk-Gu, Seoul 136-791, Republic of Korea Department of Food and Nutrition, Eulji University, 212 Yangji-Dong, Sujueong-Gu, Seongnam-Si, Gyeonggi-Do 461-713, Republic of Korea

a r t i c l e

i n f o

Article history: Received 15 September 2011 Received in revised form 17 April 2012 Accepted 18 April 2012 Available online 1 May 2012 Keywords: Polycyclic aromatic hydrocarbons (PAHs) Carcinogenicity Microarray Molecular signature

a b s t r a c t Assessing the potential carcinogenicity of human toxins represents an ongoing challenge. Chronic rodent bioassays predict human cancer risk with limited reliability, and are expensive and time-consuming. To identify alternative prediction methods, we evaluated a transcriptomics-based human in vitro model to classify carcinogens by their modes of action. The aim of this study was to determine the transcriptomic response and identify specific molecular signatures of polycyclic aromatic hydrocarbons (PAHs), which can be used as predictors of carcinogenicity of environmental toxins in human in vitro systems. We found that characteristic molecular signatures facilitate identification and prediction of carcinogens. To evaluate the change in gene expression levels, human hepatocellular carcinoma (HepG2) cells were exposed to nine different PAHs (benzo[a]pyrene, dibenzo[a,h]anthracene, 3-methylcholanthrene, naphthalene, chrysene, phenanthrene, benzo[a]anthracene, benzo[k]fluoranthene, and indeno[1,2,3-c,d]pyrene) for 48 h. Unsupervised gene expression analysis resulted in a characteristic molecular signature for each toxin, and a supervised analysis identified 31 outlier genes as distinct molecular signatures distinguishing carcinogens from noncarcinogens. Further analysis and multi-classification revealed 430 genes as surrogate markers for predicting carcinogenic potencies of each PAH with 100% accuracy. Our results suggest that these expression signatures can be used as predictable and discernible surrogate markers for detecting and predicting PAH exposure, and their carcinogenic potential. Furthermore, the use of these markers can be more widely applied in combination with traditional techniques for assessing and predicting toxic exposure to PAHs. © 2012 Elsevier Ireland Ltd. All rights reserved.

1. Introduction The human population is experiencing increased exposure to existing and novel chemicals and environmental toxins. Among the various environmental toxins, polycyclic aromatic hydrocarbons (PAHs) are an important class of environmental pollutant, and many PAHs are known or suspected carcinogens. PAHs form in the environment through incomplete combustion of organic fuels and cover a wide range of compounds. Exposure to some of these compounds is known to cause cancer in mice and is suspected to be carcinogenic in humans (Wang and Busby, 1993; Nesnow et al., 1998; Boffetta et al., 1997). Although PAHs have similar structural properties, their carcinogenic potency is variable. Despite the fact that PAHs exist at low concentrations in the environment, exposure to a combination of environmental PAHs may be carcinogenic. Based on their mechanism of action, PAHs are classified as known

∗ Corresponding author. Tel.: +82 2 958 5070; fax: +82 2 958 5059. E-mail addresses: [email protected], [email protected] (J.-C. Ryu). 0378-4274/$ – see front matter © 2012 Elsevier Ireland Ltd. All rights reserved. http://dx.doi.org/10.1016/j.toxlet.2012.04.013

or probable human carcinogens (groups 1 and 2A, respectively), possibly carcinogenic in humans (group 2B), and not a classifiable carcinogen (group 3) by the International Agency for Research on Cancer (IARC). Details are available on the IARC Web site (http://monographs.iarc.fr/ENG/Classification/index.php). Prediction of carcinogenicity is essential to assess chemical toxicity, and the effects of environmental toxic substances are still a major issue for predicting potential human health risks. The carcinogenic potential of environmental toxins is currently evaluated using chronic rodent bioassays. This strategy, however, generates a high false-positive rate and is therefore not reliable, while also being expensive and time-consuming since it requires the use of many animals and large quantities of the test compound during a typical 2-year study period. A wide range of alternative approaches are being explored and developed which are aimed at providing mechanistic information and at predicting the carcinogenic potential of a new chemical compound without conducting the 2-year bioassay in two rodent species. Cell transformation assays with cultured cells, which are thought to mimic some stages of in vivo multistep carcinogenesis, may be used to predict a carcinogenic

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

potential in general. In vivo transgenic mouse models, which should improve the sensitivity of the bioassay, are being employed in a 6 months design (Spalding et al., 2000). In case of several carcinogens, diverse parameters deduced from the suspected mechanisms have been estimated, including the potential for immunosuppression (Cohen, 2004), prechronic liver lesions to induce chronic hyperplasia (Tennant et al., 1991), induction of cell proliferation, or the measurement of DNA synthesis during acute studies (Cunningham, 1996). Results from examinations of these parameters were encouraging but sometimes ambiguous. Thus, although many of these methods support the estimation of a compound’s carcinogenic potential, none has been sufficiently evaluated to replace the 2-year bioassay. In addition an increased mechanistic understanding is required to be able to evaluate the human relevance of rodent bioassay outcomes, for which gene expression analysis may be helpful (Cohen, 2004; Jacobs and Jacobson-Kram, 2004). Novel techniques, such as multiple gene expression analysis using microarrays, allows for a more comprehensive study on the effects and mechanisms of environmental pollutants. Microarray tools have been widely used for comprehensive gene expression analysis, as well as mutations and single nucleotide polymorphism detection. In particular, large-scale microarray analysis of gene expression can detect significant changes in thousands of genes (Kushida et al., 2006; Te Pas et al., 2007). The advance in highthroughput technology provides an opportunity to identify genes associated with cellular toxicity resulting from toxin exposure (Ellinger-Ziegelbauer et al., 2008). These techniques are efficient since many genes can be studied in one experiment, and it can characterize pathways and mechanisms, as opposed to extrapolating information based on the assessment of a few genes. In addition, multiple gene expression levels can be applied as ‘fingerprint’ biomarkers. Since the liver is the main organ that metabolizes compounds (including procarcinogens), and also represents a major target organ for chemical carcinogens in vivo, in vitro models based on hepatocyte systems are frequently used to predict toxicity. Among these, primary hepatocytes, precision cut liver slices, and hepatic cell lines are well established in vitro models (Dambach et al., 2005; Butterworth et al., 1989; Groneberg et al., 2002). Combining these in vitro models with powerful genomic-based methods may provide high-throughput screening methods to predict toxin risks in humans. The application of DNA microarray technology enables examination of differential gene expression for many genes simultaneously. These DNA microarray technologies have the potential to both enhance our understanding of the mechanism underlying a compound’s carcinogenic effects and identify a characteristic set of genes from a database of reference profiles, which may allow for the prediction of an unknown compound’s mode of action (Choi et al., 2011; Song et al., 2011a,b). Because PAHs have different features according to their carcinogenicity, each PAH may induce a distinct gene expression profile, which can then be used for mechanism-based classification of unknown compounds as carcinogenic or noncarcinogenic (Hu et al., 2004; Van Delft et al., 2004). The aim of this study was to investigate changes in gene expression patterns and identify a characteristic signature to predict carcinogenicity of environmental toxins. Therefore, we exposed human hepatoma cells (HepG2) to nine PAHs selected based on their ubiquitous distribution throughout the environment and their range in carcinogenic potency. The compounds evaluated were BP (group 1), dibenzo[a,h]anthracene (DBA; group 2A), 3-methylcholanthrene (3MC; not classified), naphthalene (NP; group 2B), chrysene (CHR; group 2B), phenanthrene (PH; group 3), benzo[a]anthracene (BA; group 2B), benzo[k]fluoranthene (BF; group 2B), and indeno[1,2,3c,d]pyrene (IND; group 2B). We performed global gene expression analysis of PAH-exposed HepG2 cells and identified the consequent

19

changes in gene expression profiles, predicting different carcinogenic potencies of PAHs by genomic expression analysis. Thus, identifying the global genomic response through transcriptome changes in in vitro system may provide an important tool for identifying indicators of carcinogenicity using a noninvasive method. 2. Materials and methods 2.1. Chemicals and reagents BP, NT, 3MC, dimethylsulfoxide (DMSO) and 3-(4,5-dimethylthaizol-2-yl)-2,5diphenyltetrazolium bromide (MTT) were purchased from Sigma–Aldrich Co. (USA) and BA and BF were purchased from Acros organics (USA). DBA and PH were the products of TCI-EP (Japan) and CHR and IND were purchased from Riedel-dehaën (Germany) and WAKO (Japan), respectively. Culture media and buffer solutions were purchased from GIBCOTM (USA): Dulbecco’s Modified Eagle Medium (DMEM), Dulbecco’s Phosphate Buffered Saline (PBS), Fetal Bovine Serum (FBS) and antibiotics (penicillin and streptomycin). All other chemicals used were of analytical grade or the highest grade available. 2.2. Cell culture The human hepatocellular carcinoma cell line, HepG2 cells, was purchased from Korean Cell Line Bank (KCLB, Korea) and was maintained under a humidified atmosphere of 5% CO2 and 95% air at 37 ◦ C. The culture medium was 90% Dulbecco’s modified Eagle’s medium (GIBCO BRL, USA) supplemented with 10% fetal bovine serum (GIBCO BRL), sodium bicarbonate (Sigma, USA), sodium pyruvate (GIBCO BRL) and penicillin and streptomycin (GIBCO BRL). The medium was refreshed every 2–3 days. 2.3. Determination of cell viability To determine the cytotoxicity and effects on cell growth, the MTT cell proliferation assay was performed using the modifications described by Mosmann (Mosmann, 1983). HepG2 cells were seeded at a seeding density of 80 × 104 cells/ml on a well in 500 ␮l of media in 24-well plate. And cells were exposed to various concentrations of PAHs in culture medium at 37 ◦ C for 48 h exposure time. After exposure, the cells were incubated for 3 h with 4 mg/ml MTT in PBS. To quench the reaction, the medium was removed and DMSO was added and transferred to 96-well plate. The optimal density (O.D) of the purple formazan product was measured at a wavelength of 540 nm. The 20% inhibitory concentration (IC20 ) of cell proliferation was defined as the concentration that causes a 20% reduction in the cell viability vs. the solvent DMSO treated control. The IC20 value was calculated according to the linear equation of each curve (linear equations were derived from dose-dependent result of MTT assay of each compounds). The MTT assay was performed in triplicate for each sample. 2.4. RNA extraction Total RNA was extracted from the PAHs-treated HepG2 cells using Trizol (Invitrogen, USA) and purified using an RNeasy mini kit (Qiagen, USA), according to the manufacturer’s instructions. Genomic DNA was removed using an RNase-free DNase set (Qiagen) during RNA purification. The amount of each RNA sample was measured by the NanoDrop ND 1000 spectrophotometer (NanoDrop Technologies Inc., USA) and its quality was checked by ExperionTM (Bio-Rad Laboratories, USA). Triplicate analysis was performed for each chemical, simultaneously. 2.5. Oligonucleotide microarray hybridization Gene expression analysis was conducted using a 44 K whole human genome microarray (Agilent Technologies, USA). Labeling and hybridization were performed using a FairPlay microarray labeling kit (Stratagene, USA), followed by the coupling of the Cyanine (Cy3) dye for the control (DMSO) or Cyanine (Cy5) dye for the treated samples. Each total RNA sample (15 ␮g) was labeled with Cy3 or Cy5 conjugated dCTP (Amersharm, Piscataway, NJ) via a reverse transcription reaction. The labeled cDNA mixture was then concentrated via ethanol precipitation. The concentrated Cy3 and Cy5-labeled cDNAs were resuspended in 10 ␮l of hybridization solution. Hybridization was performed in a hybridization oven at 62 ◦ C for 12 h. After a series of washes (2× SSC/0.1% SDS for 2 min at 58 ◦ C, 1× SSC for 3 min at room temperature, 0.2× SSC for 2 min at room temperature), the slides were dried by centrifugation at 800 rpm for 3 min at room temperature. The hybridized slides were scanned using a GenePix 4000B microarray scanner (Axon Instruments, USA) and the images were analyzed using GenePix 4.1 software (Axon Instruments) to obtain gene expression ratios. 2.6. Microarray quality control and data processing Substandard spots, identified by visual examination of each slide, were flagged and excluded from further analysis. Spots harboring dust artifacts or spatial defects

20

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

were manually flagged and excluded. In an effort to filter out unreliable data, spots with signal-to-noise (signal-background SD) ratios below 10 were not included in the data analysis. Data were normalized by Global, LOWESS, print-tip, and scaled normalization for data reliability. Residual dye bias after normalization was removed by taking the average of forward and reverse fluorescent log2-ratios for each of the samples where both forward and reverse labeling were done. The spot-specific dye bias was estimated from these samples by taking half of the average difference between forward and reverse log2-ratios. For some samples, reverse labeling was not successful, so the estimate of the dye bias was subtracted from the log2-ratio. Scatter-plot analysis was performed using Microsoft Excel 2000 (Microsoft Corp., USA). A significance analysis of microarray (SAM) was performed for genes with significant changes in expression (Tusher et al., 2001). Computing a p-value for each gene assessed the statistical significance of the differential expression of genes. To determine the p-value, a permutation procedure was used and for each permutation, two-sample t statistics were computed for each gene. Genes were considered differentially expressed when logarithmic gene expression ratios in three independent hybridizations were more than 1.5 or less than 0.66, i.e., 1.5-fold difference in expression level, and when the q-value was <5. These data were loaded into GeneSpring GX V11.0.1 software (Agilent Technologies, CA) for all further analyses. To avoid an influence of technical variance on the following gene ranking and classification exercises due to experiment of all data used over a long time period, all expression data from the treated samples were normalized to their matched vehicle control expression data. Also, a distinction between so-called “Present” and “Absent” genes was not made, thus all genes were included in the normalization and the following analyses. For data visualization, hierarchical clusters were constructed with the statistically significant (p < 0.05) genes. Principle component analysis (PCA) provides a means to reduce high-dimensional gene expression data into few principal components. Functional categorization and pathway construction were performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://david.abcc.ncifcrf.gov/kegg.jsp). 2.7. Statistical analysis Unsupervised hierarchical clustering, principal component analysis (PCA) were performed using GeneSpring GX V11.0.1 software. Data reduction to define the genes most useful in class distinction was performed using a variety of metrics. Genes selected by the various metrics were used in supervised learning algorithms to build classifiers that could identify the carcinogen subgroups. We utilized an algorithm based on the Euclidean distance and average linkage to separate genes exhibiting similar patterns. To identify significantly differentially expressed genes between treatments an ANOVA with a Benjamini & Hochberg FDR multiple testing correction was used. For class prediction, the support vector machine implemented in the Class Prediction tool of GeneSpring GX was used to classify the data. 2.8. Quantitative real time RT-PCR Messenger RNA expression levels for the genes of interest were analyzed via quantitative real time reverse transcription polymerase chain reaction (RT-PCR) using a Bio-Rad iCycler system (Bio-Rad). Total RNA was reverse-transcribed into cDNA using an Omniscript RT kit (Qiagen). Primer specificity was tested by running a regular PCR for 40 cycles (95 ◦ C for 20 s and 60 ◦ C for 1 min), followed by agarose gel electrophoresis. Real time RT-PCR was performed using a SYBR supermix kit (Bio-Rad). Samples were subjected to 45 cycles of 95 ◦ C for 20 s and 60 ◦ C for 1 min. PCR efficiency was determined by running serial dilutions of template cDNA and melting curve data were collected to assure PCR specificity. Each cDNA sample was analyzed in triplicate and the corresponding no-RT mRNA sample was included as a negative control. A glyceraldehyde-3-phosphate dehydrogenase (GAPDH) primer was included in every plate as an internal loading control. The mRNA level of each sample for each gene was normalized against that of GAPDH mRNA. The relative mRNA level was determined as 2[(Ct/GAPDH − Ct/gene of interest)] . All data are presented as the mean ± standard deviation (SD) of three separate experiments. The sequences of each primers used for the quantitative real time RT-PCR are listed in Table 1.

3. Results 3.1. Cytotoxicity of PAH in HepG2 cells The relative survival of HepG2 cells following exposure to a range of concentrations of nine PAHs was determined using the MTT assay. The survival percentage, relative to the solvent control (DMSO), was determined as a percentage of the control optical density measured after treatment. Based on the results of the MTT assay, the 20% cell-viability inhibitory concentrations (IC20 ) of each compound were calculated (A linear regression fit to the data was not shown). Dose-dependent cell viability curves were

Table 1 Primer sequences for the quantitative Real time RT-PCR used in this study. Gene

Accession No.

IER3

NM 003897

TIMP1

NM 003254

LGALS1

NM 002305

RHOC

NM 175744

NDC80

NM 006101

PAH

NM 000277

GAPDH

NM 002046

Primer sequence (5 →3 ) F R F R F R F R F R F R F R

GAACTGCGGCAAAGTAGGAGAAGA TTGAGATGCTGGAGGATGCAGTCA AGCGAGGAGTTTCTCATTGCTGGA AGCCAACAGTGTAGGTCTTGGTGA TCGCCAGCAACCTGAATCTCAAAC TTGTTGCTGTCTTTGCCCAGGTTC TAAGAAGGACCTGAGGCAAGACGA ACTCAAGGTAGCCAAAGGCACTGA ACATTGTGGCAGCCTTAGTTTGGC TCAAAGCTGTCGGCACCACTCATA CTGTTGGGACATGTGCCCTTGTTT TTGTTTGCAGAGCCCAAACTCCAC TGCACCACCAACTGCTTAGC GGCATGGACTGTGGTCATGAG

obtained after 48 h exposure of the HepG2 cells to the three PAHs, as shown in Fig. 1. The IC20 values were 8.35 ␮M, 3.79 ␮M, 2.14 ␮M, 1772.23 ␮M, 34.97 ␮M, 500 ␮M, 36.09 ␮M, 22.56 ␮M, and 227.66 ␮M for BP, DBA, 3MC, NP, CHR, PH, BA, BF, and IND, respectively. 3.2. Expression profiling of PAH provided a molecular signature categorizing each toxin We identified changes in gene expression after PAH exposure by hybridizing RNA from HepG2 cells to the Human Whole Genome 44K array. First, we evaluated the characteristic molecular profiles for PAH exposure. Total RNA from HepG2 cells was extracted and subjected to DNA microarrays containing approximately 44,000 genetic elements. Generally, transcriptomic regulation occurs through dynamic intracellular responses to multiple, apparently unrelated, stimuli and stresses; thus, characteristic gene expression profiles reflect these stimuli. To investigate whether a large-scale gene expression profile in a cellular system reflects ectopic stress, we performed an unsupervised hierarchical clustering analysis. Of the 44,000 DNA microarray probes, 11,690 genetic elements that passed minimal filtering criteria [analysis of variance (ANOVA); Benjamini and Hochberg false discovery rate (FDR) < 0.05] resulted in the transcriptome of HepG2 cells. From the clustering pattern, as shown in Fig. 2A, some valid clusters were identified, but many ambiguities were observed. Hierarchical clustering analysis of these gene sets largely resulted in three subclusters on a dendrogram, but each toxin induced a different pattern of transcriptomic regulation. As depicted in Fig. 2B, cluster 1 groups, including both carcinogens and noncarcinogens (3MC, DBA, NP, and PH exposure groups), were distinguished from cluster 2 and cluster 3 groups, which included only carcinogens (cluster 2: BA, IND, and BF exposure groups; cluster 3: BP and CHR exposure groups). Moreover, the tree has a cascade shape, which is indicative of a poor correlation, and no clear distinction is made between the main groups of lesions. This suggests that each PAH induces a characteristic molecular profile in the cellular system transcriptomes, and that this profile reflects the differential gene expression based on different toxins. On the basis of this result, we identified outlier genes that can be used to identify different environmental toxins. As shown in Fig. 3A, a supervised analysis using ANOVA (Benjamini and Hochberg FDR < 0.05) and a classification algorithm resulted in 31 outlier genes that were divided exactly into two groups (noncarcinogens vs. carcinogens) based on expression profiles (Fig. 3B). Differences in carcinogenicity-related responses are also visualized by principal components analysis (PCA). Comparable with clustering data, this again shows that 31 genes could clearly separate the carcinogen-exposed groups from the noncarcinogen-exposed control group (Fig. 3C). These outlier genes were further validated by

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

21

Fig. 1. Cell viability based on MTT assay. HepG2 cells were treated for 48 h with BP, DBA, 3MC, NP, CHR, PH, BA, BF, and IDN, and the dose–response curve was plotted. BP, benzo[a]pyrene; DBA, dibenzo[a,h]anthracene; 3MC, 3-methylcholanthrene; NP, naphthalene; CHR, chrysene; PH, phenanthrene; BA, benzo[a]anthracene; BF, benzo[k]fluoranthene; IND, indeno[1,2,3-c,d]pyrene.

prediction confidence analysis, which is classified by the support vector machine (SVM) algorithm using Leave-One-Out Cross Validation (LOOCV). A summary of the assignment frequencies using a high-accuracy classifier (31 genes) is provided in Fig. 3D, with 100% accuracy for the two classes (carcinogens and noncarcinogens). 3.3. Identification of PAH molecular markers to predict carcinogenicity We identified molecular markers that could distinguish between PAHs with different carcinogenicities. To obtain the maximum number of genes that accurately distinguish between each toxin, we performed two different multi-classification algorithms (the SAM method and an ANOVA; Benjamini and Hochberg FDR significance level = 0.005) using the commercial analytical program GeneSpring GX v11.0, as reported by Kim et al. (2011). We then used the Kruskal–Wallis H-test to identify genes that classified each toxin with 100% accuracy, and finally combined the outlier genes from these two methods (Fig. 4A). Initially, we applied SAM with a high stringent cutoff value (q < 5- and 1.5-fold change;

GeneSpring package tool) to each toxin to identify large-scale gene expression changes following exposure to group 1, group 2A, group 2B, group 3, and nonclassified group toxins. Next, we identified genes that successfully classified the five different groups using the Kruskal–Wallis H-test. Then we identified maximum number of genes that satisfied with 100% accuracy by using LOOCV. This selection criterion resulted in 8697 genes as PAH classifiers. Again, we identified each toxin’s characteristic gene signature using another algorithm and one-way ANOVA (Benjamini and Hochberg FDR < 0.005), and identified 995 PAH classifier genes using the Kruskal–Wallis H-test. By combining these outliers, we found that the maximum number of genes that distinguished between each toxin within the PAH exposure group was 430 (Fig. 4A). As shown in Fig. 4B, hierarchical clustering of the expression of 430 genes classified the five carcinogenic toxin groups on the dendrogram based on their own expression profiles. PCA for the 430 genes confirmed that they could clearly separate carcinogenic from noncarcinogenic groups (Fig. 4B). Furthermore, a validation test for the class assignment potential of this prediction marker by SVM using LOOCV resulted in 100% accuracy for all five classes (Fig. 4C).

22

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

Fig. 2. Unsupervised hierarchical clustering analysis and expression profiling of PAH-exposed HepG2 cells. (A) A two-dimensional diagram of the 11,690 genes selected using minimal criteria. Rows represent HepG2 cell gene expression profiles in the presence of various environmental toxins. Columns represent probe measurements. Color saturation reflects differences in expression between PAH-treated RNA and control RNA. (B) Dendrogram derived from clustering using the 11,690 gene set. BP, benzo[a]pyrene; DBA, dibenzo[a,h]anthracene; 3MC, 3-methylcholanthrene; NP, naphthalene; CHR, chrysene; PH, phenanthrene; BA, benzo[a]anthracene; BF, benzo[k]fluoranthene; IND, indeno[1,2,3-c,d]pyrene.

3.4. Validation of molecular markers For the purpose of validating the carcinogen class discrimination models with the final classifier sets, gene expression data were generated for five additional 2B group compounds (not PAHs toxicant; ethylbenzene, dichloromethane, toxaphen, mirex and hexachlorobenzene) and the vehicle control (DMSO) for exposure periods of 48 h. PCA for the 430 classifier genes confirmed that all the independent triplicate treatments of all 2B group compounds

were classified as 2B carcinogen and classified correctly with a predicted test probability of 100% (Fig. 5A andB). 3.5. Functional categorization by gene ontology (GO) and KEGG pathway analysis of discernible characteristic molecular markers that are differentially expressed in PAH-exposed HepG2 cells To clarify the molecular mechanisms of the molecular markers that are differentially expressed in PAH-exposed HepG2 cells, we

Fig. 3. Identifying highly accurate gene classifiers. High-accuracy classifiers that precisely distinguish between PAHs and controls were identified using analysis of variance (Benjamini and Hochberg false discovery rate < 0.05). (A) Characteristic expression profiles of the 31 classifier genes. (B) Results of a principal components analysis (PCA) for the set of 31 genes presented as a two-dimensional scatterplot showing two different classes (carcinogens and noncarcinogens). (C) Accuracy was tested using the LeaveOne-Out Cross Validation (LOOCV) method, and resultantly all nine samples were categorized into two different classes. BP, benzo[a]pyrene; DBA, dibenzo[a,h]anthracene; 3MC, 3-methylcholanthrene; NP, naphthalene; CHR, chrysene; PH, phenanthrene; BA, benzo[a]anthracene; BF, benzo[k]fluoranthene; IND, indeno[1,2,3-c,d]pyrene.

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

23

Fig. 4. Identification of discernible molecular markers for environmental toxins. (A) Combination of the SAM method and a multi-classification analysis test. (B) Results of the principal components analysis (PCA) for the set of 430 genes presented as a three-dimensional scatterplot showing the categories of the five different classes (groups 1, 2A, 2B, 3, and a nonclassified group). (C) Accuracy was tested using the LOOCV method, and resultantly all nine samples were categorized into five different classes. BP, benzo[a]pyrene; DBA, dibenzo[a,h]anthracene; 3MC, 3-methylcholanthrene; NP, naphthalene; CHR, chrysene; PH, phenanthrene; BA, benzo[a]anthracene; BF, benzo[k]fluoranthene; IND, indeno[1,2,3-c,d]pyrene.

evaluated the molecular mechanism of the 430 genes. The use of the GO statistical analysis tool allows for increased biological insight into functional processes. Web-based GO analysis tools were used to identify the GO category. The data were generated with the 430 molecular markers, which were categorized based on GO terms with p-values ≤0.05. Table 2 shows the GO term categories of biological processes, cellular components, and molecular functions obtained from the DAVID 2008 annotation tool. The major biological process ontology categories were cell cycle, proteolysis, protein localization, macromolecular complex subunit organization, oxidation, reduction, and others. The cellular components of the 430 classifier genes were linked to intracellular non-membrane-bound organelles, cytoskeleton, chromosomes, proteasome complexes, and others. Regarding molecular function, transcripts were mainly

linked to peptidase activity, tubulin binding, drug binding, and others (Table 2). The KEGG pathways were compiled from multiple literature sources and integrated from individual components into a unified pathway (Kanehisa et al., 2006). Therefore, the KEGG pathway database characterized the enrichment of specific pathway components into functionally regulated gene groups. Fig. 6 shows the number of genes with altered expression patterns in the presence of individual toxins in each KEGG pathway. The major genes were involved in proteasome, spliceosome, glutathione metabolism, drug metabolism, and metabolism of xenobiotics by cytochrome P450, and phenylalanine metabolism. GO and KEGG pathway database analyses provided a valuable mechanistic insight into PAH exposure in human cellular systems (Table 3).

24

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

Fig. 5. Validation of the class prediction model with five additional compounds. (A) Results of a principal components analysis (PCA) for the set of 430 genes presented as a three-dimensional scatterplot showing the classification of the five test compounds. (B) Accuracy was tested using the Leave-One-Out Cross Validation (LOOCV) method, and resultantly all five test compounds were categorized as 2B carcinogen.

3.6. Validation of molecular markers that are differentially expressed in PAH-exposed HepG2 cells To validate microarray gene expression data, and to confirm the transcriptional levels of differentially expressed genes, we selected genes that were significantly affected by the different toxins (BP,

DBA, 3MC, NP, CHR, PH, BA, BF, and IND) based on the microarray data, and performed quantitative RT-PCR with SYBR Green fluorescent dye. As shown in Fig. 7, transcriptional levels of four genes (IER3, TIMP1, LGALS1, and RHOC) that were significantly upregulated in all nine toxin exposure groups, and another two genes (NDC80 and PAH) that were significantly downregulated in all nine

Fig. 6. Functional classification of 430 predictive molecular markers using the KEGG pathway database.

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

25

Fig. 7. Comparison of DNA microarray and quantitative real-time RT-PCR data for the PAH classifier genes. mRNA levels were analyzed by quantitative real-time RTPCR. BP, benzo[a]pyrene; DBA, dibenzo[a,h]anthracene; 3MC, 3-methylcholanthrene; NP, naphthalene; CHR, chrysene; PH, phenanthrene; BA, benzo[a]anthracene; BF, benzo[k]fluoranthene; IND, indeno[1,2,3-c,d]pyrene.

toxin groups, were compared to their corresponding controls. The GAPDH gene was used as an endogenous control for all genes. The expression of IER3, TIMP1, LGALS1, and RHOC genes based on microarray data were upregulated by 1.52–22.72-fold, while NDC80 and PAH were downregulated by 0.043–0.624-fold. The quantitative RT-PCR analyses of these genes confirmed the underor overexpression, and corroborated the microarray data. 4. Discussion This study illustrated the successful classification of RNA samples derived from the liver cells of various carcinogenic chemical exposure. Also, this study also shows the development of a rapid, highly accurate, and widely applicable method to predict chemical carcinogenicity based on gene expression profiles in the presence of different chemicals. The success of our predictions demonstrates that if one can reduce the dimension of their data set with various statistical models, then it is possible to perform predictions for unknown compounds that are similar to compounds present in database. We are now using the data set to develop more formal and rigorous models for prediction to make this task less arbitrary in future studies. The use of biological markers to directly measure the levels of suspected environmental chemicals in human tissues and fluids, and linking these exposures with disability or disease, is becoming more common (Perera, 1997). Toxicogenomics, which is used interchangeably with the term “environmental genomics” (Aardema and MacGregor, 2002; Irwin et al., 2004; Nuwaysir et al., 1999), is an emerging toxicology discipline that enables scientists to identify and characterize the molecular signatures of environmental toxins. Multiple in vitro models have been used to identify a toxicogenomics-based test assessing the toxicity or carcinogenicity of toxins as an alternative to rodent models. A previous report by Eun et al. (2008) suggested that characteristic molecular signatures of toxins in liver tissues could reveal the

molecular basis of hepatotoxicity after 72 h exposure. This study found that the transcriptomic response to chemical exposure was specific to the toxin, and this characteristic molecular signature could be used as an early prognostic biomarker. Hepatic models are frequently used for toxicogenomic studies since the liver plays a major role in the metabolism of many compounds and also represents an important target organ in systemic toxicity (Davila et al., 1998). A preferred model is the human hepatoma cell line HepG2. HepG2 cells biotransform xenobiotic compounds and are capable of activating mutagens and carcinogens. They also carry no p53 mutations, which enables the cells to activate DNA damage response, induce growth arrest, and initiate apoptosis (Hsu et al., 1993; Knasmuller et al., 1998; Wilkening et al., 2003). This cell line has frequently been applied in toxicogenomics studies (Burczynski et al., 2000; Harries et al., 2001; Hong et al., 2003; Van Delft et al., 2005), as HepG2 cells are easy to handle compared to, for example, primary human hepatocytes. Microarrays use immobilized cDNA or oligonucleotide probes to simultaneously monitor the expression of thousands of genes and obtain information on global gene expression (i.e., a view of all mRNA transcripts expressed by a cell is known as the transcriptome) (Staudt, 2003; Staudt and Brown, 2000) and are becoming increasingly common in toxicology (Hamadeh et al., 2002a,b). We hypothesized that microarrays could identify changes in gene expression, and that certain genes could act as surrogate molecular signatures of exposure, revealing the early effects of environmental toxins such as PAHs, and provide information on the carcinogenic potencies and toxicity mechanisms. In this study, we applied transcriptomic expression technology to identify and characterize a molecular signature in a cellular system, and to identify discernible and predictive characteristic markers induced by PAHs with different carcinogenicities. We found that exposure to PAHs altered global gene expression in HepG2 cells. Transcriptomic responses in HepG2 cells appeared to be toxin-specific and displayed characteristic molecular

26

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

Table 2 GO functional category analysis of 430 classifier genes. Term (A) Biological Process (top 10 out of 104 total) Cell cycle Proteolysis Macromolecule catabolic process Protein localization Macromolecular complex subunit organization Positive regulation of macromolecule metabolic process Macromolecular complex assembly Oxidation reduction RNA processing Negative regulation of apoptosis (B) Cellular component (top 10 out of 24 total) Intracellular non-membrane-bounded organelle Cytosol Cytoskeleton Ribonucleoprotein complex Chromosomal part Microtubule organizing center Proteasome complex Chromosome, centromeric region Spindle Condensed chromosome kinetochore (C) Molecular function (top 10 out of 14 total) Peptidase activity Enzyme activator activity Translation factor activity, nucleic acid binding Tubulin binding Drug binding Ribonucleoprotein binding Ribosome binding Oxidoreductase activity, acting on the CH–CH group of donors, NAD or NADP as acceptor Galactoside binding Glutamate–cysteine ligase activity

Table 3 Functional classification and gene list of molecular signature of PAHs in HepG2 cells.

Genesa

p-Value*

[KEGG pathway] Gene symbol

GenBank accession no.

Description

41 37 31 31 27

5.63E−09 3.56E−04 1.62E−04 0.0012 9.16E−04

[Proteasome] PSMD1

NM 002807

PSMB5

NM 002797

PSMC1

NM 002802

26

0.0199 PSMB3

NM 002795

25 21 19 15

0.0017 0.0183 0.0156 0.0073

PSMA1

NM 148976

PSMC2

NM 002803

69

0.0026

PSMA5

NM 002790

49 37 18 17 13 9 9 9 8

9.84E−06 0.0326 0.0201 0.0031 0.0034 1.93E−05 0.0026 0.0072 1.11E−04

PSMD4

NM 002810

PSMB1

NM 002793

Proteasome (prosome, macropain) 26S subunit, non-ATPase, 1 Proteasome (prosome, macropain) subunit, beta type, 5 Proteasome (prosome, macropain) 26S subunit, ATPase, 1 Proteasome (prosome, macropain) subunit, beta type, 3 Proteasome (prosome, macropain) subunit, alpha type, 1 Proteasome (prosome, macropain) 26S subunit, ATPase, 2 Proteasome (prosome, macropain) subunit, alpha type, 5 Proteasome (prosome, macropain) 26S subunit, non-ATPase, 4 Proteasome (prosome, macropain) subunit, beta type, 1

[Spliceosome] SFRS7

NM 001031684

SFRS1

NM 006924

19 13 7

0.0225 0.0240 0.0102

CDC40

NM 015891

PRPF38A

NM 032864

6 5 4 3 3

0.0396 0.0218 0.0253 0.0438 0.0484

SFRS3

NM 003017

SF3B3

NM 012426

PHF5A

NM 032758 NM 001011724

2 2

0.0371 0.0371

[Glutathione metabolism] RRM2 NM 001034

DAVID 2008 functional annotation bioinformatics microarray analysis software was used to obtain the GO functional categories. a Some genes are counted in more than one annotation category. * Statistical significant GO terms were only considered (p-value ≤ 0.05).

signatures (Fig. 2), which indicates that transcriptomic regulation in human cellular systems reflects xenobiotic stresses to the body. Also, all compounds can be divided into three groups: cluster 1 consisting of 3MC, DBA, NP, and PH; cluster 2 consisting of BA, IND, and BF; and cluster 3 consisting of CHR and BP. The first group includes both carcinogenic and noncarcinogenic PAHs, while groups 2 and 3 contain only carcinogenic or suspected carcinogenic PAHs. Thus, these three groups do not perfectly correspond with toxin carcinogenicity based on studies carried out in HepG2 cells. In addition, a simple nonparametric analysis resulted in 31 genes that distinguished carcinogens from noncarcinogens, regardless of the toxin (Fig. 3). All carcinogenic PAHs modified genes involved with oxidative stress (which implies an increased antioxidant production and decreased ROS formation), whereas the least carcinogenic PAHs, fluoranthene (FA) and 1-methylphenanthrene (1-MPA), did not. This suggests that affecting oxidative stress may be an important characteristic or classifier for carcinogenic PAHs (Staal et al., 2007). Also, even if individual PAH exposure is low (not carcinogenic), the combination of different PAHs may lead to the induction and increased activity of enzyme systems and increase oxidative stress (Alessio, 1996). Thus, the 31 outlier genes may represent a common genetic and epigenetic response associated with carcinogen-specific stress. This list of genes may provide a molecular signature that could discriminate carcinogenic from noncarcinogenic PAHs.

GCLM

NM 002061

GCLC

NM 001498

GSR GSTA4 GSTA3 GSTA2 [Drug metabolism] MAOA MAOB GSTA4 ADH5

NM NM NM NM

000637 001512 000847 000846

NM NM NM NM

000240 000898 001512 000671

Splicing factor, arginine/serine-rich 7, 35 kDa Splicing factor, arginine/serine-rich 1 (splicing factor 2, alternate splicing factor) Cell division cycle 40 homolog (S. cerevisiae) PRP38 pre-mRNA processing factor 38 (yeast) domain containing A Splicing factor, arginine/serine-rich 3 Splicing factor 3b, subunit 3, 130 kDa PHD finger protein 5A Heterogeneous nuclear ribonucleoprotein A1-like Ribonucleotide reductase M2 polypeptide Glutamate–cysteine ligase, modifier subunit Glutamate–cysteine ligase, catalytic subunit Glutathione reductase Glutathione S-transferase A4 Glutathione S-transferase A3 Glutathione S-transferase A2

Monoamine oxidase A Monoamine oxidase B Glutathione S-transferase A4 Alcohol dehydrogenase 5 (class III), chi polypeptide NM 000847 Glutathione S-transferase A3 GSTA3 NM 000846 Glutathione S-transferase A2 GSTA2 [Metabolism of xenobiotics by cytochrome P450] NM 001353 Aldo-keto reductase family 1, AKR1C1 member C1 (dihydrodiol dehydrogenase 1; 20-alpha (3-alpha)-hydroxysteroid dehydrogenase) NM 001512 Glutathione S-transferase A4 GSTA4 NM 000671 Alcohol dehydrogenase 5 (class III), ADH5 chi polypeptide Glutathione S-transferase A3 GSTA3 NM 000847 GSTA2 NM 000846 Glutathione S-transferase A2 [Phenylalanine metabolism] NM 000240 Monoamine oxidase A MAOA NM 003491 ARD1 homolog A, ARD1A N-acetyltransferase (S. cerevisiae) Phenylalanine hydroxylase PAH NM 000277 NM 000898 Monoamine oxidase B MAOB

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

In the field of environmental toxicology, the potential application of toxicogenomics to predict toxicity and carcinogenicity of unknown compounds remains largely unexplored, although studies showing the advantages of microarrays to evaluate and assess toxin exposure are steadily increasing. The evaluation of gene expression upon chemical exposure provides knowledge on the toxicity mechanism and reveals a “molecular signature” based on the gene expression changes both in vitro (Burczynski et al., 2000; Waring et al., 2001) and in vivo (Hamadeh et al., 2002b). The development of genomic signatures may allow for early screening of suspected toxins based on their similarity to other known compounds. Previous studies distinguished between carcinogens and noncarcinogens, or genotoxic carcinogens and non-genotoxic carcinogens, using gene expression profiling (Matsumoto et al., 2009; Mathijs et al., 2009). However, only a few studies have used this approach to examine compounds with carcinogenic potency according to IARC classification. To investigate this issue, we proposed two different serial multi-classification analyses and combined them to evaluate the toxicity of each PAH class. The combination of a one-way ANOVA and SAM analysis resulted in 430 genes that exactly annotated five different classes with 100% accuracy (Fig. 4). This suggests that HepG2 cells are a suitable model to classify the carcinogenic potency of PAHs based on gene expression profiling. However, we tested only nine compounds, so a larger number of compounds may improve the classification of PAHs in human cellular systems. Furthermore, the underlying molecular pathways and interaction of molecular markers are very important because the development of toxicity or disease is a complex process involving many interactive molecules and pathways. As expected, the underlying molecular pathways in the 430 prediction markers mainly included metabolism related pathways, which are involved in intracellular signaling pathways in liver cells and the carcinogenic effect of toxins (Fig. 6). Based on these results, we propose that pathway analyses can be used to elucidate the mechanism of various environmental toxins. Our results demonstrate that characteristic molecular signatures of PAHs exist in human liver cells, and that a multi-classification of this PAH-specific signature is able to predict carcinogenic potency. This molecular signature could be used as a predictable and discernible surrogate marker to determine transcriptomic responses to PAH exposure in the environment. In conclusion, toxicogenomics applied to in vitro carcinogenicity assays identified a gene signature that is likely predictive of chemical carcinogens. Although only a small number of chemicals were tested, our observations warrant further study to validate this alternative testing strategy. Also, further characterization of these gene signatures is required to provide mechanistic insight into the toxicological effects of PAHs in human cellular systems. Conflict of interest statement Authors declare that there are no conflicts of interest. Acknowledgement This work was supported by the Korea Research Foundation grants from Korea Ministry of Environment (grant number 212101-002) as “The Converging Technology”. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.toxlet. 2012.04.013.

27

References Aardema, M.J., MacGregor, J.T., 2002. Toxicology and genetic toxicology in the new era of “toxicogenomics”: impact of “-omics” technologies. Mutation Research 499, 13–25. Alessio, L., 1996. Multiple exposure to solvents in the workplace. International Archives of Occupational and Environmental Health 69, 1–4. Boffetta, P., Jourenkova, N., Gustavsson, P., 1997. Cancer risk from occupational and environmental exposure to polycyclic aromatic hydrocarbons. Cancer Causes and Control 8, 444–472. Burczynski, M.E., McMillian, M., Ciervo, J., Li, L., Parker, J.B., Dunn, R.T., Hicken, I.I., Farr, S., Johnson, S.M.D., 2000. Toxicogenomics-based discrimination of toxic mechanism in HepG2 human hepatoma cells. Toxicological Sciences 58, 399–415. Butterworth, B.E., Smith-Oliver, T., Earle, L., Loury, D.J., White, R.D., Doolittle, D.J., Working, P.K., Cattley, R.C., Jirtle, R., Michalopoulos, G., Strom, S., 1989. Use of primary cultures of human hepatocytes in toxicology studies. Cancer Research 49, 1075–1084. Choi, H.S., Kim, Y.J., Song, M., Song, M.K., Ryu, J.C., 2011. Identification of hepatotoxicity related genes induced by toxaphene in HepG2 cells. Molecular & Cellular Toxicology 7, 53–60. Cohen, S.M., 2004. Human carcinogenic risk evaluation: an alternative approach to the two-year rodent bioassay. Toxicological Sciences 80, 225–229. Cunningham, M.L., 1996. Role of increased DNA replication in the carcinogenic risk of nonmutagenic chemical carcinogens. Mutation Research 365, 59–69. Dambach, D.M., Andrews, B.A., Moulin, F., 2005. New technologies and screening strategies for hepatotoxicity: use of in vitro models. Toxicologic Pathology 33, 17–26. Davila, J.C., Rodriguez, R.J., Melchert, R.B., Acosta Jr., D., 1998. Predictive value of in vitro model systems in toxicology. Annual Review of Pharmacology and Toxicology 38, 63–96. Ellinger-Ziegelbauer, H., Gmuender, H., Bandenburg, A., Ahr, H.J., 2008. Prediction of a carcinogenic potential of rat hepatocarcinogens using toxicogenomics analysis of short-term in vivo studies. Mutation Research 637, 23–39. Eun, J.W., Ryu, S.Y., Noh, J.H., Lee, M.J., Jang, J.J., Ryu, J.C., Jung, K.H., Kim, J.K., Bae, H.J., Xie, H., Kim, S.Y., Lee, S.H., Park, W.S., Yoo, N.J., Lee, J.Y., Nam, S.W., 2008. Discriminating the molecular basis of hepatotoxicity using the large-scale characteristic molecular signatures of toxicants by expression profiling analysis. Toxicology 249, 176–183. Groneberg, D.A., Grosse-Siestrup, C., Fischer, A., 2002. In vitro models to study hepatotoxicity. Toxicologic Pathology 30, 394–399. Hamadeh, H.K., Amin, R.P., Paules, R.S., Afshari, C.A., 2002a. An overview of toxicogenomics. Current Issues in Molecular Biology 4, 45–56. Hamadeh, H.K., Bushel, P.R., Jayadev, S., Martin, K., DiSorbo, O., Sieber, S., Bennett, L., Tennant, R., Stoll, R., Barrett, J.C., Blanchard, K., Paules, R.S., Afshari, C.A., 2002b. Gene expression analysis reveals chemical-specific profiles. Toxicological Sciences 67, 219–231. Harries, H.M., Fletcher, S.T., Duggan, C.M., Baker, V.A., 2001. The use of genomics technology to investigate gene expression changes in cultured human liver cells. Toxicology In Vitro 15, 399–405. Hong, Y., Muller, U.R., Lai, F., 2003. Discriminating two classes of toxicants through expression analysis of HepG2 cells with DNA arrays. Toxicology In Vitro 17, 85–92. Hsu, I.C., Tokiwa, T., Bennett, W., Metcalf, R.A., Welsh, J.A., Sun, T., Harris, C.C., 1993. p53 gene mutation and integrated hepatitis B viral DNA sequences in human liver cancer cell lines. Carcinogenesis 14, 987–992. Hu, T., Gibson, D.P., Carr, G.J., Torontali, S.M., Tiesman, J.P., Chaney, J.G., Aardema, M.J., 2004. Identification of a gene expression profile that discriminates indirectacting genotoxins from direct-acting genotoxins. Mutation Research 549, 5–27. Irwin, R.D., Boorman, G.A., Cunningham, M.L., Heinloth, A.N., Malarkey, D.E., Paules, R.S., 2004. Application of toxicogenomics to toxicology: basic concepts in the analysis of microarray data. Toxicologic Pathology 32, 72–83. Jacobs, A., Jacobson-Kram, D., 2004. Human carcinogenic risk evaluation. Part III. Assessing cancer hazard and risk in human drug development. Toxicological Sciences 81, 260–262. Kanehisa, M., Goto, S., Hattori, M., Aoki-Kinoshita, K.F., Itoh, M., Kawashima, S., Katayama, T., Araki, M., Hirakawa, M., 2006. From genomics to chemical genomics: new developments in KEGG. Nucleic Acids Research 34, 354–357. Kim, J.K., Jung, K.H., Noh, J.H., Eun, J.W., Bae, H.J., Xie, H.J., Jang, J.J., Ryu, J.C., Park, W.S., Lee, J.Y., Nam, S.W., 2011. Identification of characteristic molecular signature for volatile organic compounds in peripheral blood of rat. Toxicology and Applied Pharmacology 250, 162–169. Knasmuller, S., Parzefall, W., Sanyal, R., Ecker, S., Schwab, C., Uhl, M., MerschSundermann, V., Williamson, G., Hietsch, G., Langer, T., Darroudi, F., Natarajan, A.T., 1998. Use of metabolically competent human hepatoma cells for the detection of mutagens and antimutagens. Mutation Research 402, 185–202. Kushida, T., Takagi, T., Fukuda, K.I., 2006. Event ontology: a pathway-centric ontology for biological processes. Pacific Symposium on Biocomputing 11, 152–163. Mathijs, K., Brauers, K.J., Jennen, D.G., Boorsma, A., van Herwijnen, M.H., Gottschalk, R.W., Kleinjans, J.C., van Delft, J.H., 2009. Discrimination for genotoxic and nongenotoxic carcinogens by gene expression profiling in primary mouse hepatocytes improves with exposure time. Toxicological Sciences 112, 374–384. Matsumoto, H., Yakabe, Y., Saito, K., Sumida, K., Sekijima, M., Nakayama, K., Miyaura, H., Saito, F., Otsuka, M., Shirai, T., 2009. Discrimination of carcinogens by hepatic transcript profiling in rats following 28-day administration. Cancer Inform 13, 253–269.

28

M.-K. Song et al. / Toxicology Letters 212 (2012) 18–28

Mosmann, T., 1983. Rapid colorimetric assay for cellular growth and survival: application to proliferation and cytotoxicity assays. Journal of Immunological Methods 65, 55–63. Nesnow, S., Mass, M.J., Ross, J.A., Galati, A.J., Lambert, G.R., Gennings, C., Carter Jr., W.H., Stoner, G.D., 1998. Lung tumorigenic interactions in strain A/J mice of five environmental polycyclic aromatic hydrocarbons. Environmental Health Perspectives 106, 1337–1346. Nuwaysir, E.F., Bittner, M., Trent, J., Barrett, J.C., Afshari, C.A., 1999. Microarrays and toxicology: the advent of toxicogenomics. Molecular Carcinogenesis 24, 153–159. Perera, F.P., 1997. Environment and cancer: who are susceptible? Science 278, 1068–1073. Song, M.K., Kim, Y.J., Song, M., Ryu, J.C., 2011a. Gene expression analysis identifies potential biomarkers of phenanthrene in human hepatocytes (HepG2). Toxicology and Environmental Health Sciences 3, 30–38. Song, M.K., Kim, Y.J., Song, M., Choi, H.S., Park, Y.K., Ryu, J.C., 2011b. Polycyclic aromatic hydrocarbons induce migration in human hepatocellular carcinoma cells (HepG2) through reactive oxygen species-mediated p38 MAPK signal transduction. Cancer Science 102, 1636–1644. Spalding, J.W., French, J.E., Stasiewicz, S., Furedi-Machacek, M., Conner, F., Tice, R.R., Tennant, R.W., 2000. Responses of transgenic mouse lines p53(+/−) and Tg.AC to agents tested in conventional carcinogenicity bioassays. Toxicological Sciences 53, 213–223. Staal, Y.C., van Herwijnen, M.H., Pushparajah, D.S., Umachandran, M., Ioannides, C., van Schooten, F.J., van Delft, J.H., 2007. Modulation of gene expression and DNAadduct formation in precision-cut liver slices exposed to polycyclic aromatic hydrocarbons of different carcinogenic potency. Mutagenesis 22, 55–62. Staudt, L.M., 2003. Molecular diagnosis of the hematologic cancers. New England Journal of Medicine 348, 1777–1785.

Staudt, L.M., Brown, P.O., 2000. Genomic views of the immune system. Annual Review of Immunology 18, 829–859. Te Pas, M.F., Hulsegge, I., Coster, A., Pool, M.H., Heuven, H.H., Janss, L.L., 2007. Biochemical pathways analysis of microarray results: regulation of myogenesis in pig. BMC Developmental Biology 7, 66–80. Tennant, R.W., Elwell, M.R., Spalding, J.W., Griesemer, R.A., 1991. Evidence that toxic injury is not always associated with induction of chemical carcinogenesis. Molecular Carcinogenesis 4, 420–440. Tusher, V.G., Tibshirani, R., Chu, G., 2001. Significance analysis of microarrays applied to the ionizing radiation response. Proceedings of the National Academy of Sciences of the United States of America 98, 5116–5121. Van Delft, J.H., van Agen, E., van Breda, S.G., Herwijnen, M.H., Staal, Y.C., Kleinjans, J.C., 2004. Discrimination of genotoxic from non-genotoxic carcinogens by gene expression profiling. Carcinogenesis 25, 1265–1276. Van Delft, J.H., van Agen, E., van Breda, S.G., Herwijnen, M.H., Staal, Y.C., Kleinjans, J.C., 2005. Comparison of supervised clustering methods to discriminate genotoxic from non-genotoxic carcinogens by gene expression profiling. Mutation Research 575, 17–33. Wang, J.S., Busby Jr., W.F., 1993. Induction of lung and liver tumors by fluoranthene in a preweanling CD-1 mouse bioassay. Carcinogenesis 14, 1871–1874. Waring, J.F., Ciurlionis, R., Jolly, R.A., Heindel, M., Ulrich, R.G., 2001. Microarray analysis of hepatotoxins in vitro reveals a correlation between gene expression profiles and mechanisms of toxicity. Toxicology Letters 120, 359–368. Wilkening, S., Stahl, F., Bader, A., 2003. Comparison of primary human hepatocytes and hepatoma cell line Hepg2 with regard to their biotransformation properties. Drug Metabolism and Disposition: The Biological Fate of Chemicals 31, 1035–1042.