Toxicology 328 (2015) 29–39

Contents lists available at ScienceDirect

Toxicology journal homepage: www.elsevier.com/locate/toxicol

A novel transcriptomics based in vitro method to compare and predict hepatotoxicity based on mode of action K. Nadira De Abrew *, Gary J. Overmann, Rachel L. Adams, Jay P. Tiesman, John Dunavent, Yuqing K. Shan, Gregory J. Carr, George P. Daston, Jorge M. Naciff Mason Business Center, The Procter & Gamble Company, Cincinnati, OH 45040, USA

A R T I C L E I N F O

A B S T R A C T

Article history: Received 15 July 2014 Received in revised form 31 October 2014 Accepted 29 November 2014 Available online 2 December 2014

High-content data have the potential to inform mechanism of action for toxicants. However, most data to support this notion have been generated in vivo. Because many cell lines and primary cells maintain a differentiated cell phenotype, it is possible that cells grown in culture may also be useful in predictive toxicology via high-content approaches such as whole-genome microarray. We evaluated global changes in gene expression in primary rat hepatocytes exposed to two concentrations of ten hepatotoxicants: acetaminophen (APAP), b-naphthoflavone (BNF), chlorpromazine (CPZ), clofibrate (CLO), bis(2-ethylhexyl)phthalate (DEHP), diisononyl phthalate (DINP), methapyrilene (MP), valproic acid (VPA), phenobarbital (PB) and WY14643 at two separate time points. These compounds were selected to cover a range of mechanisms of toxicity, with some overlap in expected mechanism to address the question of how predictive gene expression analysis is, for a given mode of action. Gene expression microarray analysis was performed on cells after 24 h and 48 h of exposure to each chemical using Affymetrix microarrays. Cluster analysis suggests that the primary hepatocyte model was capable of responding to these hepatotoxicants, with changes in gene expression that appear to be mode of actionspecific. Among the different methods used for analysis of the data, a combination method that used pathways (MOAs) to filter total probesets provided the most robust analysis. The analysis resulted in the phthalates clustering closely together, with the two other peroxisome proliferators, CLO and WY14643, eliciting similar responses at the whole-genome and pathway levels. The Cyp inducers PB, MP, CPZ and BNF also clustered together. VPA and APAP had profiles that were unique. A similar analysis was performed on externally available (TG-GATES) in vivo data for 6 of the chemicals (APAP, CLO, CPZ, MP, MP and WY14643) and compared to the in vitro result. These results indicate that transcription profiling using an in vitro assay may offer pertinent biological data to support predictions of in vivo hepatotoxicity potential. ã 2014 Elsevier Ireland Ltd. All rights reserved.

Keywords: Toxicogenomics Hepatotoxicity 21st century Tox Transcriptomics Toxicity pathways

1. Introduction Changes in gene expression following exposure to a xenobiotic provide a characteristic pattern that is diagnostic of mode of action (MOA). The measurement of genome-wide mRNA changes or transcriptomics, to assess the toxicological profile of chemicals was first introduced in 1999 (Nuwaysir et al., 1999) and since then been applied to various areas of toxicological relevance (Chen et al., 2012; Daston and Naciff, 2010). Demonstration of the technology’s effectiveness and predictive power has led to a general expectation

* Corresponding author. Tel.: +1 513 626 4067. E-mail address: [email protected] (K. N. De Abrew). http://dx.doi.org/10.1016/j.tox.2014.11.008 0300-483X/ ã 2014 Elsevier Ireland Ltd. All rights reserved.

that the wider field of toxicogenomics will enhance our ability to do risk assessment, ultimately leading to better predictive toxicology (Goetz et al., 2011; NRC, 2007a; van Leeuwen et al., 2011). Pathway analysis provides a current literature-based, biologically relevant method to identify putative mode of action (van Leeuwen et al., 2011). In vivo studies have shown that well-characterized hepatotoxicants change gene expression in a manner that is consistent with what is known about the toxicity of the chemical (Bulera et al., 2001,c; Hamadeh et al., 2002a,b,c,c; Thomas et al., 2001). There is limited information (Waring et al., 2001) showing that primary hepatocytes are also responsive. Other studies have compared human hepatocyte cell lines to human primary hepatocytes (Gerets et al., 2012; Jetten et al., 2013). These studies show that

30

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

different hepatic cell lines have varying capacity to mimic primary cell behavior. Among the studies done to date a majority limited their analysis to a subset of the genes (Bulera et al., 2001; Gerets et al., 2012,c; Hamadeh et al., 2002a,b,c,c; Thomas et al., 2001). This was either due to the use of partial genome arrays (the technology available at the time) or due to the hand-picking of genes thought to play a significant part in the relevant pathways of interest. Previous studies used only a few chemicals making it difficult to ascertain the usefulness of gene expression as a tool to perform MOA based clustering for a wide range of chemicals. The use of MOA as the primary means for comparison between two chemicals is regaining widespread attention and both the National Academy of Sciences (NAS) long range vision for toxicity testing (NRC, 2007b) and the NAS report on application of toxicogenomic technologies for predictive toxicology (NRC, 2007a) outline a risk assessment paradigm that is based on “toxicity pathways” which is essentially a more detailed MOA. In addition, OECD through the QSAR Toolbox project is focusing on Adverse Outcome Pathways (AOP) as the most promising path forward for grouping chemicals for SAR (OECD, 2012). The objective of the present study was to evaluate the possibility of using primary hepatocytes as an in vitro model to predict in vivo gene expression related toxicological MOA for a specific group of chemicals. We believe that chemical mediated pathway perturbations (toxicological MOAs) can be grouped in to two main areas: (1) more specific pathway perturbations mediated through chemical binding to specific receptors or enzymes and (2) more general pathway perturbations mediated through gross binding of chemicals to macromolecules (DNA and protein). While many in vitro toxicogenomic studies for the latter group exists to date (Jetten et al., 2013; Schaap et al., 2012) the literature on the former is sparse. We chose to compare the gene expression profiles of primary rat hepatocytes exposed to 10 different hepatotoxicants within the former group.

An analysis of variance model was used to identify genes exhibiting differential expression vs. control for any of the toxicants. A gene set (pathway) enrichment analysis (using GO terms) was also performed. The results showed that a method selecting for both statistically significant gene-level differential expression and gene set enrichment provides the most robust analysis for clustering the chemicals based on MOA. The study provided strong evidence for using primary cell models as an alternate for whole animal based genomic studies. 2. Materials and methods 2.1. Chemicals and reagents Acetaminophen (99.0%, Catalog no. A7085), bis(2-ethylhexyl) phthalate (99.0%, Catalog no. D201154), valproic acid sodium salt (Catalog no. P4543), methapyrilene HCl (Catalog no. 442641), phenobarbital (Catalog no. P1636), b-naphthoflavone (95%, Catalog no. N3633), chlorpromazine hydrochloride (Catalog no. C8138), diisononyl phthalate (Catalog no. 376663), clofibrate (Catalog no. C6643) and WY14643 (Catalog no. C7081) were all purchased from Sigma–Aldrich (St. Louis, MO). 2.2. Dose selection The highest dose for each chemical treatment was obtained from existing literature values and was predominantly based on in vitro studies (Table 1). Extrapolation to cell culture followed simple calculations that tried to mimic peak concentrations that were not toxic to the cells. A total of five concentrations were used, decreasing from the highest in vitro dose to provide a reasonable dose response coverage. However, only two doses were picked for gene array analysis due to cost considerations. The mid and vhigh doses were picked for this purpose, since these two doses were

