Scholarly article on topic 'Case study on the utility of hepatic global gene expression profiling in the risk assessment of the carcinogen furan'

Case study on the utility of hepatic global gene expression profiling in the risk assessment of the carcinogen furan Academic research paper on "Biological sciences"

CC BY
0
0
Share paper
Academic journal
Toxicology and Applied Pharmacology
OECD Field of science
Keywords
{Furan / "Gene expression" / Toxicogenomics / "Indirect-acting genotoxic carcinogen" / "Non-genotoxic carcinogen" / "Risk assessment"}

Abstract of research paper on Biological sciences, author of scientific article — Anna Francina Jackson, Andrew Williams, Leslie Recio, Michael D. Waters, Iain B. Lambert, et al.

Abstract Furan is a chemical hepatocarcinogen in mice and rats. Its previously postulated cancer mode of action (MOA) is chronic cytotoxicity followed by sustained regenerative proliferation; however, its molecular basis is unknown. To this end, we conducted toxicogenomic analysis of B3C6F1 mouse livers following three week exposures to non-carcinogenic (0, 1, 2mg/kgbw) or carcinogenic (4 and 8mg/kgbw) doses of furan. We saw enrichment for pathways responsible for cytotoxicity: stress-activated protein kinase (SAPK) and death receptor (DR5 and TNF-alpha) signaling, and proliferation: extracellular signal-regulated kinases (ERKs) and TNF-alpha. We also noted the involvement of NF-kappaB and c-Jun in response to furan, which are genes that are known to be required for liver regeneration. Furan metabolism by CYP2E1 produces cis-2-butene-1,4-dial (BDA), which is required for ensuing cytotoxicity and oxidative stress. NRF2 is a master regulator of gene expression during oxidative stress and we suggest that chronic NFR2 activity and chronic inflammation may represent critical transition events between the adaptive (regeneration) and adverse (cancer) outcomes. Another objective of this study was to demonstrate the applicability of toxicogenomics data in quantitative risk assessment. We modeled benchmark doses for our transcriptional data and previously published cancer data, and observed consistency between the two. Margin of exposure values for both transcriptional and cancer endpoints were also similar. In conclusion, using furan as a case study we have demonstrated the value of toxicogenomics data in elucidating dose-dependent MOA transitions and in quantitative risk assessment.

Academic research paper on topic "Case study on the utility of hepatic global gene expression profiling in the risk assessment of the carcinogen furan"

ELSEVIER

Contents lists available at ScienceDirect

Toxicology and Applied Pharmacology

journal homepage: www.elsevier.com/locate/ytaap

Case study on the utility of hepatic global gene expression profiling in the risk assessment of the carcinogen furan^

Anna Francina Jackson a,b, Andrew Williams a, Leslie Recioc, Michael D. Watersc, Iain B. Lambertb, Carole L. Yauka'*

a Environmental Health Science and Research Bureau, Health Canada, Ottawa K1A 0K9, Canada b Department of Biology, Carleton University, 1125 Colonel By Drive, Ottawa K1S 5B6, Canada c ILS, Inc., P.O. Box 13501, Research Triangle Park, NC 27709, USA

ARTICLE INFO

ABSTRACT

Article history: Received 12 August 2013 Revised 18 October 2013 Accepted 19 October 2013 Available online 29 October 2013

Keywords: Furan

Gene expression Toxicogenomics

Indirect-acting genotoxic carcinogen Non-genotoxic carcinogen Risk assessment

Furan is a chemical hepatocarcinogen in mice and rats. Its previously postulated cancer mode of action (MOA) is chronic cytotoxicity followed by sustained regenerative proliferation; however, its molecular basis is unknown. To this end, we conducted toxicogenomic analysis of B3C6F1 mouse livers following three week exposures to non-carcinogenic (0, 1, 2 mg/kg bw) or carcinogenic (4 and 8 mg/kg bw) doses of furan. We saw enrichment for pathways responsible for cytotoxicity: stress-activated protein kinase (SAPK) and death receptor (DR5 and TNF-alpha) signaling, and proliferation: extracellular signal-regulated kinases (ERKs) and TNF-alpha. We also noted the involvement of NF-kappaB and c-Jun in response to furan, which are genes that are known to be required for liver regeneration. Furan metabolism by CYP2E1 produces cis-2-butene-1,4-dial (BDA), which is required for ensuing cytotoxicity and oxidative stress. NRF2 is a master regulator of gene expression during oxidative stress and we suggest that chronic NFR2 activity and chronic inflammation may represent critical transition events between the adaptive (regeneration) and adverse (cancer) outcomes. Another objective of this study was to demonstrate the applicability of toxicogenomics data in quantitative risk assessment. We modeled benchmark doses for our transcriptional data and previously published cancer data, and observed consistency between the two. Margin of exposure values for both transcriptional and cancer endpoints were also similar. In conclusion, using furan as a case study we have demonstrated the value of toxicogenomics data in elucidating dose-dependent MOA transitions and in quantitative risk assessment.

Crown Copyright © 2013 Published by Elsevier Inc. All rights reserved.

Abbreviations: APAP, Acetaminophen; AP-1, Activator protein 1; ASK1, Apoptosis signaling kinase; BMD, Benchmark dose; BMDL, Lower confidence limit of benchmark dose; BrdU, Bromodeoxyuridine; bw, Body weight; JNK1, c-Jun NH2-terminal kinase 1; CT, Carbon tetrachloride; BDA, cis-2-butene-1,4-dial; Cyp2E1, Cytochrome P450 2E1; EtOH, Ethanol; ERK, Extracellular signal-regulated kinase; GO, Gene ontology; GSH, Glutathione; HCA, Hepatocellular adenoma; HCC, Hepatocellular carcinoma; IPA, Ingenuity Pathway Analysis; IRIS, Integrated Risk Information System; IARC, International Agency for Research on Cancer; c-Jun, Jun proto-oncogene; MAPK, Mitogen-activated protein kinase; MPT, Mitochondrial permeability transition; MOA, Mode of action; NTP, National Toxicology Program; NF-kB, Nuclear factor kappa B; Nrf2/NFE2L2, Nuclear factor (erythroid-derived 2)-like 2; PH, Partial hepatectomy; POD, Point of departure; ROS, Reactive oxygen species; SAPK, Stress-activated protein kinase; TNF, Tumor necrosis factor; TNFR1, TNF receptor 1.

☆ This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

* Corresponding author at: Environmental Health Science and Research Bureau, ERHSD, HECSB, Health Canada, Tunney's Pasture, Bldg. 8 (P/L 0803A), 50 Columbine Driveway, Ottawa, Ontario K1A 0K9, Canada. Fax: +1 613 941 8530.

E-mail addresses: Francina.Jackson@hc-sc.gc.ca (A.F. Jackson), Andrew.Williams@hc-sc.gc.ca (A. Williams), lrecio@ils-inc.com (L. Recio), mwaters@ils-inc.com (M.D. Waters), Iain.Lambert@carleton.ca (I.B. Lambert), Carole.Yauk@hc-sc.gc.ca (C.L. Yauk).

Introduction