Table 1 Chemicals, corresponding dose treatments and relevant literature used for dose selection. Chemical

vLow

Acetaminophen

10 mM

Di(2-ethylhexyl) phthalate

Sodium valproate

Clofibrate Phenobarbital

Beta-naphthoflavone Methapyrilene HCl

Chlorpromazine

WY14643 Diisononyl phthalate

Low

Mid

High

vHigh

Notes

100 mM 500 mM 1000 mM 10,000 mM Tamura et al. (2006) evaluated 300, 1000 and 3000 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 10,000 uM in rat primary hepatocytes. 10 mM 100 mM 250 mM 500 mM 1000 mM Goyak et al. (2008) evaluated 100 uM in human primary hepatocytes. Choi et al. (2010) evaluated DEHP in HepG2 cells at 2.5, 5, 10, 25, 50, 100, and 250 uM for 24 or 48 h. 5 mM 50 mM 500 mM 1000 mM 10,000 mM Tamura et al. (2006) evaluated 400, 2000 and 10,000 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 5000 uM in rat primary hepatocytes.. 1 mM 12 mM 60 mM 300 mM 1000 mM Tamura et al. (2006) evaluated 12, 60 and 300 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 300 uM in rat primary hepatocytes. 10 mM 100 mM 300 mM 1000 mM 3000 mM Tamura et al. (2006) evaluated 300, 1000 and 3000 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 10,000 uM. in rat primary hepatocytes Goyak et al. (2008) evaluated 500 uM in human primary hepatocytes. Gerets et al. (2012) evaluated500 uM in human primary hepatocytes. Jigorel et al. (2006) evaluated 100, 500, 1000, 2000 and 3200 uM, in human primary hepatocytes. 10 nM 100 nM 1 mM 10 mM 100 mM Gerets et al. (2012) evaluated 25 uM in human primary hepatocytes. 100 nM 0.6 mM 3 mM 15 mM 100 mM Tamura et al. (2006) evaluated 0.6, 3, and 15 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 600 uM in rat primary hepatocytes.. 10 nM 100 nm 0.8 mM 4 mM 20 mM Tamura et al. (2006) evaluated 0.6, 3, and 15 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 600 uM. in rat primary hepatocytes. 100 nM 1 mM 8 mM 40 mM 200 mM Tamura et al. (2006) evaluated 8, 40 and 200 uM in rat primary hepatocytes. Kiyosawa et al. (2006) evaluated 150 uM. in rat primary hepatocytes 10 mM 100 mM 250 mM 500 mM 1000 mM Hasmall et al. (1999) evaluated at 250, 500 and 750 uM in rat and human primary hepatocytes.

Grayed columns depict doses used for gene array analysis.

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

expected to provide the most significant gene expression changes while still providing a correlation to dose (dose response) (Table 1). 2.3. Cell culture Frozen male Sprague Dawley primary rat hepatocytes were purchased from XenoTech (Lenexa, KS) and stored in the vapor phase of liquid N2, until use. Cells were thawed using supplier’s protocol and XenoTech cell isolation kit (Lenexa, KS). Briefly, one, 4.5 ml cryovial containing 17  106 rat hepatocytes was removed from liquid N2 and immediately placed in a shaking 37  C water bath for 2 min. Hepatocytes were transferred to a 50 ml conical tube containing 40 ml of solution A (XenoTech Dulbecco’s Modified Eagles Media supplemented with isotonic Percoll, no pen/strep). The cryovial was rinsed with 4.5 ml of Solution B (XenoTech Dulbecco’s Modified Eagles Media supplemented with isotonic Percoll, no pen/strep) and transferred to the 50 ml conical tube containing cells in solution A. Hepatocytes were centrifuged at 60  g (1000 rpm) using an International Equipment Company Centra-8R Centrifuge (Chattanooga, TN) for 5 min at room temperature. Following centrifugation the supernatant was aspirated and hepatocytes were resuspended in 35 ml of Solution B and centrifuged at 60  g (1000 rpm) for another 3 min at room temperature. The resulting cell pellet was resuspended in 30 ml of XenoTech resuspension media (Lenexa, KS) and cell viability was evaluated by trypan blue exclusion. Following the assessment of viability, cells were resuspended in Modified Chee’s culture media (XenoTech, Lenexa, KS) at a concentration of 0.7–1.0  106 cells/ml, and 2 ml of cell suspension was transferred to each well of 6 well collagen coated plates (Corning Corning, NY) and placed in a 37  C incubator (5% CO2, 80% humidity). Following 3–4 h of incubation to allow hepatocytes to attach, media was removed and 2 ml of ice cold modified Chee’s culture media (XenoTech, Lenexa, KS) supplemented with 0.25 mg/ml BD Matrigel, (BD Biosciences, Bedford, MA) was added to each well using pipets chilled to 4  C. Hepatocytes were incubated at 37  C, 5% CO2, media was replaced every day for the next 2 days with Modified Chee’s culture media (XenoTech, Lenexa, KS) and dosed with chemical on day 3. Fresh 100X chemical dosing solutions in DMSO or H2O were prepared 24 h before exposure for all chemical doses (Table 2). Cells were treated with two concentrations of chemical (Table 2) or vehicle

31