Furan is a liver toxicant and rodent hepatocarcinogen that is classified as possibly carcinogenic to humans (group 2B) by the International Agency for Research on Cancer ([ARC, 1995). Furan was first reported in foods over 30 years ago (Maga and Katz, 1979). It is formed during heat-treatment of food, probably through thermal decomposition of carbohydrates, and is commonly found in coffee and canned and jarred foods (including baby food) (Moro et al., 2012). It is also produced during combustion and is therefore found in engine exhaust, wood smoke and tobacco smoke (IARC, 1995). The National Toxicology Program's (NTP) two-year cancer bioassay showed that furan induces hepatocel-lular carcinoma (HCC) and adenoma (HCA) in a dose-dependent manner in B6C3F1 mice and HCC, HCA, cholangiocarcinoma, and mononuclear cell leukemia in F344 rats (NTP, 1993); however, due to mixed results in the standard battery of genotoxicity tests, furan's carcinogenic mode of action (MOA) was ambiguous. Briefly, at the time of the assessment the available data on furan indicated that it is negative in the Ames assay (± S9) but positive in the L5178Y/TK+/- mouse lymphoma assay

0041-008X/$ - see front matter. Crown Copyright © 2013 Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.! 016/j.taap.2013.10.019

( — S9). In addition, furan causes sister chromatid exchanges in vitro in Chinese hamster ovary (CHO) cells (±S9) but not in vivo in B6C3F1 mouse bone marrow cells. Finally, furan exposure induces chromosomal aberrations in both CHO cells (± S9) and B6C3F1 mouse bone marrow cells (NTP, 1993). Subsequent studies to clarify the genotoxicity of furan revealed that furan is not clastogenic, but that its metabolite (cis-2-butene-1,4-dial (BDA)) causes strand breaks in L5178Y tk+/-mouse lymphoma cells at high (often cytotoxic) doses (Kellert et al., 2008). The results of in vivo MN assays for furan in mouse and rat are typically negative (Durling et al., 2007; Leopardi et al., 2010; McDaniel et al., 2012). While some researchers caution against the dismissal of a genotoxic MOA for furan (Cordelli et al., 2010; Neuwirth et al., 2012), recent studies using male transgenic Big Blue® rats (McDaniel et al., 2012), and follow-up studies using the comet and micronucleus assays in male F344 rats (Ding et al., 2012) concluded that any genotoxic action by furan is likely to be a secondary consequence of oxidative stress produced during furan metabolism.

The MOA proposed for furan based on apical (phenotypic) data derived from female B6C3F1 mice (Fransson-Steen et al., 1997; Moser et al., 2009) (supported by a weight of evidence from furan studies conducted in other rodent models) is chronic cytotoxicity and inflammation followed by dysregulated regenerative proliferation, which is also the most common MOA for spontaneous HCC (Nakagawa and Maeda, 2012). Indeed, in the NTP's 13-week study, dose-dependent increases of histopathological markers for cytotoxicity, necrosis and cellular proliferation were observed in B6C3F1 mice and F344 rats (both genders) and, in the two-year study, there was evidence of chronic inflammation, hepatic fibrosis, hyperplasia, degeneration and necrosis in the liver. After two years the NTP researchers reported increased hepatic cancer rates in both B6C3F1 mice and F344 rats (both genders), with a higher sensitivity to furan in rats. Female mice had high hepatic cancer rates at both doses tested (34/50 and 50/50 at 8 and 15 mg/kg bw, respectively) and had a much lower spontaneous cancer rate than the corresponding male mice (26/50 and 7/50 in males and females, respectively). Moser et al. (2009) exposed female B6C3F1 mice to lower doses of furan (0, 0.5,1, 2,4 or 8 mg/kg bw) for three weeks or two years. In their three week study they observed a significant and dose-dependent increase in hepatic cytotoxicity beginning at 1 mg/kg bw (measured by serum alanine aminotransferase, ALT) and a significant increase in hepatocyte proliferation in the 8 mg/kg bw group (measured by bromodeoxyuridine, BrdU, incorporation). They observed tumorigene-sis at 4 and 8 mg/kg bw after two years. This study suggests that furan's point of departure (POD) for carcinogenicity in female B6C3F1 mice lies between 2 and 4 mg/kg bw. Furan's no-observed adverse effect level (NOAEL) in B6C3F1 mice and F344 rats (males and females) has been reported as 0.12 mg/kg bw and 0.03 mg/kg bw, respectively (Gill et al., 2010, 2011). Ultimately, apical studies in both genders of mice and rats support the concept that furan causes cancer via a cytotoxicity and regenerative proliferation MOA.

The ability of furan to induce levels of cytotoxicity sufficient to promote regenerative proliferation is related to the chemical's metabolism. Furan is metabolized into an electrophilic metabolite (BDA) by cytochrome P450 2E1 (CYP2E1) (Kedderis and Held, 1996), which is the only cytochrome P450 whose expression is not receptor-mediated. Instead, CYP2E1 is constitutively expressed and its activity is induced post-transcriptionally by stabilization in the presence of its substrate. Low molecular weight ligands, including ethanol (EtOH), acetaminophen (APAP), carbon tetrachloride (CT), chloroform and furan, increase Cyp2E1's half-life from seven to 32 h (Gonzalez, 2007). In vitro experiments have demonstrated that metabolic activation to BDA is necessary for furan-induced cytotoxicity (Kellert et al., 2008), and inhibition of CYP2E1 is sufficient to prevent cytotoxicity in female B6C3F1 mice and in mouse, rat or human microsomes (Fransson-Steen et al., 1997; Gates et al., 2012). The CYP2E1 catalytic cycle produces reactive oxygen species (ROS) (Gonzalez, 2005; Lu and Cederbaum, 2008) and BDA has been shown to deplete cellular glutathione (GSH) levels in F344 rat

hepatocytes (Carfagna et al., 1993). There is evidence of oxidative stress-induced genomic damage in response to furan exposure including 8-oxo-dG adducts (Hickling et al., 2010) and oxidized purine and pyrimidine bases (Ding et al., 2012) in rat liver. However, it is unclear whether the indirect consequences of genotoxicity mediated via ROS contribute to the MOA of furan.

Established assays for the identification of non-genotoxic carcinogens are a current gap in the standard battery of short-term carcinogenicity tests. Toxicogenomics is the study of changes in gene expression following chemical exposure; induced perturbations in gene expression are related to the MOA of carcinogens. Thus, toxicogenomics is expected to provide an effective tool for the identification of diverse MOAs and may fill an important testing gap for non-genotoxic carcinogens (Waters et al., 2010). MOA information gleaned from toxicogenomic studies of carcinogens has the added benefit of providing mechanistic information, thereby facilitating inter-species comparisons for inferring human risk. In this study we characterize global gene expression profiles in liver tissue taken from female B6C3F1 mice that were sub-chronically exposed to non-carcinogenic (0,1,2 mg/kg bw) and carcinogenic (4 and 8 mg/kg bw) doses of furan. We perform extensive bioinformatics analyses that are anchored in previously published apical (phenotypic) endpoint data (Gill et al., 2011; Moser et al., 2009; NTP, 1993) in order to elucidate the molecular mechanism by which furan causes liver cancer. More broadly, our goal is to use furan as a case study to champion the concept of using sub-chronic exposures and genomic tools to inform risk assessment of non-genotoxic carcinogens.

Methods

Chemical

Furan (CAS no. 110-00-9) (>99% pure) (Sigma-Aldrich Chemical Co., Milwaukee, WI) was mixed with Mazola corn oil to the appropriate concentration. Doses were prepared separately on a volume-to-weight ratio (v:w), were de-aerated with inert gas, and were stored in 8 mL brown glass vials (sealed with plastic closures and modified silicon septa) in the refrigerator for up to 14 days (based on previous reports of furan stability (NTP, 1993)).

Animals

5-6 week old female specific pathogen free B6C3F1 mice were purchased from Charles River Breeding Laboratories (Portage, ME) and were allowed to acclimatize for at least 7 days prior to the start of the study. Feed (NIH-07; Zeigler Brothers, Inc., Gardners, PA) and tap water were available ad libitum up until the time of necropsy. Mice were housed five per cage in polycarbonate cages in a specific pathogen free (SPF) and Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC) accredited facility. All procedures were conducted in compliance with the Animal Welfare Act Regulations (9CFR1-4). Mice were handled and treated according to the guidelines provided in the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals (ILAR, 1996; http://dels.nas.edu/ilar/).

Female mice were dosed with furan in corn oil at 0, 1, 2, 4, or 8 mg/kg bw per day by oral gavage for three weeks (n = 10 per dose). We chose to use female mice because they have a lower spontaneous tumor rate than males (Haseman etal., 1998). n = 4, 5,3,4, or 5 per dose group were treated with 0.02% Bromodeoxyuridine (BrdU; Sigma Chemical Co., St. Louis, MO, USA) in drinking water for 5 days just prior to sacrifice. Upon necropsy, there remained n = 5 mice in each non-BrdU group. Some mice were lost due to early (pre-BrdU treatment) mis-dosing or esophageal puncture. Four hours after their final dosing, mice were anesthetized by CO2 inhalation prior to euthanasia by exsanguination achieved by cutting the caudal vena cava after blood collection. One animal per group was killed and this continued until all mice had been sacrificed;

this occurred over a period of 100 min (beginning at 1 pm). Body weights were recorded. Serum enzymes were not analyzed since these endpoints have been previously examined and reported (Gill et al., 2011; Moser et al., 2009; NTP, 1993). The left, median, right posterior and right anterior lobes of the liver were cut into 3 sections. One section from each lobe was cut into 0.25-0.5 cm3 pieces. These were snap-frozen in liquid nitrogen then stored at or below — 70 °C until processed for RNA extraction. The remaining two sections of liver were fixed for 18 h or three weeks in neutral phosphate-buffered formalin, then in 100% ethanol for 72 h and stored (for use in studies to be published separately).

Immunohistochemistry for BrdU nuclear staining

Hepatocytes undergoing DNA synthesis were identified by BrdU nuclear staining using immunohistochemistry. Briefly, liver sections were stained for BrdU in an automatic staining machine (IntelliPATH, Biocare, Concord, CA). A small section of duodenum was evaluated on each slide as a positive control to ensure adequate BrdU delivery and staining. Deparaffinized sections were incubated with hydrochloric acid for 20 min at 37 °C for antigen retrieval, neutralized in borate buffer, and immersed in trypsin for 3 min at 37 °C for antigen retrieval. The sections were exposed to 10% normal rabbit serum for 20 min at room temperature to inhibit nonspecific binding. The sections were then incubated for 1 h at room temperature with rabbit monoclonal anti-BrdU antibodies (Accurate Chemical and Scientific Corporation, Westbury, NY) diluted 1:1000. Sections were incubated with hydrogen peroxide for 15 min followed by incubation with avidin and biotin (Vector, Burlingame, CA) at room temperature for 15 min each. Afterwards the sections were incubated with rabbit anti-rat IgG (Vector, Burlingame, CA) at a 1:500 dilution for 30 min at room temperature to capture primary antibody. Sections were then incubated with Vector ABC label for 30 min at room temperature. The antigen-antibody reaction was visualized with betazoid diaminobenzidine (Biocare, Concord, CA) for 5 min. The sections were counterstained with CAT hematoxylin (Biocare, Concord, CA) for 2 min.

A liver section from each mouse was microscopically evaluated for hepatocyte cell proliferation. Random fields in the liver were scored without knowledge of animal identity for BrdU incorporation by light microscopy. At least 1000 cells per slide were counted. The percent labeling index was calculated by the total number of immunoreactive BrdU labeled hepatic nuclei divided by the total number of nuclei counted times 100. Representative images from each dose group were captured on an Olympus BX41 microscope and Olympus DP25 camera.

A generalized linear model was applied to the data with the binomial error distribution using the R software (R-Development-Core-Team, 2010). The glht function in the multcomp library (Hothorn et al., 2008) for testing general linear hypotheses was used to test for differences between exposed and controls. The P-values from Wald's test were adjusted using the Dunnett's multiple comparison procedure (Dunnett, 1955). Estimates were back-transformed to obtain fold change estimates and the delta method was applied to obtain the estimate for the standard errors.

RNA extraction

RNA was extracted from ~100 mg frozen liver tissue using the RNeasy Midi RNA Extraction kit (Qiagen, Mississauga, ON, Canada). An Omni tissue homogenizer with a disposable 7 mm Omni generator probe was used (Omni #34750, Omni International, Marietta, GA). RNA was quantified using a NanoDrop Spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA) and qualified using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Mississauga, ON, Canada) and stored at — 80 °C.

Microarray

RNA was extracted from each of the non-BrdU mice for all furan doses (0, 1, 2, 4, 8 mg/kg bw; n = 5 per dose group). Sample RNA (200 ng) was used together with a mouse universal reference RNA (Stratagene by Agilent Technologies Inc., Mississauga, ON, Canada) to synthesize, amplify and label cRNA using the Low Input Quick Amp Labeling Kit (Agilent Technologies Inc., Mississauga, ON, Canada). Labeled cRNA was purified using the RNeasy Mini Kit (Qiagen, Mississauga, ON, Canada). Amplification and labeling efficiency of cRNA were quantified using a NanoDrop spectrophotometer. Hybridization mixes were prepared using the Hi-RPM Gene Expression Hybridization Kit (Agilent Technologies Inc., Mississauga, ON, Canada). 300 ng of Cy3-labeled reference RNA and 300 ng Cy5-labeled sample cRNA were hybridized on SurePrint G3 Mouse GE 8 x 60 K microarrays (Agilent Technologies Inc., Mississauga, ON, Canada) at 65 °C for 17 h at 10 rpm. Slides were washed according to the manufacturer's specifications with Gene Expression Wash Buffers 1 and 2 (Agilent Technologies Inc., Mississauga, ON, Canada), scanned using an Agilent G2505B scanner at 5 |am resolution. Data were extracted using Agilent Feature Extraction Software, version 11. The complete dataset is available through the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) accession number GSE48644.

A block design, treating the slide as the blocking effect was employed to determine global differential gene expression. Median signal intensities were normalized using LOWESS (Bolstad et al., 2003) in R (R-Development-Core-Team, 2010). Probes with technical replicates were then averaged. Differential gene expression was determined using the microarray analysis of variance (MAANOVA) library (Wu et al., 2003). The Fs statistic (Cui et al., 2005), a shrinkage estimator, was used to determine gene-specific treatment effects, and the associated P values were estimated using the permutation method (30,000 permutations with residual shuffling). The P-values were adjusted for multiple comparisons using the false discovery rate (FDR) approach (Benjamini and Hochberg, 1995). Fold changes were estimated using least square means of each pairwise comparison. Genes having an FDR-adjusted P < 0.05 and a fold change > ±1.5 were deemed differentially expressed. Upon removal of outliers (arrays with high background), the final sample sizes used for gene expression analysis were n = 5, 4, 5, 4, and 5 for 0, 1, 2, 4, and 8 mg/kg bw furan dose groups, respectively.

Bioinformatics

Ingenuity Pathway Analysis (IPA). IPA (http://www.ingenuity.com/ products/ipa) was used to identify molecular pathways that were affected by furan treatment and to predict activated upstream regulators. The expression and significance cut-offs applied to the data were fold change > ±1.5 and FDR-adjusted P < 0.05. Enrichment of a canonical pathway is determined in IPA based on the number of differentially expressed genes in the data set that also appear in the pathway. The significance threshold for canonical pathways is P < 0.05 where P is calculated using a right-tailed Fisher's exact test by IPA. Significance of predicted upstream regulators was calculated by IPA and quantified using a Z-score (we considered Z > 1.9 as activated and Z < —1.9 as inhibited), which takes into account the number of downsteam genes that are differentially expressed in the gene list.

Disease prediction and chemical profile comparison. Gene expression changes induced by 8 mg/kg bw furan were mined against public geno-mic data repositories for chemical- and disease-induced changes in gene expression using NextBio (http://nextbio.com). Pairwise gene signature correlations and rank-based enrichment statistics are employed in the NextBio score calculations, where the chemical or disease with the highest similarity will be given a score of 100 and the rest are being normalized accordingly (Kupershmidt et al., 2010).

PCA in R

Principle component analysis was conducted on the relative signal intensities for the significant probes based on the correlation matrix in R using the prcomp() function (Becker et al., 1998; Mardia et al., 1979; Venables and Ripley, 1999). The first three components were then plotted using the scatterplot3d() function in the scatterplot3d library (Ligges and Machler, 2003).

Cluster analysis

Data was obtained from Gene Expression Ominibus (http://www. ncbi.nlm.nih.gov/geo/) and European Bioinformatics Institute (http:// www.ebi.ac.uk/arrayexpress/). All data sets that used the Affymetrix platform (E-MEXP-82; GSE13149; GSE18858; GSE20427; GSE26538) were normalized using Robust Multi-array Average, RMA (Irizarry et al., 2003), using the ReadAffy() function in the affy library in R (Gautier et al., 2004). Two-color platform (Furan; BaP) background subtracted signal intensities were normalized using the LOWESS approach and background subtracted signal intensities for GSE35934 and GSE4874 (with a dye adjustment) were normalized using cyclic LOWESS (Bolstad et al., 2003).

Normalized probe intensities with common Gene Symbols were then collapsed using the median normalized signal intensity. Based on Gene Symbol, all normalized data were merged together yielding 3190 common Gene Symbols. A heatmap using these data was produced using 1-Pearson correlation as the distance metric with average linkage.

Benchmark dose (BMD) calculations

Apical endpoint data. BMD and BMDL (lower confidence limit of BMD) values were modeled for previously reported apical endpoint measurements (Gill et al., 2011; NTP, 1993) using the United States Environmental Protection Agency's (US EPA) Benchmark Dose Software (BMDS) version 2.3.1 (http://www.epa.gov/ncea/bmds/) (Davis et al., 2011). Data were pre-screened for homogeneity of variance. Continuous data were run against five models (Exponential, Hill, Linear, Polynomial and Power) and dichotomous data were run against nine models (Gamma Multi, Logistic, LogLogistic, LogProbit, Multistage, Multistage Cancer, Probit, Weibull and Quantal-Linear). The best model for each data set was chosen based on visual inspection of curve fitting, lowest Akaike's Information Criterion (AIC, which measures the relative goodness of fit) value, goodness of fit P > 0.05, and scaled residuals between + 2.0 and — 2.0. The AIC was used as the deciding factor between equally well fitting models. These selection criteria were applied when values generated from each model were within 3-fold of each other, otherwise the lowest value was used.

Gene expression data. BMDExpress (Yang et al., 2007) was applied to calculate BMD (mean and median) and BMDL (mean and median) values for genes, IPA pathways and GO terms. Datasets were pre-filtered with FDR P < 0.05 and fold change > 1.5 in at least one dose before importing into BMDExpress. Five models were compared (Hill, Power, Linear, Polynomial 2° and Polynomial 3°) and the best was chosen based on: (1) a nested Chi-square test, with cut-off of 0.05; (2) the lowest AIC; (3) maximum iterations of250; (4) confidence level of 0.95; and (5) benchmark response (BMR) of 1.349 (number of standard deviation defining BMD), which corresponds to an excess risk of 10%. The power model had a power restriction of > 1. Selection of Linear and Polynomial 2° was based on choosing a model that describes the data with the least complexity. A Hill model was excluded if the "k" parameter of the model was less than 1/3 of the lowest positive dose as per Black et al. (2012). The datasets were mapped to Ingenuity Core Pathways (http://www.ingenuity.com) using the Defined Category Analyses feature. The pathway dataset was downloaded from Ingenuity Pathway

Analysis (December 13, 2012) and array annotation dataset, MM_ WG_GPL10333, was downloaded from NCBI Gene Expression Omnibus (October 18, 2011 version).

Margin of exposure

Margin of exposure (MoE) is a metric used to infer public health risk. Typically a MoE of 10,000 or greater indicates that a compound is not a public health concern (Benford et al., 2010; EFSA, 2005); however 1000 is the benchmark used for furan since it is a relatively well characterized chemical and therefore has a lower uncertainty factor associated with it (IRIS, 2012). In this study MoEs were calculated by dividing the relevant BMDL (^g/kg bw) by the actual environmental exposure (^g/kg-bw/ day). Average (50th percentile) and high (90th percentile) environmental exposure estimates for USA or Europe, and for infants or combined child/adult populations were taken from Carthew et al. (2010). Apical MoEs were calculated for HCC and HCA, with tumor BMDL values calculated using cancer rates published by Moser et al. (2009) and NTP (1993). Transcriptional MoEs were calculated for the NRF2 Oxidative Stress Response IPA pathway and for a 'General transcriptional response,' the latter was calculated using the average BMDL for all enriched pathways with at least five molecules and P < 0.05. We use the BMDL value obtained from the best model (as described above) for MoE calculations.

Results

Cellular proliferation (BrdU), cytotoxicity and cancer

Treatment-related effects of 0,1,2,4, or 8 mg/kg bw furan for three weeks on mice (including measurement of serum enzymes and histo-pathological markers of cytotoxicity and inflammation) have been previously reported (Gill et al., 2011; Moser et al., 2009; NTP, 1993) and therefore were not repeated in this study. Body weight was not affected by furan exposure (Suppl. Table S1). BrdU incorporation was observed to increase with increasing furan dose. Labeling indices for BrdU of 2.3, 2.4, 2.9, 5.2, and 7.6% at 0,1, 2,4 and 8 mg/kg bw furan, respectively, were observed, with a significant increase in BrdU incorporation at 8 mg/kg bw (P = 0.00413) (Fig. 1a). Similar results for BrdU incorporation were previously reported in Moser et al. (2009).

Previous studies have reported increases in cytotoxicity and proliferative histopathological markers (Gill et al., 2011; Moser et al., 2009; NTP, 1993) as well as cancer rates (Moser et al., 2009; NTP, 1993) in response to furan. We modeled dose responses for these apical endpoints and our own BrdU incorporation data using BMDS (Table 1) and observed that BMD/BMDL variability between studies was low for cancer endpoints, but was quite large for cytotoxicity and proliferation markers.

Gene expression analysis

Following furan treatment, the number of differentially expressed genes increased in a dose-dependent manner with 2, 17, 27 and 339 genes (2,18,43 and 441 probes) in the 1,2,4, or 8 mg/kg bw dose groups, respectively (Fig. 1b). While fold changes ranging between +20.65 and — 6.60 were observed, the majority of significantly changed probes had fold changes of ± 1.5-2.0. There was a complete overlap between the differentially expressed genes in the 4 and 8 mg/kg bw dose groups, and only one shared gene between the 2 and 4 mg/kg dose groups thus indicating a clean transition in the nature of the transcriptional response. A PCA analysis revealed a dose-related clustering of hepatic expression profiles (Fig. 1c). We saw up-regulation of a number of genes consistent with our proposed MOA for furan including genes involved in oxidative stress response, xenobiotic metabolism, inflammation, cell cycle arrest/cell death and cell survival and growth (Table 2; Suppl. Table S2). Using BMDExpress, we calculated median BMD/BMDL values for individual

a 12 10

0) 8 >

■<ö 6 o

control * 1 mg/kg 2 mg/kg

1.0 2.0 4.0

Furan dose (mg/kg bw)

ra 300 TD

œ <n <n œ

Si 250

1 200 c

4 mg/kg S mg/kg

2.0 4.0

Furan dose (mg/kg bw)

-30 -20 -10

10 20 30 40

Fig. 1. Dose response due to furan. (a) Dose-dependent increase in liver cell proliferation upon exposure to furan based upon BrdU incorporation (+/—stdev; *P = 0.00413). (b) Number of furan-responsive differentially expressed genes in each dose group and Venn diagrams represent the overlap of individual furan-responsive genes in each dose group. (C) Principal component analysis (PCA) of furan dose groups, where circles represent individual mice and colors indicated dose groups (royal blue, 8 mg/kg bw; green, 4 mg/kg bw; pink 2 mg/kg bw; red, 1 mg/kg bw; aqua blue, 0 mg/kg bw).

probes. We observed that many BMD/BMDL values for genes thought to be important for the furan MOA were close to the 2-4 mg/kg transitional dose range, and this trend persisted regardless of the magnitude of the fold change of the probes used in the modeling (Table 2, Suppl. Table S3).

Cluster analysis

Cluster analysis was performed using gene lists from 1, 2, 4 and 8 mg/kg bw furan dose groups and gene expression data sets from GEO (http://www.ncbi.nlm.nih.gov/geo/) from spontaneous hepatocel-lular carcinoma, liver regeneration following partial hepatectomy, liver inflammation, and chemical exposure to hepatocarcinogens, non-hepatocarcinogens and Cyp2E1 ligands (APAP and CT) (Fig. 2, Suppl. Table S4). Furan dose groups clustered closely together, with a split between the lower dose groups and the 8 mg/kg dose group. The 8 mg/kg bw furan dose group clustered most closely with a hepatocellular carcinoma (HCC) data set, followed by liver regeneration 38 and 48 h post-PH. The 1, 2 and 4 mg/kg furan dose groups clustered with a dataset representing NF-kappaB-mediated liver inflammation. These cluster together with a large group of hepatocarcinogens and a small group of Cyp2E1 substrates.

Pathways and upstream regulators of transcription analyses

Pathway analysis (using IPA) was employed to provide information regarding the furan molecular MOA. The number of enriched pathways increased in a dose-dependent manner. There were no enriched pathways at the lowest dose (1 mg/kg bw) and pathway enrichment at the next two doses (2 and 4 mg/kg bw) was always based on three or fewer genes, therefore we focused on pathways

enriched at the highest dose (8 mg/kg bw). Pathways relevant to the furan MOA included NRF2 Oxidative Stress Response, Glutathione Metabolism, Xenobiotic Metabolism, Death Receptor Signaling, TNF-alpha/ASK1/JNK1 Apoptosis Signaling (within the 14-3-3 Signaling pathway), Erk/MAPK Signaling and Breast Cancer Regulation by Stathminl (Fig. 3A). In addition, enrichment of the p53 Signaling pathway indicated that the DNA damage response could be involved. Other pathways, including the Xenobiotic Metabolism pathway and the Aryl Hydrocarbon Receptor pathway, were also enriched in response to furan. While the Xenobiotic Metabolism pathway comprises CAR, PXR, AhR, PPAR and Nrf2 signaling, Nrf2 was the only one among these that was appreciably activated. Closer inspection of the Aryl Hydrocarbon pathway showed that enrichment of this pathway was almost exclusively the result of an increase in the expression of the estrogen receptor (Esrl), and not the AhR. These findings are consistent with previous studies showing that furan does not undergo receptor-mediated metabolism. BMD modeling for individual pathways demonstrated that BMDL values fell around the 2-4 mg/kg transitional dose range (Fig. 3A). For all IPA pathways (that could be modeled and that had at least five genes with significant changes in expression) the average pathway BMDL mean was 2.20 mg/kg bw (range 1.504.40 mg/kg bw).

Upstream regulators of transcription were predicted using IPA software. There were no predicted upstream regulators for the 1 or 2 mg/kg bw dose groups. TNF-alpha was the only predicted upstream regulator in the 4 mg/kg bw dose group. There were many predicted upstream regulators for the 8 mg/kg dose group; in particular, activation of reactive oxygen species, hydrogen peroxide, other Cyp2E1 ligands (CT, APAP, EtOH), cytokines (including TNF-alpha), growth factors,

Table 1

BMD modeling of previously published apical endpoint data for furan exposures in B6C3F1 mice. For cytotoxicity and proliferative markers Gill et al. (2011) exposed male and female mice to 0,0.03,0.12,0.5,2 and 8 mg/kg furan for 13 weeks and the NTP (1993) exposed female mice to 0,2,4,8,15,30 and 60 mg/kg bw furan for 13 weeks. For cancer studies, Moser et al. (2009) exposed female mice 0, 0.5, 1, 2, 4 and 8 mg/kg bw furan for two years and the NTP (1993) exposed female mice to 0, 8 and 15 mg/kg furan for two years; cancer incidences in these studies are in brackets.

Apical endpoint BMD BMDL

(mg/kg bw) (mg/kg bw)

Cytotoxicity markers (Gill et al., 2011)

Serum alanine transaminase (ALT) 1.84 1.35

Hepatocyte apoptosis (caudate lobe) 0.43 0.11

Hepatocyte apoptosis (other lobes) 1.37 0.28

Kuffer cell pigmentation and inflammation 1.20 0.28

Cytotoxicity markers (NTP, 1993)

Hepatocyte degeneration 19.85 14.67

Hepatocyte necrosis 23.50 14.44

Kupffer cell pigmentation 23.50 14.44

Cholangiofibrosis 27.16 17.57

Proliferation markers (this study)

BrdU labeling index 1.10 0.67

Proliferation markers (Gill et al., 2011)

Hepatocyte basophilia 0.60 0.30

Biliary tract hyperplasia 5.81 1.05

Proliferation markers (NTP, 1993)

Biliary tract hyperplasia 24.58 14.54

Hepatocyte cytomegaly 19.85 14.67

Cancer (NTP, 1993)

Hepatocellular adenoma (5/50, 31/50, 48/50) 2.78 0.92

Hepatocellular carcinoma (2/50, 7/50, 27/50) 6.61 5.30

Hepatocellular adenoma or carcinoma 6.88 4.20

(7/50,34/50,50/50)

Cancer (Moser et al., 2009)

Hepatocellular adenoma 2.34 1.34

(3/36,4/72,4/53,4/41,11/36,25/39)

Hepatocellular carcinoma 2.60 1.57

(0/36,4/72, 2/53,1/41, 2/36,11/39)

Hepatocellular adenoma or carcinoma 5.14 4.22

(3/36,8/72, 6/53, 5/41,12/36,29/39)

mitogen-activated protein kinases (MAPKs, including SAPKs and ERKs), transcription factors (including Nrf2 (Nfe2l2) and NF-kB1), and hormones (including estrogen and thyroid hormone) were predicted (Fig. 3A).

Disease and chemical prediction

Meta-analyses comparing the changes in gene expression observed in the 8 mg/kg bw dose group to other gene expression studies were conducted using the NextBio Human Disease Atlas and revealed similarities to the following: Injury to Liver (97.4), Liver Regeneration (91.7), Hepatic Fibrosis (82.25), Hepatocellular Dysplasia (80.9), Liver Cancer (73.0), Inflammatory Disease of the Liver (68.5), and Cirrhosis of the Liver (52.07). In addition, using the NextBio Chemical Atlas, similarities were observed with chloroform (47.42), CT (43.02), hydrogen peroxide (40.75), and APAP (39.12). Interestingly, 8 mg/kg bw furan was found to produce gene expression changes that were most similar to the insecticide malathion (100) (Suppl. Tables S5 and S6).

GO terms

There were many, redundant GO terms enriched for that were relevant to furan 's primary MOA. Interestingly, mean BMDL values approached the expected ~2-4 mg/kg range, including: kinase activity

(GO: 0016301; BMDL 1.89 mg/kg bw), regulation of transcription, DNA dependent regulation of transcription (GO: 0006355; BMDL 2.14 mg/ kg bw), regulation of protein modification process (GO: 0031399; BMDL 1.90 mg/kg bw), regulation of apoptosis (GO: 0042981; BMDL 1.90 mg/kg bw), negative regulation of apoptosis (GO: 0043066; BMDL 2.18 mg/kg bw), and regulation of gene expression (GO: 0010468; BMDL 2.04 mg/kg bw). One gentoxicity-relevant GO term, response to DNA damage stimulus (GO: 0006974; BMDL 2.06 mg/kg bw), was also enriched for (full list: Suppl. Table S7).

Margin of exposure (MoE)

Two MoEs each were calculated for apical (HCC and HCA) and transcriptional (NRF2 Oxidative Stress Response pathway and General Transcriptional Response) endpoints. Estimates of risk calculated for apical or transcriptional endpoints were generally quite consistent (Table 3).

Discussion

Global transcriptomic analyses on liver tissue collected from mice exposed sub-chronically to non-carcinogenic (0, 1,2 mg/kg bw) and carcinogenic (4 and 8 mg/kg bw) doses of furan were carried out to provide qualitative and quantitative mechanistic insight into how furan causes liver cancer. The most common MOA for the development of liver cancer is chronic cytotoxicity and inflammation leading to regenerative proliferation (Nakagawa and Maeda, 2012). This MOA is also hypothesized to underlie carcinogenicity induced by other chemicals including the group 2B non-genotoxic carcinogens CT (Manibusan et al., 2007) and chloroform (Larson et al., 1994; Templin et al., 1996). Furan is believed to cause cancer by this MOA (Fransson-Steen et al., 1997; Moser et al., 2009); however, it is not clear how furan stimulates this process on a molecular level. Our case study supports the hypothesis that furan's primary MOA is cytotoxicity followed by cellular proliferation and regeneration. The analysis also supported potential genotoxicity occurring at the high dose via indirect mechanisms (i.e.: ROS generation). Our discussion focuses primarily on the non-genotoxic MOA; an important outcome of this case study is the demonstration that sub-chronic in vivo toxicogenomics studies can be used to discern MOA information for non-genotoxic carcinogens, and that these data can be used to inform quantitative human health risk assessment. This is an important finding because non-genotoxic carcinogens aren't detected by short-term mutagenicity tests, and the length and cost of the two-year rodent bioassay is prohibitive.

In this study we observed an increase in cellular proliferation after three weeks of furan exposure at both of the carcinogenic doses examined (Fig. 1a), which confirmed the induction of the furan MOA and was consistent with previously published results (Moser et al., 2009). Although no measurements for cytotoxicity were taken in our study, this in vivo response to furan is well-characterized elsewhere (Gill et al., 2011; Moser et al., 2009; NTP, 1993), as is the in vitro cytotoxicity of BDA (Kellert et al., 2008). We also noted that all of the differentially expressed genes at 4 mg/kg bw were contained within the 8 mg/kg bw gene list and were distinct from those observed at 2 mg/kg bw, indicating that a shift from an adaptive/stress response to an adverse/cancer response had presumably occurred (Fig. 1b). The molecular MOA of the adverse response to furan was best characterized by relying on the carcinogenic 8 mg/kg bw dose group, since it had the most populated list of differentially expressed genes, it clustered with liver regeneration and liver cancer datasets, and clustered separately from the lower dose groups (Fig. 1c; Fig. 2). This separate clustering suggests that the 8 mg/kg bw livers had progressed further toward the HCC and HCA disease molecular phenotype during the three week exposure than the (also carcinogenic) 4 mg/kg bw dose group. Another possibility is that the MOA at 4 mg/kg is different from that

Genes that are differentially expressed (fold change > ±1.50; P < 0.05) in response to 8 mg/kg bw furan and relevant to the furan mode of action. BMD/BMDL values for individual genes are listed (unless the software was unable to model the dose response of their gene expression, or the value exceeded 8 mg/kg).

Probe ID GenBank Gene symbol Gene name Fold change BMD BMDL

accession (8 mg/kg furan) (mg/kg bw) (mg/kg bw)

Oxidative stress response and xenobiotic metabolism

A_52_P664506 NM_018811 Abhd2 Abhydrolase domain containing 2 1.56 1.42

A_55_P2177233 NM_026179 Abhd5 Abhydrolase domain containing 5 1.61 2.16 1.65

A_52_P273821 NM_026179 Abhd5 Abhydrolase domain containing 5 1.68 1.38 0.70

A_55_P2163098 NM_134066 Akr1c18 Aldo-keto reductase family 1, member C18 2.78 4.37 3.00

A_55_P2154416 NM_030558 Car15 Carbonic anhydrase 15 2.10 2.05 1.57

A_51_P455647 NM_009801 Car2 Carbonic anhydrase 2 2.12 3.08 2.25

A_66_P108152 NM_007620 Cbr1 Carbonyl reductase 1 1.63 2.92 2.15

A_51_P481159 NM_173047 Cbr3 Carbonyl reductase 3 2.53 3.14 2.29

A_55_P2124712 NM_145603 Ces2 Carboxylesterase 2 3.08 2.68 2.00

A_52_P318361 NM_145603 Ces2 Carboxylesterase 2 3.62 2.72 2.02

A_51_P179919 NM_172759 Ces5 Carboxylesterase 5 1.91 7.36 4.51

A_55_P1954718 NM_007805 Cyb561 Cytochrome b-561 1.88

A_55_P2140107 NM_009995 Cyp21a1 Cytochrome P450, family 21, subfamily a, polypeptide 1 1.52

A_55_P2007326 NM_009997 Cyp2a4 Cytochrome P450, family 2, subfamily a, polypeptide 4 1.51

A_51_P137452 NM_013809 Cyp2g1 Cytochrome P450, family 2, subfamily g, polypeptide 1 -1.80 7.87 4.55

A_51_P482051 NM_007820 Cyp3a16 Cytochrome P450, family 3, subfamily a, polypeptide 16 1.60 1.05

A_55_P2001780 NM_001105159 Cyp3a41b Cytochrome P450, family 3, subfamily a, polypeptide 41B 1.67

A_52_P164161 NM_020010 Cyp51 Cytochrome P450, family 51 1.62 0.97 0.16

A_55_P2002577 NM_010145 Ephx1 Epoxide hydrolase 1, microsomal 1.70 2.79 2.07

A_55_P2002578 NM_010145 Ephx1 Epoxide hydrolase 1, microsomal 1.94 3.24 2.00

A_51_P365019 NM_010295 Gclc Glutamate-cysteine ligase, catalytic subunit 1.61 5.40 3.95

A_55_P2032946 NM_008181 Gsta1 Glutathione S-transferase, alpha 1 (Ya) 15.28 4.98 3.42

A_55_P2170454 NM_008182 Gsta2 Glutathione S-transferase, alpha 2 (Yc2) 3.63 3.25 1.70

A_55_P2062190 NM_010358 Gstm1 Glutathione S-transferase, mu 1 1.66 6.33 4.97

A_55_P1957038 NM_181796 Gstp2 Glutathione S-transferase, pi 2 1.56 7.74 5.79

A_51_P263965 NM_010442 Hmox1 Heme oxygenase (decycling) 1 2.09 7.61 4.20

A_55_P2029687 NM_010442 Hmox1 Heme oxygenase (decycling) 1 1.66 6.09 3.84

A_51_P424338 NM_008706 Nqo1 NAD(P)H dehydrogenase, quinone 1 1.71 3.08 2.25

A_51_P161354 NM_144907 Sesn2 Sestrin 2 1.80 2.63 1.97

A_51_P243755 NM_011388 Slc10a2 Solute carrier family 10, member 2 3.61 1.09

A_55_P1969506 NM_009198 Slc17a1 Solute carrier family 17 (sodium phosphate), member 1 1.52

A_52_P286520 NM_018824 Slc23a2 Solute carrier family 23 (nucleobase transporters), member 2 - 1.52 4.61 3.13

A_55_P1964752 NM_194333 Slc23a3 Solute carrier family 23 (nucleobase transporters), member 3 2.27 2.55 1.91

A_51_P514405 NM_019741 Slc2a5 Solute carrier family 2 (facilitated glucose transporter), member 5 - 1.92 4.88 3.27

A_52_P193194 NM_001001321 Slc35d2 Solute carrier family 35, member D2 1.66 7.84 4.66

A_55_P2033120 NM_029688 Srxn1 Sulfiredoxin 1 homolog (Saccharomyces cerevisiae) 1.73 4.65 2.95

A_51_P353895 NM_026935 Sult1c2 Sulfotransferase family, cytosolic, 1C, member 2 - 1.81 3.40 2.45

A_55_P2129449 NM_020565 Sult3a1 Sulfotransferase family 3A, member 1 - 2.55 4.16 3.31

A_51_P320614 NM_001042523 Txnrd1 Thioredoxin reductase 1 1.52 1.76 0.66

A_55_P2006236 NM_009466 Ugdh UDP-glucose dehydrogenase 1.82 7.42 4.16

A_55_P2057577 NM_145079 Ugt1a6a UDP glucuronosyltransferase 1 family, polypeptide A6A 1.51 5.23 3.44

A_51_P163578 NM_172881 Ugt2b35 UDP glucuronosyltransferase 2 family, polypeptide B35 2.26 7.66 4.79

A_55_P1967350 NM_172881 Ugt2b35 UDP glucuronosyltransferase 2 family, polypeptide B35 2.69 7.69 5.11

Inflammation

A_55_P2016462 NM_021274 Cxcl10 Chemokine (C-X-C motif) ligand 10 2.55 1.67 1.29

Cell cycle arrest/cell death

A_55_P2145804 NM_026531 Aen Apoptosis enhancing nuclease 2.46 1.04 0.82

A_55_P2141860 NM_026531 Aen Apoptosis enhancing nuclease 2.34 1.02 0.81

A_55_P2002849 NM_175178 Aifm3 Apoptosis-inducing factor, mitochondrion-associated 3 - 1.62 3.23 2.35

A_55_P2137406 NM_007527 Bax BCL2-associated X protein 1.73 2.45 1.84

A_51_P117794 NM_007546 Bik BCL2-interacting killer - 1.68 5.20 3.43

A_55_P2029106 NM_138313 Bmf BCL2 modifying factor -2.30 6.34 4.08

A_51_P363947 NM_007669 Cdkn1a/p21 Cyclin-dependent kinase inhibitor 1a (p21) 4.64 2.16 1.31

A_55_P1969131 NM_178373 Cidec Cell death-inducing DFFA-like effector c 3.08 4.51 3.08

A_52_P311853 NM_030143 Ddit4l DNA-damage-inducible transcript 4-like 2.45 4.31 2.28

A_51_P336385 NM_175551 Dido1 Death inducer-obliterator 1 1.63 3.15 2.30

A_51_P189361 NM_027950 Osgin1 Oxidative stress induced growth inhibitor 1 2.38 3.62 2.58

A_55_P2027836 NM_020275 Tnfrsf10b/DR5 Tumor necrosis factor receptor superfamily, 10b/death receptor 5 1.92 2.33 1.76

A_55_P2027879 NM_172571 Fbf1 Fas (TNFRSF6) binding factor 1 2.28 3.72 2.64

Both (cell cycle arrest/cell death and cell survival/growth)

A_52_P452689 NM_007498 Atf3 Activating transcription factor 3 2.47 3.87 2.64

A_51_P502614 NM_026268 Dusp6 Dual specificity phosphatase 6 2.59 2.61 1.95

A_51_P263246 NM_008748 Dusp8 Dual specificity phosphatase 8 1.94 2.57 1.93

A_55_P2158990 NM_010591 Jun/c-Jun Jun oncogene 4.16 2.69 1.69

Cell survival/growth

A_55_P1979728 NM_009716 Atf4 Activating transcription factor 4 1.58 2.66 1.98

A_55_P1983773 NM_001012273 Birc5/survivin Baculoviral IAP repeat-containing 5/survivin 1.54* 6.32 3.94

A_66_P111562 NM_007631 Ccnd1 Cyclin D1 1.53

A_52_P612803 NM_009831 Ccng1 Cyclin G1 2.00 1.63 0.48

(continued on next page)

Table 2 (continued)

Probe ID GenBank Gene symbol Gene name Fold change BMD BMDL

accession (8 mg/kg furan) (mg/kg bw) (mg/kg bw)

A_51_P372550 NM_026770 Cgref1 Cell growth regulator with EF hand domain 1 1.65 4.14 2.88

A_52_P296913 NM_031396 Cnnm1 Cyclin M1 1.61 7.77 4.94

A_55_P1976127 NM_007900 Ect2 Ect2 oncogene 1.51 5.91 3.76

A_52_P237077 NM_007956 Esr1/ER1a Estrogen receptor 1 (alpha) 1.52 1.06

A_52_P235347 NM_020013 Fgf21 Fibroblast growth factor 21 3.25 7.63 5.23

A_55_P1960735 NM_011819 Gdf15 Growth differentiation factor 15 4.13 1.97 1.51

A_51_P220806 NM_008110 Gdf9 Growth differentiation factor 9 3.14 4.72 2.94

A_51_P447545 NM_008341 Igfbp1 Insulin-like growth factor binding protein 1 1.86 7.88 4.93

A_51_P260300 NM_010547 Ikbkg/NEMO Inhibitor of kappab kinase gamma 1.80 4.05 2.83

A_55_P2050652 NM_001136067 Ikbkg/NEMO Inhibitor of kappab kinase gamma 1.76 2.78 2.06

A_52_P108346 NM_010849 Myc/c-Myc Myelocytomatosis oncogene 1.84 4.48 3.06

A_55_P2007470 NM_008808 Pdgfa Platelet derived growth factor, alpha 1.89 1.03 0.72

A_55_P2125588 NM_008808 Pdgfa Platelet derived growth factor, alpha 1.88 1.80 1.39

A_55_P2092909 NM_019713 Rassf1 Ras association (ralgds/AF-6) domain family member 1 1.66 1.91 1.47

A_51_P111164 NM_172612 Rnd1 Rho family gtpase 1 — 1.59 5.70 3.59

A_51_P404377 NM_009708 Rnd2 Rho family gtpase 2 1.89 1.13 0.45

A_55_P2184009 NM_009708 Rnd2 Rho family gtpase 2 1.73 4.24 2.93

A_55_P2031045 NM_011489 Stat5b Signal transducer and activator of transcription 5B —1.51 3.93 2.76

A_55_P2068663 NM_019641 Stmn1 Stathmin 1 1.53 4.94 3.30

A_51_P131408 NM_013749 Tnfrsf12a/TweakR Tumor necrosis factor receptor superfamily, 12a/TNF related weak inducer of apoptosis receptor 2.54 3.70 2.35

A_55_P2116496 NM_011654 Tuba1b Tubulin, alpha 1B 1.65 3.86 2.72

A_55_P2108837 NM_009448 Tuba1c Tubulin, alpha 1C 2.40 1.42 0.54

A_55_P2064547 BC022182 Tuba1c Tubulin, alpha 1C 2.36 1.53 0.55

A_52_P676271 NM_009447 Tuba4a Tubulin, alpha 4A 1.71 6.35 3.95

A_51_P490023 NM_009450 Tubb2a Tubulin, beta 2A 4.30 2.69 2.00

A_55_P2077783 NR_003964 Tubb2a-ps2 Tubulin, beta 2a, pseudogene 2 2.48 2.76 2.05

A_55_P2153292 NM_146116 Tubb2c Tubulin, beta 2C 2.86 1.41 1.10

A_55_P1959703 XM_001473123 Tubb2c-ps1 Tubulin, beta 2c, pseudogene 1 2.71 2.02 1.54

A_51_P421140 NM_026473 Tubb6 Tubulin, beta 6 1.71 1.96 1.51

* P = 0.0776.

at 8 mg/kg; however, there is limited evidence from our data or others to support this hypothesis. We propose a molecular mechanism for furan that accounts for furan-induced increases in liver cytotoxicity, inflammation and proliferation, leading to tumors, and is anchored in the furan, Cyp2E1, liver regeneration, and spontaneous or chemically induced HCC literature. Briefly, the molecular initiating event is furan metabolism by CYP2E1. The resultant production of BDA and ROS brings about changes in the cell that trigger cytotoxicity and oxidative stress, which leads to inflammation. Regenerative proliferation follows liver injury; however, when paired with chronic elevation of NRF2 oxidative stress response and chronic inflammation, a transition point is created that tips the balances toward sustained cellular proliferation. Tumorigenicity occurs when the dose of furan results in a level of BDA and ROS production that exceeds the rate at which the cell can neutralize them, thereby requiring chronic elevation of these (normally protective) cellular pathways. The result is sustained and dysregulated cellular proliferation that ultimately leads to cancer (Figs. 3 and 4).

In the larger context of spontaneous and chemically-induced liver cancer, many potential molecules and signaling pathways are hypothesized to be responsible for carcinogenic transformation during chronic liver injury. Prominent candidates, that are also implicated in the present study, include tumor necrosis factor alpha (TNF-alpha, which is a cytokine), c-Jun NH2-terminal kinase 1 (jNK1, which is a SAPK), and nuclear factor kappa B (NF-kB) and c-Jun proto-oncogene/activator protein 1 (c-Jun/AP-1), which are transcription factors. Indeed, previous studies of liver regeneration have also identified TNF-alpha, JNK1, c-Jun and NF-kB as required for liver regeneration. In this study we obtained support for their involvement, as well as many other gene products and pathways including death receptor 5 (DR5), NRF2 oxidative stress response, extracellular signal-regulated kinase (ERK) signaling, and c-Myc and Statmin1 oncogenes in the development of HCC and HCA in sub-chronically exposed mouse liver. All of these are thoroughly explored below.

Furan metabolism to BDA by CYP2E1 has been reported to produce ROS and deplete GSH (Carfagna et al., 1993; Gonzalez, 2005; Lu and Cederbaum, 2008). We thus examined our transcriptomic data for evidence of induction of the oxidative stress response. First, we noted that ROS and hydrogen peroxide were predicted upstream regulators. Second, ROS producing CYP2E1 ligands (APAP, CT, and EtOH) were also predicted to be activated (Fig. 3A), which indicates that the downstream effects of their metabolism are likely comparable to those of furan. In addition, meta-analysis of gene expression data shows similarity in furan-induced changes in genes expression to those induced by hydrogen peroxide, CT and APAP (Fig. 2; Suppl. Table S6). Third, Enrichment of the Glutathione Metabolism pathway (Fig. 3A) and increased expression of relevant genes including glutamate-cysteine ligase, catalytic subunit (Gclc) and glutathione-s-transferases (Gsta1 was the most highly up-regulated gene; Table 2), all indicated that ROS production, and GSH synthesis and depletion are occurring. Finally, the most important indication of oxidative stress was the enrichment of the Nrf2 Oxidative Stress Response pathway and prediction of Nrf2 activation (Fig. 3). ROS can directly modify protein activity by oxidizing disulfhydryl (R-(SH)2) groups on redox-sensitive proteins to disulfides (R-S2). A canonical example of this is the oxidation of Kelch-like ECH protein 1 (KEAP1), which can only bind to and block Nrf2 in its reduced form. Reduced KEAP1 holds Nrf2 in the cytoplasm with the Cullin-3 (CUL3) ubiquitin ligase, which targets Nrf2 for proteosomal degradation. When KEAP1 is oxidized, Nrf2 enters the nucleus and alters gene expression of anti-oxidant and metabolic genes. We see up-regulation of a plethora of downstream genes committed to redox homeostasis and xenobiotic metabolism (Table 2). Therefore our toxicogenomics data clearly support that hepatic, ROS-associated toxicity occurs in the mouse model following exposure to furan.

Clearing ROS from the cell is a top priority as, in addition to modifying protein function by redox, ROS can damage cellular macromolecules (such as DNA). This damage, when severe, leads to cellular necrosis. Two cellular systems for neutralizing ROS are thioredoxin and peroxiredoxin, which are involved in redox repair of oxidized cellular

Non-Hepatocarcinogens Щ Hepatocarcinogens Щ Other (cancer, regeneration, inflammation)! I

0 0 I_

Cyp2E1 ligands

High dose furan and liver cancer

Low dose furan and Inflammation

Hepatocarcinogens

Non Hepato-carcinogens

APAP (300mg/kg)_ 12 hr_G5 C T( 1 ml/kg12 hr_G5 APAP (300moftg)_ 24 hr^G5 BaP(150mg/kg|_24hr_G7 BaP(300mj/VgL24hr_G7 PG8E(1200mg/Xg)_13wk_G1 MECL(2000mg/kg)_13wk_G1 LR_24f*_G4 TC PN( 40mgfVg)_ 13wk_G 1 EtOx (100mg/kg)_ 13wk_G 1 VANP(2mgikg)_ 13wk_G1 CT(1m№g)_6hr„G5 APAP(300mg/kg)_6hf_G5 LR_38hr_G4

LR_48hr_G4 Fu ran(8 mg/kg |_3wk_G6 liver cancer_G2 Livef hliainiiBtJon_G3 Fu ran(2 rn^Vg (_3wk_G6 Fu ran(1 in(y'k,g |_3wk_G6 Fu ran(4 mg/kg t_3wk_G6 BaP(5mg/kg)_24hf_G7 Ba P(50mg/kg |_24 hr_G7 MAC R( 50mg/Vgi_ 13wk_G 1 TFEl(1250mg/kg)_13wk_G1 TFE A(50000mg/kg)_ 13wk_G 1 DCBZ (300mgJVg)_ 13wk_G 1 TDPP( 1000mg/kg)_ 13wk_G 1 AD6Q(20Q00mg/Kg)_ 13wk_G 1 CT(1 ml/kg >_24hr_G5 TCFM(3925mg/kg)_ 13wk_G1 IOOO(94mg/kg)_ 13wk_G1 COUM (200mg/kg)_ 13wk_G 1 EtDA(2000mg/kg)_ 13wk_G1 8ZFU(2000mg/kg)_ 13wk_G1 PCnB(2000mg/kg)_ 13wk_G1 MAL A(16000mg/kg)_ 13wk_G1 N AAC (10000mg/kg)_ 13wk_G1 DIAZ(200mg/kg)_ 13wk_G1 BB M P< 1250mg/kg)_ 13wk_G 1 CMPH (250mg/kg)_ 13wk_G1 BENZ( 100mg/kg)_ 13wk_G1 D6E T( 62mg/kg)_ 13wk_G1

0 2 _I_

04 -j_

06 _I—

0 8 _I_

Fig. 2. Cluster analysis comparing furan gene expression data with publically available datasets that reflect liver cancer (HCC), liver regeneration (LR), Cyp2E1 ligand-, and hepatocarcinogen- and non-hepatocarcinogen-induced changes in gene expression in mouse liver. Datasets were downloaded from GEO, where G1 is GSE18858, G2 is GSE26538, G3 is GSE35934, G4 is GSE20427, G5 is GSE4874, G6 is this study (GSE48644), and G7 is from another study conducted by our laboratory (Yauk et al., unpublished data; manuscript in preparation). Additional information for individual studies can be found in Suppl. Table S3.

proteins and detoxification of hydrogen peroxide, respectively. NF-kB and AP-1 are both redox-sensitive proteins that require reduction by thioredoxin for DNA binding (Arner and Holmgren, 2000). Up-regulation of the NRF2-target-genes thioredoxin reductase (Trxnd1) and sulfiredoxin (Srxn1) in our furan-exposed samples suggest that both of these redox pathways are active (since TRXND1 and SRXN1 reduce oxidized thioredoxin and peroxiredoxin, respectively) (Fig. 3B). NRF2 control over maintenance of cellular redox state is normally cytoprotective; however, chronic elevation of NRF2-mediated gene expression has been shown to be detrimental and even to confer a survival advantage to pre-cancer and cancer cells (Hayes and McMahon, 2006;

Lau et al., 2008). We propose that this dual nature of NRF2 activation might be a tipping point between normal cytotoxicity/regeneration cycles and dysregulated ones. Moreover, chronic exposure to a ROS-producing agent, such as furan, could tip the balance by causing chronic KEAP1-inactivation/NRF2-activation and sustained maintenance of downstream effects (Fig. 4).

Cytotoxicity and liver injury following chemical exposure occur first in furan's proposed MOA. Our transcriptomic data indicates that cellular damage is occurring in our furan-exposed liver tissues. Meta-analysis of gene expression produced similar changes as the following disease states: Injury to Liver, Hepatic Fibrosis and Cirrhosis of the Liver

Growth factors

MAPK/ERK Signaling

Caspase cascade

Apoptosis

NF-KB .

Pro-inflammatory genes

- Anti-apoptotic genes

Bax C-Jlin/AP-l Proliferation genes

OÖO0^

Apoptotic genes (oxidized cell) Proliferation genes (reduced cell

Ras/Raf

MEKl/2

I Oxidation

ERKl/2 ( J DUSP6

Substrate (c-Jun/AP-1,

JJ, c-Myc, Stathminl...)

Proliferation

Fig. 3. Furan MOA. (A) Summary of predicted upstream regulators (orange), enriched molecular pathways (green) with BMDL values indicated, and selected differentially expressed genes (bold text) for the 8 mg/kg bw dose group. (B) A more detailed image of the proposed redox switch and other molecular signals implicated in the furan MOA.

Margin of exposure (MoE) values for furan-induced molecular pathways and cancer endpoints. A MoE threshold of 1000 has been previously used for furan (IRIS, 2012). BMDLs for HCC and HCA were modeled using previously published data from Moser et al. (2009) and NTP (1993) for female B6C3F1 mice. BMDL mean for pathway data were calculated using transcriptional data from this study in BMDExpress. Estimated human exposure levels were taken from Carthew et al. (2010).

NRF2 oxidative stress response pathway All-pathway-averagea Hepatocellular carcinoma Hepatocellular adenoma

(This study) (This study) Moser et al. (2009) NTP (1993) Moser et al. (2009) NTP (1993)

BMDL (|0g/kg bw): 2250 2200 4218 5303 1567 923

Exposure estimates for child and adult (>2 years)

Average (50%): 0.3 |og/l<g-bw/day 7500 7333 14,060 17,677 5223 3077

High level (90%): 0.6 |g/kg-bw/day 3750 3667 7030 8838 2612 1538

Europe

Average (50%): 0.8 |g/kg-bw/day 2813 2750 5273 6629 1959 1154

High level (90%): 1.75 |g/kg-bw/day 1286 1257 2410 3030 895 527

Exposure estimates for infants (0-1 years)

Average (50%): 0.4 |g/kg-bw/day 5625 5500 10,545 13,258 3918 2308

High level (90%): 1 |g/kg-bw/day 2250 2200 4218 5303 1567 923

Europe

Average (50%): 0.3 |g/kg-bw/day 7500 7333 14,060 17,677 5223 3077

High level (90%): 1 |g/kg-bw/day 2250 2200 4218 5303 1567 923

a Calculated using the average BMDL for all enriched pathways (enriched by 8 mg/kg bw furan, with at least 5 differentially expressed genes in each pathway).

(Suppl. Table S5). In addition we see enrichment of GO terms for Apo-ptosis and Regulation of Programmed Cell Death (Suppl. Table S7). Enrichment of the Arachidonic Acid Metabolism pathway suggests that Cyp2E1-catalyzed lipid peroxidation of polyunsaturated fatty acids likely contributes to cytotoxicity through apoptosis (Chen et al., 1998). Upstream regulator prediction is a powerful tool that makes predictions regarding molecule activation or inhibition based on downstream gene expression. TNF-alpha was the only predicted upstream regulator in response to 4 mg/kg bw furan and was among many predicted for the 8 mg/kg bw dose group (Fig. 3). As such, we believe that TNF-alpha must be among the earliest cytokines to respond to furan exposure. TNF-alpha is a cytokine that signals through its receptor, TNF receptor 1 (TNFR1). During EtOH metabolism by CYP2E1, oxidative stress promotes cell death by converting the (normally mitogenic) TNF-alpha signal to a cytotoxic one and inducing mitochondrial permeability transition (MPT) (Lu and Cederbaum, 2010; Pastorino and Hoek, 2000). MPT occurs when the MPT pores in the inner mitochondrial membrane open, leading to membrane depolarization and uncoupling that causes ATP depletion and eventually necrosis. MPT pore opening can also result in mitochondrial swelling, leading to mitochondrial rupture, cytochrome c release, caspase activation and apoptosis. MPT is triggered by a variety of stimuli (e.g., various xenobiotics, Bax overexpression, disulfide formation, increased mitochondrial Ca+, GSH depletion, TNF-alpha signaling, and ROS) and the outcome (necrosis versus apoptosis) is dependent on the biological context but is usually a combination of the two (Jaeschke et al., 2002). Further evidence to support the important role of the TNF superfamily (SF) in the cellular response to furan is the enrichment of the Death Receptor (DR) pathway based on up-regulation of Death Receptor 5 (DR5, aka TnfrsflOb; Table 2), which signals for cell death through TNFR-associated death domain (TRADD), Fas-associated death domain (FADD), and the caspase cascade.

SAPK phosphorylation cascades can also signal for MPT (Ashkenazi, 2002; Chu, 2013; Yuan and Kaplowitz, 2009). During APAP metabolism by CYP2E1, ROS production and GSH depletion lead to the activation of apoptosis signaling kinase (ASK1). The ensuing ASK1, MAP2K4, JNK1 phosphorylation cascade leads to activation of Bcl-2-associated X protein (Bax), thereby causing hepatocyte cell death by MPT (Seki et al., 2012). Interestingly, when the cell regains redox homeostasis, reduced thioredoxin can block (ASK1) (Saitoh et al., 1998), and thereby inhibit this cell death pathway (Fig. 3B). Enrichment of 14-3-3 signaling (which has a ASK1, JNK1 component), prediction of activation of the MAPK group and MAP2K4 (Fig. 3), enrichment of GO terms for protein kinase activity and regulation of protein phosphorylation (Suppl.

Table S7), and up-regulation of Bax and Dusp8 (a JNK-specific phosphatase) (Table 2, Fig. 3) all indicate that cytotoxicity from increased MPT (as per Seki et al. (2012)) is relevant to furan. Another JNK1 substrate, c-Jun, is a monomer of the heterodimeric transcription factor AP-1. AP-1 has been described as a 'double-edged sword' in tumorigenesis because, depending on dimer composition and abundance, and the oxidation state of the cell, it has been shown to have either oncogenic or tumor suppressor capabilities (Eferl and Wagner, 2003; Hess et al., 2004). c-Jun's contribution is typically understood to be pro-proliferative; however, during periods of CYP2E1 activity and oxidative stress, c-Jun and JNK1 are also known to be involved in hepatocyte cell death (Amir et al., 2012; Liu et al., 2002; Singh and Czaja, 2007; Singh et al., 2009). Indeed, loss of c-Jun has been shown to prevent the CYP2E1-induced TNF-alpha switch from proliferative to apoptotic signaling (Liu et al., 2002). The exact mechanisms by which c-Jun functions during periods of oxidative stress compared to normal cellular conditions are under investigation (Amir et al., 2012; Karin and Shaulian, 2001). The combination of Cyp2E1 activation, ROS production and up-regulation of c-Jun and Dusp8 suggests that the oxidative stress- and JNK1/c-Jun-mediated cell death mechanism (as per Amir et al. (2012)) is likely to also be relevant to furan.

Inflammation accompanies cell death and is triggered by proinflammatory cytokines, released by Kupffer cells (liver-specific macrophages), and molecular signals released from dying cells (Nakagawa and Maeda, 2012). Although acute inflammation is a normal part of the innate immune response, chronic inflammation is thought to link the cytotoxicity/regeneration MOA with HCC. This link is thought to be the result of an increase in pro-proliferative cytokines, pro-angiogenic factors and transcription factors with pro-proliferative target genes, particularly NF-kB (Berasain et al., 2009; DiDonato et al., 2012; Karin, 2009). TNF-alpha can signal to activate NF-kB, which is important for initiation of both inflammation and hepatocyte proliferation (Karin, 2006; Nakagawa and Maeda, 2012). NF-kB is a dimeric transcription factor that is typically retained in the cytoplasm by its inhibitor (IkB); however, up-regulation of IkB kinase (iKBk, which has three sub-units: a, (3, y) allows NF-kB to enter the nucleus and modify gene expression. NF-kB activation in Kupffer cells leads to production of proinflammatory cytokines, growth and survival signals, and angiogenic factors (Karin, 2006). Kupffer cell-specific NF-kB inactivation (by iKBk-Y knock-out) has been shown to reduce HCC tumor load, presumably by reducing production of inflammatory signals (Berasain et al., 2009; Maeda et al., 2005). NF-kB activation in hepatocytes leads to up-regulation of cell-cycle and anti-apoptotic genes (Karin, 2006) and its inhibition prevents tumor development and triggers apoptosis of

Adaptive outcome (Regeneration) __Cytotoxicity signalling

■=> I

Furan < 2.0 mg/kg bw

Cell death

NRF2 oxidative stress response

Inflammation (NF-kB)

' Acute_Chronlc_ ^

Proliferative signaling

Regeneration

Healthy liver

Adverse outcome (Cancer)

Furan > 2.0 mg/kg bw

Cytotoxicity signalling (TNF, ROS, SAPK, DR5...)

Cell death

BMDL: 2.25 mg/kg bw

Proliferative signaling (ERK, TNF, NF-kB, c-Jun,

oncogenes...)

Regeneration

Liver cancer

Fig. 4. Schematics for the adaptive (upper panel) and adverse (lower panel) outcomes in response to non-carcinogenic (upper) and carcinogenic (lower) doses of furan. We propose that NRF2 activation and inflammation could represent important tipping points in the MOA and that, when chronically activated, drive the cancer outcome.

initiated hepatocytes (Pikarsky et al., 2004). In the presence of ROS, NF-kB and JNK1 are known to antagonize one-another (Nakano et al., 2006; Papa et al., 2004; Zhang and Chen, 2004). In addition, NF-kB itself is a redox sensitive protein (Yuan and Kaplowitz, 2009). While the details of this molecular cross-talk are still being investigated, it seems reasonable to assume that if NF-kB activity were to get the 'upper-hand' it could tip the balance away from JNK1-mediated cytotoxicity and toward NF-KB-mediated inflammation and proliferation.

We expect that NF-kB is transcriptionally active in response to furan since we see up-regulation of¡hBk-y (Table 2) and the NF-kB monomer, NF-KB1, is predicted to be activated (Fig. 3). In addition, our 1, 2 and 4 mg/kg dose groups each clustered most closely with an NF-KB-mediated liver inflammation dataset (GSE20427; Fig. 2) and our 8 mg/kg bw dose group is correlated with the 'Inflammatory Disease of the Liver and Inflammatory Disorder' phenotypes (Suppl. Table S5). These data are consistent with the notion that NF-KB-mediated inflammation in pre-cancer cells is an important part of setting the stage for malignant transformation in furan exposed liver (but may not be sufficient for transformation on its own since 1 and 2 mg/kg bw doses furan do not cause HCC). We observed up-regulation of the pro-inflammatory chemokine C-X-C motif ligand (Cxcl10; Table 2), and prediction of activation of other cytokines (IL-1a, IL-1b, IL-6, IL-8, TNFSF11, TNF-alpha, CXCL12, CCL5, IFN-y) and of Triggered Receptor Expressed on Myeloid Cells 1 (TREM1) (Fig. 3A). Kupffer cell expression of TREM1 has been shown to be required for HCC development; further, loss of TREM1 results in down-regulation of transcription of pro-inflammatory and pro-proliferative cytokines (IL-6, IL-b1, TNF-alpha, CCL2 and CXCL10) and loss of activation of ERk1/2,JNK and NF-kB, ultimately preventing malignant transformation (Wu et al., 2012). Our data support the idea that inflammation precedes and facilitates carcinogenic transformation in chemically induced HCC. These findings are consistent with previous studies that have demonstrated the chronic inflammation phenotype in furan-exposed mice and rats (Moser et al., 2009; NTP, 1993).

Regenerative proliferation is the penultimate stage to tumor formation in the furan MOA. Molecular mechanisms that underlie hepatocyte proliferation have been studied using liver development, liver regeneration and spontaneous HCC models. NF-kB and c-Jun are both required

for liver development. Loss of NF-kB (by knock-out of lKBk-p> or ¡KEk-y) is embryonic lethal because of extensive hepatocyte apoptosis and liver degeneration (Li et al., 1999; Rudolph et al., 2000). c-Jun-null mice are also embryonic lethal as a result of impaired liver and heart development (Behrens et al., 2002; Eferl et al., 1999; Jochum et al., 2001). Liver regeneration has been extensively studied in animals following two-thirds partial hepatectomy (PH). Both NF-kB and c-Jun are not normally expressed/activated in adult liver, but are up-regulated following PH in a TNF-alpha/TNFR1-dependent manner, and are required for hepatic regeneration (Alcorn et al., 1990; Chaisson et al., 2002; Cressman et al., 1994; Diehl et al., 1994; FitzGerald et al., 1995; Hilberg et al., 1993; Yamada et al., 1997). Rats lacking NF-kB display impaired liver regeneration due to hepatocyte apoptosis (Schrum et al., 2000) and mice lacking AP-1 also show impaired liver regeneration, often leading to death (Behrens et al., 2002). Similarly, in HCC TNF-alpha signals activation of both c-Jun/AP-1 and NF-kB, which up-regulate expression of pro-survival and pro-proliferation genes. Increased NF-kB expression in HCC cells is important for the transition from liver injury and inflammation to cancer (Karin, 2009; Luedde and Schwabe, 2011) and for HCC cell growth (Pikarsky et al., 2004; Wang et al., 2003). Recently, up-regulation of c-Jun gene expression in the early stages of liver cancer has been shown to be required for HCC initiation and development (Min et al., 2012). Both c-Jun and NF-kB have been proposed as molecular targets for liver tumor prevention (Young et al., 2003).

Meta-analysis of our transcriptomic data demonstrates similarity in gene expression changes produced by 8 mg/kg bw furan to those generated during liver regeneration, hepatocellular dysplasia, and liver cancer (Suppl. Table S5). Our 8 mg/kg bw dose group clustered most closely with the spontaneous HCC and liver regeneration datasets (Fig. 2). This is remarkable since our exposure was sub-chronic (therefore cells had presumably not yet undergone malignant transformation). We observed up-regulation of growth factors (Pdgf and Gdf9/15; Table 2) and prediction of activation of others (VEGF, HGF, IGF1, EGF, FGF, BMP2, BMP6; Fig. 3A) in response to furan. Growth factors are secreted proteins that bind their tyrosine receptor kinase and relay signals through the cell via extracellular signal-regulated kinases (ERKs). The MAPK phosphorylation cascade follows the order of RAS (GTPase),

RAF (MAP3K), MEK1/2 (MAP2K), and ERK1/2 (MAPK), all of which are predicted upstream activators as is the MAPK/ERK Signaling pathway (Fig. 3). ERKs were first identified on the basis of being required for the G1/S transition (Meloche and Pouyssegur, 2007) and are known to activate at least 160 substrates (Min et al., 2011), which go on to alter target gene expression. One target gene is dual-specific phosphatase 6 (Dusp6), which is up-regulated in response to furan (Table 2) and feeds back to reset the pathway by dephosphorylating ERK1/2. During periods of chronic oxidative stress, the ERK1/2 pathway is protective against hepatocyte death (Singh and Czaja, 2007); this could be due to the fact that DUSPs are oxidation-inactivated (Patterson et al., 2009). We suggest that ERK signaling is responsible for promoting hepatic cellular proliferation during liver regeneration.

In our samples, ERK signaling might promote cellular proliferation via activation of protein substrates including: c-Jun/AP-1, c-Myc, and Stathmin1. c-Jun has been shown to promote HCC by preventing hepatocyte apoptosis through inhibition of p53 tumor suppressor (Eferl et al., 2003) and SIRT6 (Sirtuin-6; a stress-responsive deacetylase). For the latter, c-Jun blocks up-regulation of Sirt6, which results in increased levels of Survivin (aka: Birc5) thereby blocking apoptosis (Min et al., 2012). Survivin is marginally up-regulated by 8 mg/kg bw furan. Two other AP-1 monomers were up-regulated by furan: activating transcription factors 3 and 4 (Atf3/4; Table 2). ATF3 is up-regulated during the stress response (Hai et al., 1999) and has been implicated in oncogenesis when bound to c-Jun (Thompson et al., 2009). ATF4 is up-regulated in cancer and is known to regulate genes involved in oxidative stress response, differentiation, metastasis and angiogenesis (Ameri and Harris, 2008). C-Myc oncogene is a transcription factor that is associated with many cancers (Dang, 2012) and has been shown to contribute to HCC progression during chronic EtOH exposure and CYP2E1 activation (Tsuchishima et al., 2013) and to up-regulate Nrf2 expression during tu-morigenesis (Denicola et al., 2011). Finally, the Breast Cancer Regulation by Stathminl pathway is enriched in both the 4 and 8 mg/kg bw furan dose groups (Fig. 3A). Stathmin1 (Stmn1; aka Oncoprotein18, Op18) is a microtubule remodeling protein that is not actually specific to breast tissue/breast cancer (Stmnl and Tubs are up-regulated by furan; Table 2). Indeed, it has been shown to be involved in both liver regeneration (Okazaki et al., 1993; Rowlands et al., 1995) and HCC. Up-regulation of Stathmin1 in HCC is correlated with poor prognosis and tumor reoccurrence (Hsieh et al., 2010; Yuan et al., 2006). Silencing of Stathminl in an HCC cell line (HCCLM3) suppressed proliferation, induced apoptosis and impaired metastasis (Gan etal., 2010). ERKs likely direct hepatic cellular proliferation via activation of other substrates as well; however we believe that c-Jun, c-Myc, and Stathmin1 likely play important roles during this process.

Taken together, we propose that all of the above molecules and pathways (summarized in Fig. 3) are important contributors to furan-induced regenerative proliferation leading to HCC. We also acknowledge that, in order for HCA and HCC to develop, genomic damage must occur. We propose that furan-induced genomic damage occurs indirectly as a by-product of its metabolism (particularly the oxidative stress and cytotoxicity aspects of the MOA). This hypothesis is based on our observed induction of oxidative stress responses in combination with the mixed genotoxicity findings (discussed in the introduction; occurring at high doses). In support of the occurrence of DNA damage at high doses, we saw enrichment of the p53 Signaling pathway in response to 8 mg/kg bw furan (but not 4 mg/kg bw). This pathway enrichment was due to the up-regulation of the following genes: Cdknla, Cyclin Gl, Bax, DR5, and Teapl (p53 inducible nuclear protein 1), which are all downstream targets of p53 that are suggestive of DNA damage response.

In addition to applying our transcriptomics data to derive MOA information for furan, another objective of this work was to demonstrate the practical application of toxicogenomics in quantitative risk assessment. We accomplished this by calculating BMD and MoE values for apical and transcriptional endpoints since demonstration of a strong

correlation between transcriptomic and apical endpoints will facilitate the implementation of toxicogenomics in routine risk assessments (Thomas et al., 2011). Strikingly, the majority of the transcriptional BMD modeling that was conducted placed most values near to the expected 2-4 mg/kg bw POD range (Table 2, Fig. 3A, Suppl. Tables S3 and S7) as did modeling of HCC and HCA cancer endpoints reported by Moser etal. (2009) and the NTP (1993) (Table 1). It is worth noting that, while the transcriptional BMD/BMDL values were often modeled using genes that had fold changes of less than ± 2.0, the most significantly enriched pathways and functions tended to also contain genes with larger fold changes. In addition, the magnitude of the fold change of the individual probes did not appear to impact the magnitude of the BMD/BMDL values (Table 2). Modeling of previously-studied cytotoxicity and proliferation apical endpoints showed more variability. BMD values calculated from apical data for male and female B6C3F1 mice published by Gill et al. (2011) were often below the 2 mg/kg bw cancer threshold and resembled their NOAEL value (0.12 mg/kg bw) for furan. Alternatively, BMD values calculated using apical data published in the NTP's (1993) rodent cancer bioassay were much higher than the 2-4 mg/kg bw range. The most sensitive cytotoxicity apical endpoint was hepatocyte apoptosis in the caudate lode (BMDL: 0.11 mg/kg bw), and the most sensitive proliferation apical endpoint was hepatocyte ba-sophilia (BMDL: 0.30 mg/kg bw) (Table 1). Both of these underestimate the known cancer POD. The BMDLs for liver cancer more closely approximated cancer POD with the HCA and HCC BMDLs equal to 0.92 and 1.57 mg/kg bw, respectively (Table 1). Transcriptional (pathway) BMDLs were within the same range of the cancer BMDLs (1.50-4.40 mg/kg bw), with an average pathway BMDL of 2.20 mg/kg bw. The two pathways that we believe are very important to the furan MOA, NRF2 Oxidative Stress Response and ERK/ MAPK Signaling, have very similar mean BMDLs of 2.25 and 2.62 mg/kg bw, respectively (Fig. 3). Strikingly, after only a three week exposure, both the individual and the all-pathway-average BMDL values perfectly reflect the 2-4 mg/kg bw cancer POD range.

We used the cancer and pathway BMDL values to estimate human health risk using MoE values. Pathway MoE values were calculated using either the NRF2 Oxidative Stress Response Pathway BMDL, or the all-pathway-average BMDL. Pathway and HCA MoEs calculated for transcriptional or apical endpoints were generally equivalent (Table 3). While it is still early days and this work should be repeated using other chemicals, we feel that this concordance between individual and all-pathway-average BMDL (and therefore MoE) is an interesting and important observation and tentatively suggest that an all-pathway-average MoE might be a useful screening tool for initial estimation of chemical risk.

In summary, using the MOA and the risk assessment information gathered in this study we assembled a simplified adverse outcome pathway for furan (Fig. 4). We believe that chronic exposure to > 2 mg/kg bw furan hijacks the normal response to manageable levels of oxidative stress-induced damage thereby tipping the balances in favor of chronic regenerative proliferation. The molecular initiating event is Cyp2E1 ligand binding and furan metabolism to produce ROS. The key events that follow are: (1) cytotoxicity, (2) activation of the NRF2 oxidative stress response, (3) inflammation, and (4) regenerative proliferation, where (2) and (3) represent tipping points between the adaptive and the adverse outcomes. There is a large body of literature supporting the idea that chronic NRF2 activation is carcinogenic, and it is known that the POD for furan is between 2 and 4 mg/kg bw. Therefore, we propose that the NRF2 oxidative stress response pathway BMDL (2.25 mg/kg bw) would denote the best transcriptional benchmark dose for the POD between the adaptive and adverse responses to furan. Above this dose the cytoprotective and defensive roles of NRF2 are hijacked in premalignant and malignant cells to confer a survival advantage, and this is exacerbated by chronic inflammation and up-regulation of other oncogenes (including c-Jun, c-Myc and

Stathminl). While our molecular MOA should be validated by additional experiments, we feel that it provides realistic first insight into how furan causes cancer. We believe that our toxicogenomic approach is illustrative of how gene expression data can provide MOA information for indirectly acting and non-genotoxic carcinogens that is reflective of apical cancer data. In addition, we believe that the consistency between apical and genomic data demonstrates the utility of this approach for the investigation of other non-genotoxic carcinogens that do not have apical or two-year rodent bioassay data associated with them.

Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.taap.2013.10.019.

Conflict of interest

The authors declare no conflicts of interest.

Acknowledgments

We thank Byron Kuo for completing gene and pathway BMD analyses in BMDExpress. We thank the ILS animal care group for performing the furan mouse exposures. This research was funded by the Health Canada Genomics R&D Initiative. A.F.J. was supported by the Natural Science and Engineering Research Council of Canada.

References

Alcorn, J.A., Feitelberg, S.P., Brenner, D.A., 1990. Transient induction of c-Jun during hepatic regeneration. Hepatology 11,909-915.

Ameri, K., Harris, A.L., 2008. Activating transcription factor 4. Int. J. Biochem. Cell Biol. 40, 14-21.

Amir, M., Liu, K., Zhao, E., Czaja, M.J., 2012. Distinct functions of JNK and c-Jun in oxidant-induced hepatocyte death. J. Cell. Biochem. 113, 3254-3265.

Arner, E.S.J., Holmgren, A., 2000. Physiological functions of thioredoxin and thioredoxin reductase. Eur. J. Biochem. 267, 6102-6109.

Ashkenazi, A., 2002. Targeting death and decoy receptors of the tumour-necrosis factor superfamily. Nat. Rev. Cancer 2,420-430.

Becker, RA, Chambers, J.M., Wilks, A.R., 1998. The New S Language. Chapman and Hall/ CRC.

Behrens, A., Sibilia, M., David, J., Möhle-Steinlein, U., Tronche, F., Schütz, G., Wagner, E.F., 2002. Impaired postnatal hepatocyte proliferation and liver regeneration in mice lacking c-Jun in the liver. EMBO J. 21,1782-1790.

Benford, D., Bolger, P.M., Carthew, P., Coulet, M., DiNovi, M., Leblanc, J.C., Renwick AG., Setzer, W., Schlatter, J., Smith, B., Slob, W., Williams, G., Wildemann, T., 2010. Application of the Margin of Exposure (MOE) approach to substances in food that are genotoxic and carcinogenic. Food Chem. Toxicol. 48 (Suppl. 1), S2-S24.

Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R Stat. Soc. Series B (Methodological) 57, 289-300.

Berasain, C., Castillo, J., Perugorria, M.J., Latasa, M.U., Prieto, J., Avila, M.A., 2009. Inflammation and liver cancer: new molecular links. Ann. N. Y. Acad. Sci. 1155,206-221.

Black M.B., Budinsky, RA, Dombkowski, A., Cukovic, D., LeCluyse, E.L., Ferguson, S.S., Thomas, R.S., Rowlands, J.C., 2012. Cross-species comparisons of transcriptomic alterations in human and rat primary hepatocytes exposed to 2,3,7,8-tetrachlorodibenzo-p-dioxin. Toxicol. Sci. 127,199-215.

Bolstad, B.M., Irizarry, RA, Astrand, M., Speed, T.P., 2003. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19,185-193.

Carfagna, M.A., Held, S.D., Kedderis, G.L., 1993. Furan-induced cytolethality in isolated rat hepatocytes: correspondence with in vivo dosimetry. Toxicol. Appl. Pharmacol. 123, 265-273.

Carthew, P., DiNovi, M., Woodrow Setzer, R., 2010. Application of the Margin of Exposure (MOE) approach to substances in food that are genotoxic and carcinogenic: example: CAS no: 105650-23-5 PhIP (2-amino-1-methyl-6-phenylimidazo[4,5-b]pyridine). Food Chem. Toxicol. 48 (Suppl. 1), S98-S105.

Chaisson, M.L., Brooling, J.T., Ladiges, W., Tsai, S., Fausto, N., 2002. Hepatocyte-specific inhibition of NF-kB leads to apoptosis after TNF treatment, but not after partial hepatec-tomy. J. Clin. Invest. 110,193-202.

Chen, Q., Galleano, M., Cederbaum, A.I., 1998. Cytotoxicity and apoptosis produced by ar-achidonic acid in HepG2 cells overexpressing human cytochrome P-4502E1. Alcohol. Clin. Exp. Res. 22, 782-784.

Chu, W.M., 2013. Tumor necrosis factor. Cancer Lett. 328, 222-225.

Cordelli, E., Leopardi, P., Villani, P., Marcon, F., Macri, C., Caiola, S., Siniscalchi, E., Conti, L., Eleuteri, P., Malchiodi-Albedi, F., Crebelli, R., 2010. Toxic and genotoxic effects of oral administration of furan in mouse liver. Mutagenesis 25,305-314.

Cressman, D.E., Greenbaum, L.E., Haber, B.A., Taub, R., 1994. Rapid activation of post-hepatectomy factor/nuclear factor kB in hepatocytes, a primary response in the regenerating liver. J. Biol. Chem. 269,30429-30435.

Cui, X., Hwang, J.T., Qiu, J., Blades, N.J., Churchill, GA, 2005. Improved statistical tests for differential gene expression by shrinking variance components estimates. Biostatis-tics 6, 59-75.

Dang, C.V., 2012. MYC on the path to cancer. Cell 149,22-35.

Davis, JA, Gift, J.S., Zhao, QJ., 2011. Introduction to benchmark dose methods and U.S. EPA's benchmark dose software (BMDS) version 2.1.1. Toxicol. Appl. Pharmacol. 254,181-191.

Denicola, G.M., Karreth, F.A., Humpton, T.J., Gopinathan, A., Wei, C., Frese, K., Mangal,

D., Yu, K.H., Yeo, C.J., Calhoun, E.S., Scrimieri, F., Winter, J.M., Hruban, R.H., lacobuzio-Donahue, C., Kern, S.E., Blair, I.A., Tuveson, D.A., 2011. Oncogene-induced Nrf2 transcription promotes ROS detoxification and tumorigenesis. Nature 475,106-110.

DiDonato, J.A., Mercurio, F., Karin, M., 2012. NF-kappaB and the link between inflammation and cancer. lmmunol. Rev. 246, 379-400.

Diehl, A.M., Yin, M., Fleckenstein, J., Yang, S.Q., Lin, H.Z., Brenner, DA, Westwick J., Bagby, G., Nelson, S., 1994. Tumor necrosis factor-alpha induces c-Jun during the regenerative response to liver injury. Am. J. Physiol. 267, G552-G561.

Ding, W., Petibone, D.M., Latendresse, J.R., Pearce, M.G., Muskhelishvili, L., White, G.A., Chang, C., Mittelstaedt, RA, Shaddock J.G., McDaniel, L.P., Doerge, D.R., Morris, S.M., Bishop, M.E., Manjanatha, M.G., Aidoo, A., Heflich, R.H., 2012. In vivo genotoxicity of furan in F344 rats at cancer bioassay doses. Toxicol. Appl. Pharmacol. 261,164-171.

Dunnett, C.W., 1955. A multiple comparison procedure for comparing several treatments with a control. J. Am. Stat. Assoc. 50 (1096-1096-1121).

Durling, L.J., Svensson, K., Abramsson-Zetterberg, L., 2007. Furan is not genotoxic in the micronucleus assay in vivo or in vitro. Toxicol. Lett. 169,43-50.

Eferl, R., Ricci, R., Kenner, L., Zenz, R., David, J., Rath, M., Wagner, E.F., 2003. Liver tumor development: c-Jun antagonizes the proapoptotic activity of p53. Cell 112,181-192.

Eferl, R., Sibilia, M., Hilberg, F., Fuchsbichler, A., Kufferath, I., Guertl, B., Zenz, R., Wagner,

E.F., Zatloukal, K., 1999. Functions of c-Jun in liver and heart development. J. Cell Biol. 145,1049-1061.

Eferl, R., Wagner, E.F., 2003. AP-1: a double-edged sword in tumorigenesis. Nat. Rev. Cancer 3, 859-868.

EFSA 2005. EFSA — Scientific Opinion of the Scientific Committee: Opinion of the Scientific Committee on a Request from EFSA Related to a Harmonised Approach for Risk Assessment of Substances Which Are Both Genotoxic and Carcinogenic (2013).

FitzGerald, M.J., Webber, E.M., Donovan, J.R., Fausto, N., 1995. Rapid DNA binding by nuclear factor kB in hepatocytes at the start of liver regeneration. Cell Growth Differ. 6,417-427.

Fransson-Steen, R., Goldsworthy, T.L., Kedderis, G.L., Maronpot, RR, 1997. Furan-induced liver cell proliferation and apoptosis in female B6C3F1 mice. Toxicology 118,195-204.

Gan, L., Guo, K., Li, Y., Kang, X., Sun, L., Shu, H., Liu, Y., 2010. Up-regulated expression of stathmin may be associated with hepatocarcinogenesis. Oncol. Rep. 23,1037-1043.

Gates, L.A., Lu, D., Peterson, L.A., 2012. Trapping of cis-2-butene-1,4-dial to measure furan metabolism in human liver microsomes by cytochrome P450 enzymes. Drug Metab. Dispos. 40, 596-601.

Gautier, L., Cope, L., Bolstad, B.M., Irizarry, RA, 2004. Affy — analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 20, 307-315.

Gill, S., Bondy, G., Lefebvre, D.E., Becalski, A., Kavanagh, M., Hou, Y., Turcotte, A.M., Barker, M., Weld, M., Vavasour, E., Cooke, G.M., 2010. Subchronic oral toxicity study of furan in Fischer-344 rats. Toxicol. Pathol. 38, 619-630.

Gill, S., Kavanagh, M., Barker, M., Weld, M., Vavasour, E., Hou, Y., Cooke, G.M., 2011. Sub-chronic oral toxicity study of furan in B6C3F1 mice. Toxicol. Pathol. 39, 787-794.

Gonzalez, F.J., 2007. The 2006 Bernard B. Brodie award lecture. Cyp2e1. Drug Metab. Dispos. 35,1-8.

Gonzalez, F.J., 2005. Role of cytochromes P450 in chemical toxicity and oxidative stress: studies with CYP2E1. Mutat Res. Fundam. Mol. Mech. Mutagen. 569,101-110.

Hai, T., Wolfgang, C.D., Marsee, D.K., Allen, A.E., Sivaprasad, U., 1999. ATF3 and stress responses. Gene Expr. 7,321-335.

Haseman, J.K., Hailey, J.R., Morris, R.W., 1998. Spontaneous neoplasm incidences in Fischer 344 rats and B6C3F1 mice in two-year carcinogenicity studies: a National Toxicology Program update. Toxicol. Pathol. 26, 428. http://dx.doi.org/10.1177/ 019262339802600318.

Hayes, J.D., McMahon, M., 2006. The double-edged sword of Nrf2: subversion of redox homeostasis during the evolution of cancer. Mol. Cell 21, 732-734.

Hess, J., Angel, P., Schorpp-Kistner, M., 2004. AP-1 subunits: quarrel and harmony among siblings. J. Cell Sci. 117, 5965-5973.

Hickling, K.C., Hitchcock J.M., Oreffo, V., Mally, A., Hammond, T.G., Evans, J.G., Chipman, J.K., 2010. Evidence of oxidative stress and associated DNA damage, increased prolif-erative drive, and altered gene expression in rat liver produced by the cholangiocarcinogenic agent furan. Toxicol. Pathol. 38, 230-243.

Hilberg, F., Aguzzi, A., Howells, N., Wagner, E.F., 1993. c-Jun is essential for normal mouse development and hepatogenesis. Nature 365,179-181.

Hothorn, T., Bretz, F., Westfall, P., 2008. Simultaneous inference in general parametric models. Biom. J. 50, 346-363.

Hsieh, S.Y., Huang, S.F., Yu, M.C., Yeh, T.S., Chen, T.C., Lin, Y.J., Chang, C.J., Sung, C.M., Lee, Y.L., Hsu, C.Y., 2010. Stathmin1 overexpression associated with polyploidy, tumor-cell invasion, early recurrence, and poor prognosis in human hepatoma. Mol. Carcinog. 49, 476-487.

lARC, 1995. Dry cleaning, some chlorinated solvents, and other industrial chemicals. 3407 (63), 3194.

IRIS, 2012. Furan (CASRN 110-00-9). http://www.epa.gov/iris/subst/0056.htm (2013).

Irizarry, RA, Hobbs, B., Collin, F., Beazer-Barclay, Y.D., Antonellis, K.J., Scherf, U., Speed, T.P., 2003. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 4, 249-264.

Jaeschke, H., Gores, G.J., Cederbaum, A.I., Hinson, J.A., Pessayre, D., Lemasters, J.J., 2002. Mechanisms of hepatotoxicity. Toxicol. Sci. 65,166-176.

Jochum, W., Passegué, E., Wagner, E.F., 2001. AP-1 in mouse development and tumorigen-esis. Oncogene 20,2401-2412.

Karin, M., 2009. NF-kappaB as a critical link between inflammation and cancer. Cold Spring Harb. Perspect Biol. 1, a000141.

Karin, M., 2006. Nuclear factor-kappaB in cancer development and progression. Nature 441,431-436.

Karin, M., Shaulian, E., 2001. AP-1: linking hydrogen peroxide and oxidative stress to the control of cell proliferation and death. IUBMB Life 52,17-24.

Kedderis, G.L., Held, S.D., 1996. Prediction of furan pharmacokinetics from hepatocyte studies: comparison of bioactivation and hepatic dosimetry in rats, mice, and humans. Toxicol. Appl. Pharmacol. 140,124-130.

Kellert, M., Brink, A., Richter, I., Schlatter, J., Lutz, W.K., 2008. Tests for genotoxicity and mutagenicity of furan and its metabolite cis-2-butene-1,4-dial in L5178Y tk+/— mouse lymphoma cells. Mutat Res. 657,127-132.

Kupershmidt, I., Su, Q.J., Grewal, A., Sundaresh, S., Halperin, I., Flynn, J., Shekar, M., Wang, H., Park, J., Cui, W., Wall, G.D., Wisotzkey, R., Alag, S., Akhtari, S., Ronaghi, M., 2010. Ontology-based meta-analysis of global collections of high-throughput public data. PLoS One 5. http://dx.doi.org/10.1371/journal.pone.0013066.

Larson, J.L., Wolf, D.C., Butterworth, B.E., 1994. Induced cytolethality and regenerative cell proliferation in the livers and kidneys of male B6C3F1 mice given chloroform by gavage. Fundam. Appl. Toxicol. 23, 537-543.

Lau, A., Villeneuve, N.F., Sun, Z., Wong, P.K., Zhang, D.D., 2008. Dual roles of Nrf2 in cancer. Pharmacol. Res. 58, 262-270.

Leopardi, P., Cordelli, E., Villani, P., Cremona, T.P., Conti, L., De Luca, G., Crebelli, R., 2010. Assessment of in vivo genotoxicity of the rodent carcinogen furan: evaluation of DNA damage and induction of micronuclei in mouse splenocytes. Mutagenesis 25, 57-62.

Li, Q., Van Antwerp, D., Mercurio, F., Lee, K.F., Verma, I.M., 1999. Severe liver degeneration in mice lacking the IkappaB kinase 2 gene. Science 284, 321 -325.

Ligges, U., Mächler, M., 2003. Scatterplot3d — an R package for visualizing multivariate data. J. Stat. Softw. 8, 1 -20.

Liu, H., Jones, B.E., Bradham, C., Czaja, M.J., 2002. Increased cytochrome P-450 2E1 expression sensitizes hepatocytes to c-Jun-mediated cell death from TNF-alpha. Am. J. Physiol. Gastrointest. Liver Physiol. 282, G257-G266.

Lu, Y., Cederbaum, A.I., 2010. CYP2E1 potentiation of LPS and TNFalpha-induced hepato-toxicity by mechanisms involving enhanced oxidative and nitrosative stress, activation of MAP kinases, and mitochondrial dysfunction. Genes Nutr. 5, 149-167.

Lu, Y., Cederbaum, A.I., 2008. CYP2E1 and oxidative liver injury by alcohol. Free Radic. Biol. Med. 44, 723-738.

Luedde, T., Schwabe, R.F., 2011. NF-kB in the liver-linking injury, fibrosis and hepatocellular carcinoma. Nat. Rev. Gastroenterol. Hepatol. 8,108-118.

Maeda, S., Kamata, H., Luo, J.L., Leffert, H., Karin, M., 2005. IKKbeta couples hepatocyte death to cytokine-driven compensatory proliferation that promotes chemical hepatocarcinogenesis. Cell 121, 977-990.

Maga, J.A., Katz, I., 1979. Furans in foods. Crit Rev. Food Sci. Nutr. 11, 355-400.

Manibusan, M.K., Odin, M., Eastmond, D.A., 2007. Postulated carbon tetrachloride mode of action: a review. J. Environ. Sci. Health C Environ. Carcinog. Ecotoxicol. Rev. 25, 185-209.

Mardia, K.V., Kent, J.T., Bibby, J.M., 1979. Multivariate Analysis. Academic Press, London.

McDaniel, L.P., Ding, W., Dobrovolsky, V.N., G., S.J.Jr, Mittelstaedt, RA, Doerge, D.R., Heflich, R.H., 2012. Genotoxicity of furan in Big Blue rats. Mutat. Res. 742, 72-78.

Meloche, S., Pouyssegur, J., 2007. The ERK1/2 mitogen-activated protein kinase pathway as a master regulator of the G1- to S-phase transition. Oncogene 26,3227-3239.

Min, L., He, B., Hui, L., 2011. Mitogen-activated protein kinases in hepatocellular carcinoma development. Semin. Cancer Biol. 21,10-20.

Min, L., Ji, Y., Bakiri, L., Qiu, Z., Cen, J., Chen, X., Chen, L., Scheuch, H., Zheng, H., Qin, L., Zatloukal, K., Hui, L., Wagner, E.F., 2012. Liver cancer initiation is controlled by AP-1 through SIRT6-dependent inhibition of survivin. Nat. Cell Biol. 14,1203-1211.

Moro, S., Chipman, J.K., Wegener, J.W., Hamberger, C., Dekant, W., Mally, A., 2012. Furan in Heat-treated Foods: Formation, Exposure, Toxicity, and Aspects of Risk Assessment. Molecular Nutrition & Food Research.

Moser, G.J., Foley, J., Burnett, M., Goldsworthy, T.L., Maronpot, R., 2009. Furan-induced dose-response relationships for liver cytotoxicity, cell proliferation, and tumorigenic-ity (furan-induced liver tumorigenicity). Exp. Toxicol. Pathol. 61,101-111.

Nakagawa, H., Maeda, S., 2012. Molecular mechanisms of liver injury and hepatocarcinogenesis: focusing on the role of stress-activated MAPK. Pathol. Res. Int. 2012,172894.

Nakano, H., Nakajima, A., Sakon-Komazawa, S., Piao, J.H., Xue, X., Okumura, K., 2006. Reactive oxygen species mediate crosstalk between NF-kappaB and JNK. Cell Death Differ. 13, 730-737.

Neuwirth, C., Mosesso, P., Pepe, G., Fiore, M., Malfatti, M., Turteltaub, K., Dekant, W., Mally, A., 2012. Furan Carcinogenicity: DNA Binding and Genotoxicity of Furan in Rats In Vivo. Molecular Nutrition & Food Research.

NTP, 1993. NTP Technical Report on the Toxicology and Carcinogenesis Studies of Furan in F344/N Rats and B6C3F1 Mice (Gavage Studies). 402.

Okazaki, T., Himi, T., Peterson, C., Mori, N., 1993. Induction of stathmin mRNA during liver regeneration. FEBS Lett. 336,8-12.

Papa, S., Zazzeroni, F., Pham, C.G., Bubici, C., Franzoso, G., 2004. Linking JNK signaling to NF-kappaB: a key to survival. J. Cell Sci. 117,5197-5208.

Pastorino, J.G., Hoek J.B., 2000. Ethanol potentiates tumor necrosis factor-alpha cytotoxic-ity in hepatoma cells and primary rat hepatocytes by promoting induction of the mi-tochondrial permeability transition. Hepatology 31,1141-1152.

Patterson, K.I., Brummer, T., O'Brien, P.M., Daly, R.J., 2009. Dual-specificity phosphatases: critical regulators with diverse cellular targets. Biochem. J. 418,475-489.

Pikarsky, E., Porat, R.M., Stein, I., Abramovitch, R., Amit, S., Kasem, S., Gutkovich-Pyest, E., Urieli-Shoval, S., Galun, E., Ben-Neriah, Y., 2004. NF-kappaB functions as a tumour promoter in inflammation-associated cancer. Nature 431,461-466.

R-Development-Core-Team, 2010. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (http://www.lsw. uni-heidelberg.de/users/christlieb/teaching/UKStaSS10/R-refman.pdf).

Rowlands, D.C., Harrison, R.F., Jones, N.A., Williams, A., Hubscher, S.G., Brown, G., 1995. Stathmin is expressed by the proliferating hepatocytes during liver regeneration. J. Clin. Pathol. Clin. Mol. Pathol. 48, M88-M92.

Rudolph, D., Yeh, W., Wakeham, A., Rudolph, B., Nallainathan, D., Potter, J., Elia, A.J., Mak, T.W., 2000. Severe liver degeneration and lack of NF-kB activation in NEMO/IKK 7-deficient mice. Genes Dev. 14, 854-862.

Saitoh, M., Nishitoh, H., Fujii, M., Takeda, K., Tobiume, K., Sawada, Y., Kawabata, M., Miyazono, K., Ichijo, H., 1998. Mammalian thioredoxin is a direct inhibitor of apoptosis signal-regulating kinase (ASK) 1. EMBOJ. 17, 2596-2606.

Schrum, L.W., Black, D., Iimuro, Y., Rippe, R.A., Brenner, D.A., Behrns, K.E., 2000. c-Jun does not mediate hepatocyte apoptosis following NFkB inhibition and partial hepatecto-my. J. Surg. Res. 88,142-149.

Seki, E., Brenner, D.A., Karin, M., 2012. A liver full of JNK: signaling in regulation of cell function and disease pathogenesis, and clinical approaches. Gastroenterology 143,307-320.

Singh, R., Czaja, M.J., 2007. Regulation of hepatocyte apoptosis by oxidative stress. J. Gastroenterol. Hepatol. 22 (Suppl. 1), S45-S48.

Singh, R., Wang, Y., Schattenberg, J.M., Xiang, Y., Czaja, M.J., 2009. Chronic oxidative stress sensitizes hepatocytes to death from 4-hydroxynonenal by JNK/c-Jun overactivation. Am. J. Physiol. Gastrointest. Liver Physiol. 297, G907-G917.

Templin, M.V., Jamison, K.C., Sprankle, C.S., Wolf, D.C., Wong, B.A., Butterworth, B.E., 1996. Chloroform-induced cytotoxicity and regenerative cell proliferation in the kidneys and liver of BDF1 mice. Cancer Lett. 108, 225-231.

Thomas, R.S., Clewell III, H.J., Allen, B.C., Wesselkamper, S.C., Wang, N.C.Y., Lambert, J.C., Hess-Wilson, J.K., Zhao, QJ., Andersen, M.E., 2011. Application of transcriptional benchmark dose values in quantitative cancer and noncancer risk assessment. Toxicol. Sci. 120,194-205.

Thompson, M.R., Xu, D., Williams, B.R., 2009. ATF3 transcription factor and its emerging roles in immunity and cancer. J. Mol. Med. 87,1053-1060.

Tsuchishima, M., George, J., Shiroeda, H., Arisawa, T., Takegami, T., Tsutsumi, M., 2013. Chronic ingestion of ethanol induces hepatocellular carcinoma in mice without additional hepatic insult. Dig. Dis. Sci. 1 -11.

Venables, W.N., Ripley, B.D., 1999. Modern Applied Statistics with S, 3rd ed. SpringerVerlag, New York.

Wang, J., Huang, Q., Chen, M., 2003. The role of NF-kappaB in hepatocellular carcinoma cell. Chin. Med. J. (Engl) 116, 747-752.

Waters, M.D., Jackson, M., Lea, I., 2010. Characterizing and predicting carcinogenicity and mode of action using conventional and toxicogenomics methods. Mutat. Res. Rev. Mutat. Res. 705,184-200.

Wu, H., Kerr, M., Cui, X., Churchill, G., 2003. MAANOVA: A software package for the analysis of spotted cDNA microarray experiments In the analysis of gene expression data. (http://lectures.molgen.mpg.de/Genexpression_WS0506/material/wu2002.pdf).

Wu, J., Li, J., Salcedo, R., Mivechi, N.F., Trinchieri, G., Horuzsko, A., 2012. The proinflamma-tory myeloid cell receptor TREM-1 controls Kupffer cell activation and development of hepatocellular carcinoma. Cancer Res. 72,3977-3986.

Yamada, Y., Kirillova, I., Peschon, J.J., Fausto, N., 1997. Initiation of liver growth by tumor necrosis factor: deficient liver regeneration in mice lacking type I tumor necrosis factor receptor. Proc. Natl. Acad. Sci. U. S. A. 94,1441-1446.

Yang, L., Allen, B.C., Thomas, R.S., 2007. BMDExpress: a software tool for the benchmark dose analyses of genomic data. BMC Genomics 8, 387.

Young, M.R., Yang, H., Colburn, N.H., 2003. Promising molecular targets for cancer prevention: AP-1, NF-kB and Pdcd4. Trends Mol. Med. 9,36-41.

Yuan, L., Kaplowitz, N., 2009. Glutathione in liver diseases and hepatotoxicity. Mol. Aspects Med. 30, 29-41.

Yuan, R., Jeng, Y., Chen, H., Lai, P., Pan, H., Hsieh, F., Lin, C., Lee, P.H., Hsu, H., 2006. Stathmin overexpression cooperates with p53 mutation and osteopontin overexpression, and is associated with tumour progression, early recurrence, and poor prognosis in hepa-tocellular carcinoma. J. Pathol. 209, 549-558.

Zhang, Y., Chen, F., 2004. Reactive oxygen species (ROS), troublemakers between nuclear factor-kappaB (NF-kappaB) and c-Jun NH(2)-terminal kinase (JNK). Cancer Res. 64, 1902-1905.