(0.1% DMSO) on the morning of day 3 and harvested for RNA at 24 and 48 h post treatment. Prior to RNA isolation all treatments were examined under a light microscope for signs of visible cytotoxicity (i.e., cell death). Five independent cell-cultures (1 well of 6 well plate  5) per each time point/dose combination were performed to obtain 5 separate biological replicates. 2.4. Gene expression microarray measurements Total RNA was isolated at 24 and 48 h post treatment, using TriReagent (Molecular Research Center, Inc., Cincinnati, OH), further purified using a Qiagen RNeasy kit (Valencia, CA) and RNA integrity was validated using a NanoDrop 8000 Spectrophotometer (Wilmington, DE). Labeled cRNA was synthesized from 1 mg of total RNA from the 24 h and 48 h time point samples using the MessageAmpII Biotin-Enhanced labeling kit (Ambion, Austin, TX) and Agencourt RNAClean beads according to manufacturer’s instructions (Beckman Coulter, Inc., Danvers, MA). Twenty micrograms of labeled cRNA was fragmented and hybridized to Affymetrix Rat Genome 230 2.0 array (Santa Clara, CA) for 16 h at 45  C. Hybridized arrays were washed using a GeneChip Fluidics Station 450 and scanned using a GeneChip 3000 scanner (Affymetrix). 2.5. Gene expression microarray statistical analysis For gene expression analysis, scanned output files of Affymetrix microarrays were visually inspected for hybridization artifacts and then analyzed using the Affymetrix Expression Console Software and MAS5 (Affymetrix Microarray Analysis Suite 5) Statistical algorithm (http://www.affymetrix.com/products). Arrays were scaled to an average intensity of 1500 units and analyzed independently. The Affymetrix Rat Genome 230 2.0 array (Santa Clara, CA) used in this study has more than 31,000 probe sets analyzing over 30,000 transcripts and variants from over 28,000 well substantiated rat genes. The complete gene expression data have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) (accession no.: GSE40348) Sample mRNA quality was assessed using both the overall quality measures provided by the Affymetrix software and by

Table 2 Chemicals and corresponding vehicle, dose and number of significantly altered genes for the four different time dose combinations. No. of significant genes T1  D2 (0.05 < FDR)

Chemical

Vehicle Conc. (low)

Conc. (high)

Acetaminophen

1% DMSO 1% DMSO 1% DMSO 1% DMSO 1% DMSO 1% DMSO 1% DMSO H2O H2O 1% DMSO

0.5 mM

5 mM

209

4888

4914

3673

1 mM

100 mM

648

3560

163

2871

0.8 mM

20 mM

14

3297

1018

6597

60 mM

1000 mM

627

4645

6733

5062

250 mM

1000 mM

2644

2836

228

5424

250 mM

1000 mM

2589

2635

3505

3293

482

5258

2

5146

5792 353 4199

5459 11,442 4301

5555 243 5304

10,066 6916 5339

b-Naphthoflavone Chlorpromazine Clofibrate Bis(2-ethylhexyl) phthalate (DEHP) Diisononyl phthalate (DINP) Methapyrilene Phenobarbital Sodium valproate WY14643

3250 mM 100 mM 300 mM 0.5 mM 8 mM

3000 mM 10 mM 200 mM

No. of significant genes T1  D1 (0.05 < FDR)

DMSO: dimethyl sulfoxide. Grayed column depicts time dose combination used for further analysis.

No. of significant genes T2  D1 (0.05 < FDR)

No. of significant genes T2  D2 (0.05 < FDR)

32

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

statistical outlier detection methods. No samples were excluded for quality-related issues. Initially, our goal was to identify genes exhibiting differential expression between any of the 10 different hepatotoxicants and the control at each dose–timepoint condition – low (D1) and high (D2) doses at 24 h (T1) and 48 h (T2) timepoints, resulting in 4 different comparisons: T1  D1, T1  D2, T2  D1 and T2  D2. The number of statistically-significant gene expression changes across these dose–timepoint comparisons would then be used to identify the most suitable dose–timepoint combination for predictive toxicology. The differential expression analysis was accomplished by fitting one-way analysis of variance (ANOVA) models and calculating empirical Bayes regularized F statistics and corresponding p- and q-values for each probeset (Smyth, 2004). Separate models were fit for each dose–timepoint combination. This analysis was conducted using the limma package available from the R Bioconductor repository (http://www.bioconductor. org/). Hierarchical clustering was then performed on statistically significant probesets (defined as q < 0.05) to visually assess hepatotoxicant groupings.

treatments. Gene Ontology annotations provided by Affymetrix was used to perform a geneset enrichment analysis using a mean rank geneset enrichment test (MR-GSE) – also available through R Bioconductor’s limma package (http://www.bioconductor.org/). Genes were first ranked according to the above calculated F-statistic. The highest-ranked genes are those genes exhibiting the largest expression changes between (any of the) chemical treatments and the controls. Then, for each geneset (pathway), the ranks of the genes in D1  T1 were compared to the ranks of the remaining genes using a Wilcoxon rank sum test (Michaud et al., 2008). In order to visualize the pathway enrichment, a hierarchical clustering was then performed using the genes belonging to the enriched genesets: genesets with statistically significant MR-GSE test p-values using a family-wise error rate of 5% (Bonferroni adjusted p-value < 0.05). As a final analysis to identify both biologically relevant and statistically significant genes, a hierarchical clustering was performed using the overlap of the above two assessments. Probesets that were statistically significant themselves and belonged to a statistically significant geneset were clustered.

2.6. Geneset (pathway) enrichment

2.7. Analysis of TG GATES data

Following the probeset-level analyses, pathway analysis was conducted to identify biological themes that could explain the differential expression changes observed across the toxicological

In order to validate our method and understand if the process described in this manuscript was able to cluster in vivo hepatotoxicity data in a manner similar to what was seen in the

Fig. 1. Flow chart of step wise work process for gene expression microarray data analysis.

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

33

in vitro result, we analyzed publicly available in vivo data using an approach similar to what was described in this paper. CEL files were obtained for the chemicals PB, CPX, MP, APAP, CLO and WY14643 from TG-GATEs (Toxicogenomics Project-Genomics Assisted Toxicity Evaluation system) database, a publicly available database covering 170 compounds (Hirode et al., 2009; Kiyosawa et al., 2006, 2009; Uehara et al., 2011). The phthalates and BNF used in the in vitro experiments were not available on the TG-GATES website. VPA, while available on the TG GATES website was not used for further analysis since the number of significantlyexpressed genes was substantially lower for this treatment as compared to the other 6 chemicals. CEL files for the PB, CPX, MP, APAP, CLO and WY14643, 24 h timepoint, high dose were preprocessed using robust multi-array average (RMA) and log2 transformed (Irizarry et al., 2003). A oneway ANOVA was performed on the expression data for chemical treated and control samples separately and the procedure repeated for each chemical (Partek Genomic Suite, Version 6.6; Partek, Inc., St. Louis, MO). The p-values for all contrasts were adjusted for multiple comparisons using false discovery rate (FDR; Hochberg and Benjamini, 1990). All genes with an FDR corrected pvalue < 0.05 (FDR < 0.05) were used as the criteria for determining significance and were used for further analysis Functional ontology enrichment analysis was performed on the genes that exhibited significantly altered expression for each chemical. The enrichment analysis was conducted using the GeneGo process networks in the Metacore database (Thomson Reuters, Philadelphia, PA). The enrichment p-values were calculated based on a hypergeometric distribution and significant enrichment was defined as an FDR-corrected p-value < 0.05. A hierarchical clustering was performed for the probesets identified via the above process (Fig. 5). 3. Results 3.1. Gene expression changes in primary rat hepatocytes following exposure to ten different hepatotoxicants Gene expression microarray analysis was performed on primary rat hepatocytes exposed to 10 different hepatotoxicants at two different doses for 24 h and 48 h. The number of significant genes for each of the time/dose combinations ranged from 2 to 11,442 (Table 2). The time/dose comparison T1  D2 (24 h  high dose) was used for further analysis since this time/dose combination was considered most suitable for predictive toxicology assessments. The high dose was selected based on the higher number of differentially expressed genes and the earlier time point of 24 h was selected to minimize gene/pathway enrichment due to secondary transcription events which are likely to occur by 48 h post treatment (Zhang et al., 2014).

Fig. 2. Whole genome hierarchical clustering of differentially expressed (q < 0.05) probesets for the 24 h, high dose samples. Genes differentially expressed vs. DMSO control for any of the 10 hepatotoxicants were identified using an empirical Bayes regularized F-statistic.

3.2. Hierarchical clustering based comparison of significantly altered genes for ten hepatotoxicants Hierarchical clustering was performed on statistically significant probesets (defined as q < 0.05) for T1  D2 to visually assess hepatotoxicant groupings. The results showed a roughly similar distribution in the number of up and down regulated genes for each of the hepatotoxicants even though the groups of genes that were up and down regulated were quite varied between the treatments (Fig. 2). The dendrogram exhibited a clear clustering of 2 groups (of 4 chemicals) and 2 chemicals that did not cluster with either group or with each other. The grouping to the far right (Fig. 2) included the two structurally similar compounds, DEHP and DINP along with CLO and WY14643 (group 1 chemicals) (Fig. 2). The 2nd clustering included

PB, MP, BNF and CPZ (group 2 chemicals). VPA and APAP showed minimum similarity to either of the groups or to each other. 3.3. Pathway based enrichment analysis of significantly altered genes for ten hepatotoxicants In order to understand how gene expression following chemical treatment affected cellular pathways, Gene Ontology annotations provided by Affymetrix were used to perform a geneset enrichment analysis (Supp. Table 1). The results showed that general liver function/toxicity themes were enriched when GO terms were filtered to represent only the biological pathway (BP) terms (Table 3). Following this analysis a hierarchical

34

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

Table 3 Top 25 BP GO themes sorted by p-value of parent GO term. go_id

Label

Ontology

nProbes

p-value T1  D2

GO:0042493 GO:0001889 GO:0014070 GO:0071320 GO:0071549 GO:0071392 GO:0071385 GO:0071407 GO:0031000 GO:0051591 GO:0042220 GO:0031960 GO:0051412 GO:0032355 GO:0043627 GO:0051593 GO:0051384 GO:0043278 GO:0035094 GO:0032570 GO:0048545 GO:0033574 GO:0033280 GO:0033197 GO:0006635 GO:0007584 GO:0051593 GO:0033591 GO:0033189 GO:0033280 GO:0033197 GO:0006631 GO:0019369 GO:0006635 GO:0006633 GO:0019395 GO:0001516 GO:0010033 GO:0071230 GO:0071320 GO:0071549 GO:0071392 GO:0071398 GO:0071385 GO:0071333 GO:0071363 GO:0032870 GO:0032869 GO:0071347 GO:0071222 GO:0071407 GO:0071375 GO:0071356 GO:0019221 GO:0008543 GO:0008286 GO:0046627 GO:0014075 GO:0043200 GO:0001975 GO:0031000 GO:0051591 GO:0009743 GO:0042220 GO:0031960 GO:0051412 GO:0034097 GO:0070849 GO:0032355 GO:0043627 GO:0045471 GO:0070542 GO:0051593 GO:0033762 GO:0051384 GO:0009749

Response to drug Liver development Response to organic cyclic substance Cellular response to cAMP Cellular response to dexamethasone stimulus Cellular response to estradiol stimulus Cellular response to glucocorticoid stimulus Cellular response to organic cyclic substance Response to caffeine Response to cAMP Response to cocaine Response to corticosteroid stimulus Response to corticosterone stimulus Response to estradiol stimulus Response to estrogen stimulus Response to folic acid Response to glucocorticoid stimulus Response to morphine Response to nicotine Response to progesterone stimulus Response to steroid hormone stimulus Response to testosterone stimulus Response to vitamin D Response to vitamin E Fatty acid beta-oxidation Response to nutrient Response to folic acid Response to L-ascorbic acid Response to vitamin A Response to vitamin D Response to vitamin E Fatty acid metabolic process Arachidonic acid metabolic process Fatty acid beta-oxidation Fatty acid biosynthetic process Fatty acid oxidation Prostaglandin biosynthetic process Response to organic substance Cellular response to amino acid stimulus Cellular response to cAMP Cellular response to dexamethasone stimulus Cellular response to estradiol stimulus Cellular response to fatty acid Cellular response to glucocorticoid stimulus Cellular response to glucose stimulus Cellular response to growth factor stimulus Cellular response to hormone stimulus Cellular response to insulin stimulus Cellular response to interleukin-1 Cellular response to lipopolysaccharide Cellular response to organic cyclic substance Cellular response to peptide hormone stimulus Cellular response to tumor necrosis factor Cytokine-mediated signaling pathway Fibroblast growth factor receptor signaling pathway Insulin receptor signaling pathway Negative regulation of insulin receptor signaling pathway Response to amine stimulus Response to amino acid stimulus Response to amphetamine Response to caffeine Response to cAMP Response to carbohydrate stimulus Response to cocaine Response to corticosteroid stimulus Response to corticosterone stimulus Response to cytokine stimulus Response to epidermal growth factor stimulus Response to estradiol stimulus Response to estrogen stimulus Response to ethanol Response to fatty acid Response to folic acid Response to glucagon stimulus Response to glucocorticoid stimulus Response to glucose stimulus

BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP

519 82 292 11 18 10 27 24 23 92 35 15 30 161 101 16 149 29 58 44 57 43 29 15 22 157 16 11 32 29 15 31 17 22 12 11 10 152 10 11 18 10 14 27 17 33 40 79 21 19 24 20 18 22 16 14 12 21 50 15 23 92 26 35 15 30 95 10 161 101 147 30 16 15 149 96

3.87E-10 5.96E-14 3.32E-08 0.080028 0.02321 0.996711 0.445833 0.017672 0.506574 0.000102 0.234819 0.866358 0.058602 0.408136 0.047749 0.039691 1.23E-06 0.648256 0.318268 0.235842 4.31E-05 0.141721 0.072151 0.003144 1.33E-08 3.21E-09 0.039691 0.134535 0.506413 0.072151 0.003144 1.47E-07 0.009079 1.33E-08 0.007201 0.027963 0.976887 3.35E-08 0.744586 0.080028 0.02321 0.996711 0.300306 0.445833 0.11645 0.773396 0.026249 7.31E-07 0.408391 0.045991 0.017672 0.017062 0.307099 0.598847 0.904173 0.007209 0.132433 0.10427 5.79E-07 0.674827 0.506574 0.000102 0.040422 0.234819 0.866358 0.058602 0.259761 0.252361 0.408136 0.047749 0.000599 0.058793 0.039691 0.008582 1.23E-06 0.001637

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

35

Table 3 (Continued) go_id

Label

Ontology

nProbes

p-value T1  D2

GO:0034698 GO:0070848 GO:0060416 GO:0009725 GO:0032868 GO:0034341 GO:0070555 GO:0033591 GO:0033993 GO:0032496 GO:0051597 GO:0043278 GO:0035094 GO:0014070 GO:0010243 GO:0043434 GO:0032570 GO:0032526 GO:0048545 GO:0009744 GO:0033574 GO:0033189 GO:0033280 GO:0033197 GO:0007179 GO:0006979 GO:0070301 GO:0042542 GO:0000302 GO:0043200 GO:0071230 GO:0051593 GO:0046686 GO:0051384 GO:0071549 GO:0071385 GO:0051412 GO:0006695

Response to gonadotropin stimulus Response to growth factor stimulus Response to growth hormone stimulus Response to hormone stimulus Response to insulin stimulus Response to interferon-gamma Response to interleukin-1 Response to L-ascorbic acid Response to lipid Response to lipopolysaccharide Response to methylmercury Response to morphine Response to nicotine Response to organic cyclic substance Response to organic nitrogen Response to peptide hormone stimulus Response to progesterone stimulus Response to retinoic acid Response to steroid hormone stimulus Response to sucrose stimulus Response to testosterone stimulus Response to vitamin A Response to vitamin D Response to vitamin E Transforming growth factor beta receptor signaling pathway Response to oxidative stress Cellular response to hydrogen peroxide Response to hydrogen peroxide Response to reactive oxygen species Response to amino acid stimulus Cellular response to amino acid stimulus Response to folic acid Response to cadmium ion Response to glucocorticoid stimulus Cellular response to dexamethasone stimulus Cellular response to glucocorticoid stimulus Response to corticosterone stimulus Cholesterol biosynthetic process

BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP BP

10 29 11 119 67 10 26 11 27 167 17 29 58 292 85 171 44 58 57 10 43 32 29 15 24 86 24 76 17 50 10 16 33 149 18 27 30 14

0.080146 0.04162 2.31E-05 2.02E-05 0.021409 0.55787 0.90112 0.134535 0.060443 0.095619 0.071982 0.648256 0.318268 3.32E-08 0.000211 0.015989 0.235842 0.044776 4.31E-05 0.092247 0.141721 0.506413 0.072151 0.003144 0.619929 6.64E-07 0.098892 0.003501 0.758269 5.79E-07 0.744586 0.039691 5.72E-09 1.23E-06 0.02321 0.445833 0.058602 3.35E-06

All child BP GO terms are indented below parent BP GO term. Each BP GO theme is contained within two solid lines.

clustering was performed using genes belonging to statistically significant gene sets (MR-GSE p < 0.05). While the gene set dendogram (Fig. 3) exhibited a very similar clustering to the probe set clustering (Fig. 2), the number of genes used for the geneset dendogram was significantly lower (4212 probes vs. 14,468 probes). This suggests that a few specific pathways may be responsible for driving the MOA of a chemical. Furthermore, traditional hierarchical clustering methods which use all significantly altered genes for hierarchical clustering may potentially be introducing more noise to the analysis by using genes that are not enriched via pathways. Supplementary material related to this article found, in the online version, at http://dx.doi.org/10.1016/j.tox.2014.11.008. 3.4. Hierarchical clustering of genes common to probe set enrichment and gene set enrichment Pathway enrichment provided a biological reasoning for the underlying toxicology that is taking place (pathways that are driving the MOA). The probesets provided the actual genes that are being regulated in order to activate these particular pathways. However we noted that for both of these analyses (Figs. 2 and 3) a significant number of genes whilst showing significantly altered gene expression for a couple of chemicals did not show much change for the rest of the chemicals (white areas of dendogram). Hence to get at a more robust cluster analysis, that focused on MOA that overlaps between chemicals, a final clustering was performed using only genes that were statistically significant themselves and also belonged to a statistically significant geneset (overlap of

significant probesets and genesets: Figs. 1 and 4 ). This analysis reduced the number of genes to 2478 and preserved the clustering pattern seen in Figs. 2 and 3. Obvious patterns in gene expression were visible in this dendogram (Fig. 4 compared to Figs. 2 and 3) with groups of genes with similar up/down regulation grouping together providing a more tighter clustering. We believe by clustering the overlap of steps 1 and 2 (Fig. 1) a refined number of genes were used for clustering (Fig. 4). There by limiting the analysis to a smaller critical set of genes which lead to a more weighted clustering (reduced noise). 3.5. Hierarchical clustering of genes common to probe set enrichment and gene set enrichment using publicly available (TG-GATES) in vivo data In order to truly understand the applicability of our method, we wanted to test it merits using a well-controlled externally available in vivo data set. We were able to obtain a robust data set for 6 of the 10 chemicals on the TG-GATEs database (Hirode et al., 2009; Kiyosawa et al., 2006, 2009; Uehara et al., 2008). We performed an analysis similar to what was performed for the in vitro data set and completed a final clustering for genes that were statistically significant themselves and also belonged to a statistically significant geneset (overlap of significant probesets and genesets) for the TG-GATES data (Fig. 5). Then compared the assessment of the in vivo data set to the in vitro result (Fig. 5 compared to Fig. 4). The two assessments showed good correlation for the chemicals that have receptor-mediated effects with weaker correlation for chemicals for which the MOA appears to involve covalent reactivity.

36

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

Fig. 3. Hierarchical clustering of probesets belonging to enriched genesets for the 24 h, high dose samples. Genes were ranked according to their overall differential expression vs. DMSO control using an empirical Bayes regularized F-statistic. Geneset enrichment for Gene Ontology terms was evaluated using the mean-rank geneset enrichment test. Clustering includes only probesets belonging to one or more genesets with Bonferroni-adjusted p < 0.05.

Fig. 4. Hierarchical clustering of differentially expressed probesets belonging to one or more enriched genesets for the 24 h, high dose samples. Genes differentially expressed vs. DMSO control for any of the 10 hepatotoxicants were identified using an empirical Bayes regularized F-statistic. Enriched genesets for the overall comparison vs. DMSO were identified using a mean-rank geneset enrichment test. Clustering includes only probesets that were both differentially expressed (q < 0.05) and belonged to an enriched geneset (Bonferroni-adjusted p < 0.05).

4. Discussion

clustering of gene lists from individual pathways (unsupervised clustering based on pathways). In combination with the traditional hierarchical clustering method the novel geneset/pathway based method provides increased weight of evidence for discriminating and predicting chemical toxicity, based on mode of action. This method shows that gene expression experiments using primary hepatocytes are much more data rich than previously anticipated. Since the 10 chemicals used in this study were well studied compounds, a majority of them had in vitro toxicity data. Hence, we choose to extrapolate from this data to determine the most suitable in vitro doses for our experiments (Table 1). Our results

Multiple studies have attempted to classify toxicants based on their gene expression profiles. Most such studies have been limited to animal models, run on gene expression microarrays with a limited number of probe sets, used hand-picked sets of genes for analysis (supervised clustering) or used statistical methods that are highly restrictive and conservative. In the current study we address these pitfalls and present a novel more biologically relevant method that uses primary rat hepatocytes, whole genome microarrays and an analysis method that is based on hierarchical

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

37

Fig. 5. Hierarchical clustering of differentially expressed probesets belonging to one or more enriched genesets for the 24 h, high dose samples for the TG-GATEs data set. Genes differentially expressed vs. DMSO control for any of the 6 hepatotoxicants were identified using p-value adjusted for multiple comparisons (FDR < 0.05). Functional ontology enrichment analysis was performed on the significantly altered expression for each chemical identified from the above step using the GeneGo process networks in the Metacore database (Thomson Reuters, Philadelphia, PA). The enrichment p-values were calculated based on a hypergeometric distribution and significant enrichment was defined as an FDR-corrected p-value < 0.05. Clustering was performed on probesets that were both differentially expressed (FDR < 0.05) and belonged to an enriched geneset (FDR < 0.05).

showed that the higher exposure dose resulted in the higher number of altered genes except for PB, where a slightly higher number of significantly altered genes were observed for T1  D1 vs. T1  D2 (Table 2). Overall, these total significantly altered gene number observations suggests the need to use a sufficiently high dose of the chemical when performing in vitro experiments using primary cells/cell lines in order to obtain a sufficient number of significantly altered genes to perform a robust pathway analysis. Whilst our study doses were based on calculations from in vitro data on the chemicals of interest (Table 1), this method would not be an option for novel chemicals with no in vitro data. In such cases, following an approach such as the one described in the OECD “guidance document on using cytotoxicity tests to estimate starting doses for acute oral systemic toxicity tests” (OECD, 2010) could be used. A dose equivalent to an EC50 could be calculated, this dose could be used as the highest dose and a number of serial dilutions of this dose could be adopted such that the potential human exposure concentration is flanked by the concentrations within the serial dilution scheme. In general, all 3 of the hierarchical clustering methods (probeset based hierarchical clustering, gene set based hierarchical clustering and overlap of probeset and gene set hierarchical clustering) (Fig. 1) managed to group the 10 chemicals in to the anticipated toxicological classes for T1  D2. However, clustering was not consistent and did not follow anticipated groupings based on chemical toxicity class for the other dose time combinations for all 3 of the hierarchical clustering methods (Supp. Figs. 1–3), providing further evidence for the sensitivity of transcriptomic experiments to dose and time point. Supplementary material related to this article found, in the online version, at http://dx.doi.org/10.1016/j.tox.2014.11.008. The first of these groups of chemicals, CLO, WY14643, DEHP, and DINP consisted of peroxisome proliferators (Corton and Lapinskas, 2005; Lapinskas et al., 2005; Tamura et al., 2006), the second group, PB, MP, BNF and CPZ consisted of cytochrome

P450 inducers (Graham et al., 2006; Mingoia et al., 2007; Uehara et al., 2008) and the last two chemicals that did not cluster with either group consisted of chemicals that caused cytotoxicity predominantly through covalent binding to macro molecules (Cerveny et al., 2007; Rogiers et al., 1995; Tujios and Fontana, 2011). Of the 2 groups and 2 individual chemicals, the MOA of the individual chemicals (APAP and VPA) was the least specific. A survey of the current literature shows that the exact mechanism by which these two compounds cause liver injury are not clearly defined (Cerveny et al., 2007; Rogiers et al., 1995; Tujios and Fontana, 2011). Furthermore, the origin of toxicity for each of these compounds cannot be associated with modification of a single pathway. Instead, both of these chemicals seem to impart its toxicity through multiple mechanisms. VPA has been identified as a histone deacytalase inhibitor, an inducer of certain P450s and been suggested to have a role in the activation of the constitutive androstane receptor (CAR) and/or pregnane X receptor (PXR) pathways (Cerveny et al., 2007; Rogiers et al., 1995). APAP’s toxicity is mediated through a highly reactive metabolite, N-acetyl-pbenzoquinone imine (NAPQI). Mechanisms of APAP’s toxicity known to date include interference with mitochondrial and nuclear function, generation of reactive oxygen species along with mounting evidence for a NAPQI mediated host immune response role (Tujios and Fontana, 2011). Furthermore, the reactive chemistry behind these two compounds has the potential to activate multiple pathways, hence the mechanisms are non-specific and the pathways that are enriched are not a true representation of the chemical’s unique MOA. The activation of these various nonspecific pathways and the minimal overlap of these pathways with groups 1 and 2 chemicals may explain the clustering of these chemicals to clusters of 1 (clustered independently and not with either group) (Figs. 2–4). To further gauge the potential of our method to cluster gene expression based on MOA we tested our method on a well-controlled externally available in vivo data set. TG-GATEs is

38

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39

a publicly available database covering 170 compounds (Hirode et al., 2009; Kiyosawa et al., 2006, 2009; Uehara et al., 2011). The data include time-series gene expression for in vivo experiments in Rattus norvegicus using the Affymetrix Rat Genome 230 2.0 array, which is the same array that was used in our experiments We clustered genes that were statistically significant themselves and also belonged to a statistically significant geneset for the 6 chemicals available on the TG-GATEs database (Fig. 5). Then compared the TG-GATEs result to the in vitro result (Fig. 4 compared to Fig. 5). We observed that the two, group 1 chemicals (peroxisome proliferators) in the TG-GATEs data set, WY14643 and CLO clustered together and showed good correlation to the in vitro result. However, the correlation between the in vitro and in vivo result for the three group 2 chemicals (cytochrome P450 inducers), CPZ, PB and MP were weak. CPZ and PB showed more similarity to each other than to any of the other chemicals, even though the clustering was not as strong as the in vitro result. MP showed the most deviation between the in vitro and in vivo result and clustered with APAP instead of CPZ and PB when using the TG-GATEs data set. The in vivo in vitro discrepancy of MP is well documented and has been attributed to pharmacokinetic differences between the in vitro and in vivo systems. (Schug et al., 2013). Furthermore, downstream effects of MP results in cytotoxicity (Omura et al., 2014). If genes representing pathways that mediate cytotoxicity (cytotoxicity related pathways) were dominating the expression profile of MP, this may also explain MP clustering with APAP, another chemical known to induce significant hepatic cytotoxicity. While CPZ and PB did not cluster as tightly as the in vitro experiment when using TG-GATEs data, they were the most unlike the other 4 chemicals. The significant separation of CPZ and PB from MP could be attributed to the details of their MOAs. While all three of these chemicals can be very generally categorized as CYP inducers the process by which they promote the expression of specific CYPs are quite different. PB affects are mediated via CAR (Elcombe et al., 2014), CPZ via farnesoid X receptor (FXR) (Szalowska et al., 2013) and MP via the induction of oxidative stress through multiple mechanisms that lead to DNA damage and cytotoxicity (Omura et al., 2014). The TG-GATEs result showed a pronounced clustering of the P450 inducers, that mediate their MOA via nuclear receptors (CAR, FXR) as compared to the in vitro result. Other factors contributing to the differences between the in vitro and TG-GATEs result, include the number of chemicals that were being assessed in each procedure. The lesser the number of chemicals that are being assessed the more the chemicals maybe “forced” to cluster with other chemicals in the data set, hence the number of chemicals used in the clustering (10 vs. 6) may also contribute to the differences observed. Overall, this comparison showed that, dose and time point had significant impact on in vivo data similar to its impact on in vitro data. In addition, other factors such as pharmacokinetics and number of chemicals used for an experiment (size of database) are important considerations for a method such as the above. Considering that we used an externally available in vivo data set with no control over dose and timepoint, we think that the correlation between the two data sets is substantial, however careful attention should be given to the factors listed about in order to obtain a meaningful result from a method such as ours. In order to further understand what factors may be contributing to discrepancies between the TG-GATEs data and our in vitro data we performed a deeper dive on two representative chemicals and compared the results from our experiment to the TG-GATES result. A single chemical was picked from group 1: WY14643 and group 2: MP chemicals and the respective data sets were then compared to the respective in vivo high dose data set for the 24 h time point on the TG-GATES database. For the group 1 representative, WY14643, the number of significantly altered genes were 6440 and 2727 for

our (in vitro) and TG-GATEs data set respectively. Of these, 1475 (54%) overlapped, showing good correlation. To further assess the biological significance of these genes, we performed gene ontology enrichment analysis for each data set separately and compared the top twenty most enriched pathways from our (in vitro) data to the TG-GATEs results. The comparison showed that 13 of the top 20 pathways from our results falling within the top 50 of the TG-GATEs data (Supp. Fig. 4). A similar analysis was done for the group 2 representative, MP. For this analysis the number of significantly altered genes were 7881 and 1462 for our (in vitro) and TG-GATEs data, respectively. Of these, 956 (56%) showed overlap. Pathway enrichment analysis showed that 10 pathways from the top twenty most enriched pathways from our (in vitro) data fell within the top 50 of the TG-GATEs data (Supp. Fig. 5). It is widely accepted that in vitro experiments result in a higher number of significantly altered genes when compared to performing the same experiment in vivo. This same phenomenon was observed in our analysis (Supp. Figs. 4 and 5). When this discrepancy is taken into account, we believe that the correlation we saw between the in vivo and in vitro experiments were quite significant. Supplementry material related to this article found, in the online version, at http://dx.doi.org/10.1016/j.tox.2014.11.008. Several interesting conclusions can be drawn from this study. First, we provide further evidence for the already established view that transcriptomic profiling can be used to classify and thereby predict the toxicological profile of a chemical. Most importantly, this study offers a more complete MOA based method to perform such analyses using non-supervised clustering, in an in vitro model of primary hepatocytes. In addition, both the analysis of our in vitro data set and an externally available in vivo data set showed that factors such as dose, time point and size of database should be well controlled to obtain meaningful results from a method such as this. With increased pressure to reduce and eliminate animal use, this type of approach may provide an alternative, tailored testing method that is underpinned by mechanistic understanding. Conflict of interest The authors declare that there are no conflicts of interest. Transparency document The Transparency document associated with this article can be found in the online version. Acknowledgements The authors thank Karen Blackburn, Robert Binder and Raghunandan Kainkaryam for their review of the manuscript. References Bulera, S.J., Eddy, S.M., Ferguson, E., Jatkoe, T.A., Reindel, J.F., Bleavins, M.R., De La Iglesia, F.A., 2001. RNA expression in the early characterization of hepatotoxicants in Wistar rats by high-density DNA microarrays. Hepatology 33, 1239–1258. Cerveny, L., Svecova, L., Anzenbacherova, E., Vrzal, R., Staud, F., Dvorak, Z., Ulrichova, J., Anzenbacher, P., Pavek, P., 2007. Valproic acid induces CYP3A4 and MDR1 gene expression by activation of constitutive androstane receptor and pregnane X receptor pathways. Drug Metab. Dispos.: Biol. Fate Chem. 35, 1032–1041. Chen, M., Zhang, M., Borlak, J., Tong, W., 2012. A decade of toxicogenomic research and its contribution to toxicological science. Toxicol. Sci.: Off. J. Soc. Toxicol. 130, 217–228. Choi, S., et al., 2010. Identification of toxicological biomarkers of di(2-ethylhexyl) phthalate in proteins secreted by HepG2 cells using proteomic analysis. Proteomics 10 (9), 1831–1846. Corton, J.C., Lapinskas, P.J., 2005. Peroxisome proliferator-activated receptors: mediators of phthalate ester-induced effects in the male reproductive tract? Toxicol. Sci.: Off. J. Soc. Toxicol. 83, 4–17.

K.N. De Abrew et al. / Toxicology 328 (2015) 29–39 Daston, G.P., Naciff, J.M., 2010. Predicting developmental toxicity through toxicogenomics. Birth Defects Res. Part C, Embryo: Rev. 90, 110–117. Elcombe, C.R., Peffer, R.C., Wolf, D.C., Bailey, J., Bars, R., Bell, D., Cattley, R.C., Ferguson, S.S., Geter, D., Goetz, A., Goodman, J.I., Hester, S., Jacobs, A., Omiecinski, C.J., Schoeny, R., Xie, W., Lake, B.G., 2014. Mode of action and human relevance analysis for nuclear receptor-mediated liver toxicity: a case study with phenobarbital as a model constitutive androstane receptor (CAR) activator. Crit. Rev. Toxicol. 44, 64–82. Gerets, H.H., Tilmant, K., Gerin, B., Chanteux, H., Depelchin, B.O., Dhalluin, S., Atienzar, F.A., 2012. Characterization of primary human hepatocytes, HepG2 cells, and HepaRG cells at the mRNA level and CYP activity in response to inducers and their predictivity for the detection of human hepatotoxins. Cell Biol. Toxicol. 28, 69–87. Goetz, A.K., Singh, B.P., Battalora, M., Breier, J.M., Bailey, J.P., Chukwudebe, A.C., Janus, E.R., 2011. Current and future use of genomics data in toxicology: opportunities and challenges for regulatory applications. Regul. Toxicol. Pharmacol. 61 (2), 141–153. Goyak, K.M., et al., 2008. Expression profiling of interindividual variability following xenobiotic exposures in primary human hepatocyte cultures. Toxicol. Appl. Pharmacol. 231 (2), 216–224. Graham, R.A., Tyler, L.O., Krol, W.L., Silver, I.S., Webster, L.O., Clark, P., Chen, L., Banks, T., LeCluyse, E.L., 2006. Temporal kinetics and concentration–response relationships for induction of CYP1A, CYP2B, and CYP3A in primary cultures of beagle dog hepatocytes. J. Biochem. Mol. Toxicol. 20, 69–78. Hamadeh, H.K., Amin, R.P., Paules, R.S., Afshari, C.A., 2002a. An overview of toxicogenomics. Curr. Issues Mol. Biol. 4, 45–56. Hamadeh, H.K., Bushel, P.R., Jayadev, S., DiSorbo, O., Bennett, L., Li, L., Tennant, R., Stoll, R., Barrett, J.C., Paules, R.S., Blanchard, K., Afshari, C.A., 2002b. Prediction of compound signature using high density gene expression profiling. Toxicol. Sci.: Off. J. Soc. Toxicol. 67, 232–240. 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., 2002c. Gene expression analysis reveals chemical-specific profiles. Toxicol. Sci.: Off. J. Soc. Toxicol. 67, 219–231. Hasmall, S.C., et al., 1999. Suppression of apoptosis and induction of DNA synthesis in vitro by the phthalate plasticizers monoethylhexylphthalate (MEHP) and diisononylphthalate (DINP): a comparison of rat and human hepatocytes in vitro. Arch. Toxicol. 73 (8–9), 451–456. Hirode, M., Omura, K., Kiyosawa, N., Uehara, T., Shimuzu, T., Ono, A., Miyagishima, T., Nagao, T., Ohno, Y., Urushidani, T., 2009. Gene expression profiling in rat liver treated with various hepatotoxic-compounds inducing coagulopathy. J. Toxicol. Sci. 34, 281–293. Hochberg, Y., Benjamini, Y., 1990. More powerful procedures for multiple significance testing. Stat. Med. 9, 811–818. Irizarry, R.A., Bolstad, B.M., Collin, F., Cope, L.M., Hobbs, B., Speed, T.P., 2003. Summaries of Affymetrix GeneChip probe level data. Nucl. Acids Res. 31, e15. Jetten, M.J., Kleinjans, J.C., Claessen, S.M., Chesne, C., van Delft, J.H., 2013. Baseline and genotoxic compound induced gene expression profiles in HepG2 and HepaRG compared to primary human hepatocytes. Toxicol. In Vitro: Int. J. Publ. Assoc. BIBRA 27, 2031–2040. Jigorel, E., et al., 2006. Differential regulation of sinusoidal and canalicular hepatic drug transporter expression by xenobiotics activating drug-sensing receptors in primary human hepatocytes. Drug Metab. Dispos. 34 (10), 1756–1763. Kiyosawa, N., Shiwaku, K., Hirode, M., Omura, K., Uehara, T., Shimizu, T., Mizukawa, Y., Miyagishima, T., Ono, A., Nagao, T., Urushidani, T., 2006. Utilization of a one-dimensional score for surveying chemical-induced changes in expression levels of multiple biomarker gene sets using a large-scale toxicogenomics database. J. Toxicol. Sci. 31, 433–448. Kiyosawa, N., Ando, Y., Watanabe, K., Niino, N., Manabe, S., Yamoto, T., 2009. Scoring multiple toxicological endpoints using a toxicogenomic database. Toxicol. Lett. 188, 91–97. Lapinskas, P.J., Brown, S., Leesnitzer, L.M., Blanchard, S., Swanson, C., Cattley, R.C., Corton, J.C., 2005. Role of PPARalpha in mediating the effects of phthalates and metabolites in the liver. Toxicology 207, 149–163. Michaud, J., Simpson, K.M., Escher, R., Buchet-Poyau, K., Beissbarth, T., Carmichael, C., Ritchie, M.E., Schutz, F., Cannon, P., Liu, M., Shen, X., Ito, Y., Raskind, W.H., Horwitz, M.S., Osato, M., Turner, D.R., Speed, T.P., Kavallaris, M., Smyth, G.K.,

39

Scott, H.S., 2008. Integrative analysis of RUNX1 downstream pathways and target genes. BMC Genomics 9, 363. Mingoia, R.T., Nabb, D.L., Yang, C.H., Han, X., 2007. Primary culture of rat hepatocytes in 96-well plates: effects of extracellular matrix configuration on cytochrome P450 enzyme activity and inducibility and its application in in vitro cytotoxicity screening. Toxicol. In Vitro: Int. J. Publ. Assoc. BIBRA 21, 165–173. NRC, 2007a. Applications of Toxicogenomic Technologies to Predictive Toxicology and Risk Assessment. National Academies Press, Washington, D.C. NRC, 2007b. Toxicity Testing in the 21st Century: A Vision and a Strategy. National Academies Press, Washington, D.C. Nuwaysir, E.F., Bittner, M., Trent, J., Barrett, J.C., Afshari, C.A., 1999. Microarrays and toxicology: the advent of toxicogenomics. Mol. Carcinog. 24, 153–159. OECD, 2010. OECD Series on Testing and Assessment—Guidance Document on Using Cytotoxicity Tests to Estimate Starting Doses for Acute Oral Systematic Toxicity Tests. Paris, France. OECD, 2012. Draft Template, and Guidance on Developing and Assessing the Completeness of Adverse Outcome Pathways (AOPs). Paris, France. Omura, K., Uehara, T., Morikawa, Y., Hayashi, H., Mitsumori, K., Minami, K., Kanki, M., Yamada, H., Ono, A., Urushidani, T., 2014. Detection of initiating potential of non-genotoxic carcinogens in a two-stage hepatocarcinogenesis study in rats. J. Toxicol. Sci. 39, 785–794. Rogiers, V., Akrawi, M., Vercruysse, A., Phillips, I.R., Shephard, E.A., 1995. Effects of the anticonvulsant valproate on the expression of components of the cytochrome-P-450-mediated monooxygenase system and glutathione S-transferases. Eur. J. Biochem./FEBS 231, 337–343. Schaap, M.M., Zwart, E.P., Wackers, P.F., Huijskens, I., van de Water, B., Breit, T.M., van Steeg, H., Jonker, M.J., Luijten, M., 2012. Dissecting modes of action of non-genotoxic carcinogens in primary mouse hepatocytes. Arch. Toxicol. 86, 1717–1727. Schug, M., Stober, R., Heise, T., Mielke, H., Gundert-Remy, U., Godoy, P., Reif, R., Blaszkewicz, M., Ellinger-Ziegelbauer, H., Ahr, H.J., Selinski, S., Gunther, G., Marchan, R., Sachinidis, A., Nussler, A., Oberemm, A., Hengstler, J.G., 2013. Pharmacokinetics explain in vivo/in vitro discrepancies of carcinogen-induced gene expression alterations in rat liver and cultivated hepatocytes. Arch. Toxicol. 87, 337–345. Smyth, G.K., 2004. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3 Article 3. Szalowska, E., Stoopen, G., Groot, M.J., Hendriksen, P.J., Peijnenburg, A.A., 2013. Treatment of mouse liver slices with cholestatic hepatotoxicants results in down-regulation of Fxr and its target genes. BMC Med. Genomics 6, 39. Tamura, K., Ono, A., Miyagishima, T., Nagao, T., Urushidani, T., 2006. Profiling of gene expression in rat liver and rat primary cultured hepatocytes treated with peroxisome proliferators. J. Toxicol. Sci. 31, 471–490. Thomas, R.S., Rank, D.R., Penn, S.G., Zastrow, G.M., Hayes, K.R., Pande, K., Glover, E., Silander, T., Craven, M.W., Reddy, J.K., Jovanovich, S.B., Bradfield, C.A., 2001. Identification of toxicologically predictive gene sets using cDNA microarrays. Mol. Pharmacol. 60, 1189–1194. Tujios, S., Fontana, R.J., 2011. Mechanisms of drug-induced liver injury: from bedside to bench. Nat. Rev. Gastroenterol. Hepatol. 8, 202–211. Uehara, T., Kiyosawa, N., Hirode, M., Omura, K., Shimizu, T., Ono, A., Mizukawa, Y., Miyagishima, T., Nagao, T., Urushidani, T., 2008. Gene expression profiling of methapyrilene-induced hepatotoxicity in rat. J. Toxicol. Sci. 33, 37–50. Uehara, T., Minowa, Y., Morikawa, Y., Kondo, C., Maruyama, T., Kato, I., Nakatsu, N., Igarashi, Y., Ono, A., Hayashi, H., Mitsumori, K., Yamada, H., Ohno, Y., Urushidani, T., 2011. Prediction model of potential hepatocarcinogenicity of rat hepatocarcinogens using a large-scale toxicogenomics database. Toxicol. Appl. Pharmacol. 255, 297–306. van Leeuwen, D.M., Pedersen, M., Knudsen, L.E., Bonassi, S., Fenech, M., Kleinjans, J. C., Jennen, D.G., 2011. Transcriptomic network analysis of micronuclei-related genes: a case study. Mutagenesis 26, 27–32. 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. Toxicol. Lett. 120, 359–368. Zhang, J.D., Berntenis, N., Roth, A., Ebeling, M., 2014. Data mining reveals a network of early-response genes as a consensus signature of drug-induced in vitro and in vivo toxicity. Pharmacogenomics J. 14, 208–216.