The environmental contaminant 2,3,7,8-tetrachlorodibenzo- p -dioxin (TCDD) elicits dose-dependent hepatotoxicity that includes fat accumulation, inflammation, and fibrosis that may progress to hepatocellular carcinoma. To further investigate these effects, RNA-Seq data were integrated with computationally identified putative dioxin response elements, and complementary targeted metabolomic and aryl hydrocarbon receptor (AhR) ChIP-Seq data from female C57BL/6 mice gavaged with TCDD every 4 days for 28 days. Data integration using CytoKEGG with manual curation identified dose-dependent alterations in central carbon and amino acid metabolism. More specifically, TCDD increased pyruvate kinase isoform M2 (PKM2) gene and protein expression. PKM2 has lower catalytic activity resulting in decreased glycolytic flux and the accumulation of upstream intermediates that were redirected to the pentose phosphate pathway and serine/folate biosynthesis, 2 important NADPH producing pathways stemming from glycolysis. In addition, the GAC:KGA glutaminase (GLS1) protein isoform ratio was increased, consistent with increases in glutaminolysis which serves an anaplerotic role for the TCA cycle and compensates for the reduced glycolytic flux. Collectively, gene expression, protein, and metabolite changes were indicative of increases in NADPH production in support of cytochrome P450 activity and ROS defenses. This AhR-mediated metabolic reprogramming is similar to the Warburg effect and represents a novel advantageous defense mechanism to increase anti-oxidant capacity in normal differentiated hepatocytes.
2,3,7,8-Tetrachlorodibenzo- p -dioxin (TCDD) is the prototypical ligand for a series of structurally diverse environmental contaminants, natural products, and endogenous metabolites that bind and activate the aryl hydrocarbon receptor (AhR), a basic helix-loop-helix ligand activated transcription factor ( Denison et al. , 2011 ). Binding of TCDD and related compounds to the cytosolic AhR causes chaperone protein dissociation, followed by nuclear translocation and heterodimerization with the AhR nuclear translocator (ARNT). The TCDD-AhR-ARNT complex binds dioxin response elements (DREs, also known as AHREs and XREs) consisting of the invariant core sequence 5′-GCGTG-3′, eliciting changes in gene expression, although DRE-independent interactions have also been described ( Beischlag et al. , 2008 ; Dere et al. , 2011b ; Huang and Elferink, 2012 ).
TCDD and related compounds elicit a variety of adverse effects including teratogenesis, immunosuppression, and tumor promotion ( Pohjanvirta and Tuomisto, 1994 ), and have also been implicated in the development of non-alcoholic fatty liver disease and type II diabetes (T2D) ( Taylor et al. , 2013 ). In mice, a single bolus dose of TCDD results in reversible hepatic fat accumulation ( Boverhof et al. , 2005 ; Lee et al. , 2010 ) while repeated treatment sustains lipid accumulation with progression to steatohepatitis and fibrosis ( Nault et al. , 2015a ; Pierre et al. , 2014 ). TCDD-elicited metabolic disruption is not limited to lipid metabolism, but also disrupts glycolysis, the TCA cycle, and nucleoside metabolism ( Forgacs et al. , 2012 ; Lin et al. , 2011 ; Matsubara et al. , 2012 ). However, these were short term exposure studies involving a single bolus dose, or did not include complementary ‘omic’ data.
The objective of this study was to integrate complementary AhR ChIP-Seq, hepatic RNA-Seq, and hepatic and serum metabolomic data with computationally identified putative dioxin response elements (pDREs), to further elucidate the mechanisms involved in TCDD-elicited hepatotoxicity. Changes in central carbon and amino acid metabolism were reminiscent of the Warburg effect. This included pyruvate kinase isoform switching ( Christofk et al. , 2008 ), enhanced glutaminolysis ( DeBerardinis et al. , 2007 ), and increased serine (Ser)/folate metabolism ( Locasale, 2013 ; Maddocks et al. , 2013 ). This metabolic reprogramming is consistent with increased anti-oxidant capacity and may represent an advantageous, more general response, in normal differentiated tissue following toxic insult to support cell survival.
MATERIALS AND METHODS
Postnatal day 25 (PND25) female C57BL/6 mice (Charles River Laboratories, Portage, Michigan), housed in polycarbonate cages with cellulose fiber chips (Aspen Chip Laboratory Bedding, Warrensberg, New York) at 30–40% humidity and a 12-h light/dark cycle and acclimatized for 4 days, were fed ad libitum (Harlan Teklad 22/5 Rodent Diet 8940, Madison, Wisconsin) with free access to deionized water. On PND 28 and every following fourth day animals (n = 5) were orally gavaged with 0.1 ml sesame oil or 0.01, 0.03, 0.1, 0.3, 1, 3, 10, and 30 µg/kg TCDD (Dow Chemical Company, Midland, Michigan) for a total of 28 d. Blood was collected by submandibular vein puncture into BD microtainer tubes with serum separator and centrifuged at 10 000×g for 5 min and stored at −80°C. Liver samples were frozen in liquid nitrogen, and stored at −80°C. Only 30 µg/kg samples were analyzed by ChIP-Seq at 2 h. Livers from 92 days treated mice were collected using the same treatment and collection methods. Hepatic samples from female C57BL/6 mice treated with 100 µg/kg β-naphthoflavone (βNF) were collected following daily gavage for 28 days. TCDD time course, 3,3′,4,4′,5-pentachlorobiphenyl (PCB126), 2,3,7,8-tetrachlorodibenzofuran (TCDF), 2,2′,4,4′,5,5′-hexachlorobiphenyl (PCB153) were collected as described ( Boverhof et al. , 2005 ; Kopec et al. , 2010a , b ). All procedures were approved by the All-University Committee on Animal Use and Care.
pDREs were updated using the mm10 GRCm38 genome build ( Dere et al. , 2011a ). The position weight matrix used the same mouse and rat bona fide DRE’s ( Dere et al. , 2011a ) with Gsta2 removed (ie, not present in the most recent rat genome build (rn5)). Four new functional DREs for human Bach2 and Cyp1b1 were added ( Supplementary Table S1 ). pDREs UCSC genome browser tracks are available at http://dbzach.fst.msu.edu/index.php/supplementarydata.html .
Hepatic AHR ChIP-Seq
ChIP-Seq analysis was performed on liver samples collected at 2 h following a single oral gavage with 30 µg/kg TCDD. Cross-linked DNA was immunoprecipitated with either rabbit IgG or rabbit IgG and anti-AhR as previously described ( Dere et al. , 2011b ; Lo and Matthews, 2012 ). Libraries prepared using the MicroPlex kit (Diagenode) were pooled and sequenced at a depth of approximately 30 M on an Illumina HiSeq 2500 at the MSU Research and Technology Support Facility Core. Read processing and analysis was performed using the MSU High Performance Computing Center. Quality was determined using FASTQC v0.11.2 and adaptor sequences removed using Cutadapt v1.4.1 while low-complexity reads were cleaned using FASTX v0.0.14. Reads were mapped to the mouse reference genome (GRCm38 release 76) using Bowtie 2.0.0 and alignments were converted to SAM format using SAMTools v0.1.19. Normalization and peak calling was performed using CisGenome ( Ji et al. , 2008 ) determined by comparison of IgG control and AhR enriched samples (n = 5) using a bin size (−b) of 25 and boundary refinement resolution (−bw) of 1 with default parameters.
Liver (approximately 100 mg) and serum (30 µl) metabolites were extracted using chloroform:water:methanol from animals orally gavaged every 4 days for 28 days with sesame oil vehicle, 3, 10, or 30 µg/kg TCDD. Briefly, frozen liver samples were homogenized (Polytron PT2100, Kinematica) or vortexed (serum) in HPLC grade methanol:water (5:3, 4.63 mL, −20°C) in Pyrex glass tubes. HPLC-grade chloroform was added following homogenization (methanol:water:chloroform, 5:3:5), vortexed, shaken on ice for 10 min, and centrifuged at 3000×g. The transferred aqueous phase was dried under nitrogen gas at room temperature and resuspended in HPLC-grade water prior to analysis. Hepatic extract protein layers were dried and a bicinchoninic acid assay (Sigma-Aldrich) on a Tecan Infinite 200 microplate reader (Männedorf, Switzerland) was used to normalize signal in liver.
Amino acids were carbobenzyloxy derivatized for analysis. Redissolved solutions (20 µl) were added to 80 µl of HPLC-grade methanol, and 2.5 µl of triethylamine. 0.5 µl benzylchloroformate was added and samples were vortexed then centrifuged, and the supernatant was transferred to HPLC vials for analysis.
Targeted Metabolomic Analysis
Hepatic extracts were first loaded onto a trapping column (C18, 4 × 2 mm, Phenomenex) and washed for 30 s with HPLC grade water containing 10 mM tributylamine and 15 mM acetic acid for desalting. Hepatic extracts were examined by liquid chromatography (LC) and tandem mass spectrometry (MS/MS) using a Paradigm MS4 HPLC (Michrom Bioresources, Auburn, California), and a Synergi Hydro column (4 µm particle size, 80 Å, 150 × 2 mm, from Phenomenex) ( Lunt et al. , 2015 ). HPLC was coupled with negative-mode electrospray ionization to a TSQ Vantage Triple Stage Quadrupole Mass Spectrometer (Thermo Scientific) operating in multiple reaction monitoring (MRM) mode. Serum extracts were analyzed using an Acquity UPLC system (Waters), Ascentis Express column (C18, 5 cm × 2.1 mm, 2.7 µm, Sigma-Aldrich) while MRM was performed using a Xevo TQ-S triple quadrupole mass spectrometer (Waters) (see Supplementary Materials ). Raw data and MRM parameters are deposited in www.ebi.ac.uk/metabolights (MTBLS225).
Data analysis, including peak integration was performed using MAVEN ( Clasquin et al. , 2012 ; Melamud et al. , 2010 ). Waters raw data were converted to mzxml format using msconvert in the ProteoWizard Tools ( Kessner et al. , 2008 ) prior to analysis by MAVEN. One-way ANOVA was performed in SAS 9.3 (SAS Institute, Inc., Cary, North Carolina).
Hepatic Gene Expression
Published RNA-Seq data for mice orally gavaged with sesame oil vehicle, 0, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, or 30 µg/kg TCDD every 4 days for 28 days (GSE62902) were used for the integrative analysis ( Nault et al. , 2015b ). Filtering criteria of |fold-change| ≥ 1.5 and P1( t ) ≥ 0.8 were used for all data analysis. UCSC Genome browser tracks combine all reads from vehicle or 30 µg/kg doses ( http://dbzach.fst.msu.edu/index.php/supplementarydata.html ). Quantitative real-time PCR (QRTPCR) was performed as previously described ( Boverhof et al. , 2005 , see Supplementary Materials ). Relative expression was determined using the 2 −ΔΔCT method and the geometric mean of reference genes ActB , Gapdh , and Hprt . Reference genes showed no evidence of treatment related mRNA effects. Statistical analyses of QRTPCR data were performed using a 1-way ANOVA for dose-response study designs followed by Dunnett’s posthoc test. For time-course designs a 2-way ANOVA with dose and time as factors was used followed by Tukey’s posthoc test. For chemical treatment, each study was considered independently and analyzed by Student’s t test. Statistical analyses were performed using SAS 9.3.
Liver samples from animals orally gavaged with either sesame oil or 30 µg/kg TCDD (approximately 25 mg, n = 5) were ground by mortar and pestle in radioimmunoprecipitation buffer supplemented with protease inhibitor cocktail (Sigma-Aldrich), sonicated on ice, and centrifuged. Total protein in the supernatant was quantitated by bicinchoninic acid assay (Sigma-Aldrich). Protein levels were determined using the Wes (ProteinSimple, San Jose, Califrnia) system. PKM1 (Sigma-Aldrich; 1:400), PKM2 (Cell Signaling, Danvers, Massachusetts; 1:65), PKM1/2 (Cell Signaling; 1:65), and GLS1 (ProteinTech, Chicago, Illinois; 1:30) antibodies were duplexed with β-actin (Cell Signaling; 1:65) as loading control and detected by goat anti-rabbit secondary antibody. Chemiluminescence signals were exported from Compass software (ProteinSimple) and used to calculate fold-change. Statistical analyses were performed using Student’s t test in SAS 9.3.
Enrichment Analyses and KEGG Pathway Integration
Over-represented metabolites within a pathway, as well as impact scores representing the connectivity to other metabolites or clusters were determined using MetPA ( Xia and Wishart, 2010 ). Over-represented differentially expressed genes (DEGs) compared with all mouse Kyoto Encyclopedia of Genes and Genomes (KEGG; http://genome.jp/kegg/ ) pathway genes was performed using ClusterProfiler in R 3.2.0 ( Yu et al. , 2012 ). Pathways were manually curated to include only normal metabolism pathways which were considered relevant to the expression dataset (eg, signaling, tissue specific, and cancer pathways were removed).
The CytoKEGG 0.0.5 plugin (Cytoscape 3.2.0) was used to generate pathway templates. Through manual curation for completeness, and refinement, pathways were streamlined by combining reactions with the same compounds as a single node. pDRE, AhR enrichment, differential gene expression, and altered metabolite level data represent the most significant change for each node (eg, highest MSS score or largest |fold-change|). An expandable interactive versions can be viewed at http://dbzach.fst.msu.edu/index.php/supplementarydata.html .
Alterations in Central Carbon, Amino Acid, and Lipid Metabolism
TCDD elicits significant lipidomic changes consistent with the development of fatty liver and steatohepatitis in mice ( Boverhof et al. , 2005 ; Fernandez-Salguero et al. , 1996 ). In this study, complementary targeted metabolomics was used to examine aqueous metabolite changes in hepatic extracts and serum from mice treated with 3, 10, and 30 µg/kg TCDD every 4 days for 28 days. Of the approximately 190 metabolites examined, only 124 and 84 were detected in hepatic and serum samples, respectively ( Supplementary Table S2 and S3 ) with 48 hepatic and 41 serum metabolites altered ( P ≤ .05) by treatment. This included an 11-fold increase in hepatic oxaloacetate (OAA) and >1000-fold decrease in serum xanthine. Altered metabolite levels were analyzed for pathway over-representation (enrichment) and connectivity within related metabolites (impact) using MetPA ( Xia and Wishart, 2010 ). Central carbon (glycolysis/gluconeogenesis, the TCA cycle, and the pentose phosphate pathway [PPP]), and amino acid metabolism consistently exhibited higher enrichment (−log( p )) and/or impact scores ( Figs. 1 A and B, Supplementary Tables S4 and S5 ).
RNA-Seq ( Nault et al. , 2015b ) and metabolomic data were integrated using KEGG metabolic pathways. Analysis of 3406 DEGs (|fold-change| ≥ 1.5 and P1( t ) ≥ .8) using clusterProfiler ( Yu et al. , 2012 ) identified over-represented ( P ≤ .05) pathways associated with altered glucose, amino acid, and fatty acid metabolism ( Figure 1 C). KEGG enrichment also identified drug, starch and sucrose, and pyruvate metabolism enrichment. Steroid hormone biosynthesis was the most enriched pathway, consistent with previous reports ( Angrish et al. , 2013 ; Tanos et al. , 2012b ) although our targeted metabolomic data did not include steroid metabolites. Overall, metabolite changes were consistent with the identified DEGs.
To further evaluate AhR-mediated effects, ChIP-Seq was performed 2 h post dose with 30 µg/kg TCDD. Approximately 13 936 enrichment peaks were identified, of which 5766 satisfied the FDR ≤ 0.05 criteria, with 2711 located within the coding region or 10 kb upstream region ( Supplementary Table S6 ). Comparing the 3406 DEGs with 2711 AhR enriched regions with 23 054 pDREs identified 685 genes at the intersection. Only 475 had a pDRE within an AhR enriched region ( Figure 1 D) suggesting AhR-mediated DRE-dependent regulation, while 210 genes with AhR enrichment lacked an overlapping pDRE.
Increased Hepatic Lipid Uptake but Inhibited Fatty Acid Synthesis and Oxidation
KEGG enrichment analysis identified DEGs associated with lipid transport, processing, and metabolism ( Figure 1 E). This included the down-regulation of apolipoproteins (ie, 3.3-fold Apoa1 ) involved in lipid packaging and transport. Conversely, Apol7c , Apol10d , and Fabp12 were induced 3.6-, 14-, and 46-fold, respectively. Cd36 involved in lipid uptake was also induced 3.2-fold as previously reported ( Lee et al. , 2010 ). However, AhR enrichment with a pDRE was only observed for Fabp12 . Similarly, lipid processing genes such as elongases ( Elovl2 , 3 , and 6 ) and desaturases ( Fads , Scd ) were repressed >2-fold. Lipases Lpl , Mgll , Pnliprp1 , and Pnliprp2 were induced 4.8-, 2.1-, 48-, and 1.9-fold, respectively, while Lipa , Pnpla3 , and Pnpla5 were repressed 1.5-, 8.0-, and 2.4-fold, respectively. Overall, de novo fatty acid (FA) synthesis and β-oxidation genes were repressed ( Figure 1 E and Supplementary Figure S1 ) including a 5.3-fold repression of fatty acid synthase ( Fasn ) with AhR enrichment that included a pDRE. β-Oxidation was marked by repression of acetyl-CoA acyltransferase ( Acaa1b 3.1-fold), carboxylase ( Acaca 2.5-fold and Acacb 1.5-fold), and enoyl-CoA hydratase ( Ehhadh 1.6-fold). Acsl1 , Acsm3 , Acaa1b , Acaca , and Ehhadh exhibited AhR enrichment with a pDRE, except for Acsm3 . The inhibition of β-oxidation is consistent with short acyl-CoA reductions, particularly decreases in butanoyl-CoA (2.9-fold).
Induction of Oxidative Stress and Anti-Oxidant Defenses
CytoKEGG analysis with manual curation identified the enrichment of oxidative stress (eg, P450 induction, glutathione metabolism), central carbon metabolism, and amino acid metabolism pathways ( Figs. 2 and 3 ). These pathways were found to converge on NADP + /NADPH, indicating altered production and/or consumption of reducing equivalents.
The 1250-, 807-, and 16-fold induction of classical AhR battery genes Cyp1a1 , Cyp1b1 , and Cyp1a2 as well as the 1.6- and 1.9-fold increase in NADPH oxidase ( Nox4 ) and nitric oxide synthase ( Nos2 ) represent major NADPH consumers as well as producers of reactive oxygen species (ROS; Figure 2 A). In response, glutathione peroxidase ( Gpx2 ) was induced 11-fold, an important defense mechanism that uses glutathione (GSH) to reduce lipid peroxide and hydrogen peroxide oxidative stress products to alcohols and water, respectively ( Figure 2 B). The resulting oxidized glutathione (GSSG) is then recycled to GSH (increased 1.4-fold) using NADPH by glutathione reductase ( Gsr ), which was induced 1.4-fold, to preserve cellular anti-oxidant capacity. TCDD also induced a 1.7-fold increase in glutathione synthetase ( Gss ), which catalyzes the condensation of glycine (Gly) with γ-glutamylcysteine to produce more GSH ( Figure 2 B). Accordingly, these processes increased total hepatic glutathione levels 1.3-fold ( Supplementary Figure S2 ). The absence of AhR enrichment within Gpx2 , Gsr, and Gss loci suggest these responses were mediated by oxidative stress as opposed to direct AhR regulation, and likely involve Nrf2 signaling which is known to be induced following AhR activation ( Lu et al. , 2011 ; Wang et al. , 2011b ; Yeager et al. , 2009 ).
In addition, purine metabolism, producing the anti-oxidant, uric acid, was increased 1.5-fold in the liver ( Figure 2 ). Inosine monophosphate (IMP) is converted to inosine by 5′ nucleotidase ( Nt5e ), and then to hypoxanthine by purine nucleoside phosphorylase ( Pnp ) and subsequently to xanthine and uric acid by xanthine dehydrogenase ( Xdh ). Although decreases in serum inosine (33-fold), hypoxanthine (170-fold), and xanthine (>1000-fold), are consistent with the 2.1- and 1.9-fold AhR-mediated induction of Nt5e and Xdh ( Figure 2 C), increased enzymatic activity are likely the primary factors. Interestingly, the oxidase form ( XO ) of Xdh irreversibly produced under low energy conditions is also a significant ROS contributor.
Interaction of Glucose, Amino Acid, and Glutathione Metabolism
Alterations in the catabolic oxidation of glucose were also revealed by enrichment analysis ( Figs. 3 A and B). Although hexose transporters with different kinetic properties and substrate preferences were differentially expressed ( Slc2a3 , 2a6 , 2a10 , and 2a3 were induced 1.6-, 2.0-, 1.7-, and 1.6-fold, respectively while Slc2a2, 2a4, and 2a5 were repressed 2.1-, 3.0-, and 1.8-fold, respectively), only Slc2a2 was associated with AhR enrichment. Hexokinases ( Hk ) 1, 2, and 3 were induced 2.3-, 2.5-, and 2.2-fold, respectively, while glucokinase ( Gck ) was repressed 4.2-fold. HK2 induction is consistent with the Warburg effect. Glucose-6-phosphate (G6P) levels in liver extracts were reduced 1.9-fold, consistent with a approximately 30% decrease at 168 h after a single bolus dose of TCDD ( Forgacs et al. , 2012 ). There was a negligible decrease in glucose-6-phosphate isomerase ( Gpi ) mRNA with a 1.9-fold decrease in fructose-6-phosphate levels (F6P). Liver phosphofructokinase ( Pfkl ), which catalyzes the production of fructose 1,6-bisphophate (F1,6BP) was also repressed 1.5-fold while fructose 1,6 bisphosphatase ( Fbp1 ) catalyzing the reverse reaction was induced 1.7-fold. Moreover, phosphofructokinase-2/fructose-2,6-bisphosphatase ( Pfkfb3 ) was induced 1.8 fold, which could decrease fructose-2,6-bisphosphate (F2,6BP) levels, a potent phosphofructokinase allosteric activator ( Chesney, 2006 ). Fructose bisphosphate aldolase B ( Aldob ) expression was also repressed 2.1-fold. Although cytosolic glyceraldehyde-3-phosphate dehydrogenase ( Gpd1 ) which converts DHAP to glycerol-3-phosphate to support triglyceride synthesis was unaffected, mitochondrial Gpd2 was repressed 2.8-fold which may compromise glycerol-3-phosphate re-oxidation to DHAP and FADH 2 production, limiting cytosolic NADH electron transfer to support oxidative phosphorylation. TCDD did not affect glyceraldehyde-3-phosphate dehydrogenase ( Gapdh ) or phosphoglycerate mutase ( Pgam1 and 5 ) expression while modestly inducing phosphoglycerate kinase ( Pgk1 , 1.3-fold), and repressing enolase ( Eno1 , 1.6-fold) which catalyzes the transformation of 2-phosphoglycerate to phosphoenolpyruvate (PEP). Both Pgam1 and Eno1 showed AhR enrichment.
PEP conversion to pyruvate by pyruvate kinase is a rate limiting step in glycolysis. Liver/red blood cell pyruvate kinase ( Pklr ) was repressed 6.7-fold, consistent with a reported 3.2-fold decrease ( Forgacs et al. , 2012 ) while muscle pyruvate kinase ( Pkm ) was induced 4.6-fold with AhR enrichment, as previously reported ( Forgacs et al. , 2012 ). Further analysis of Pkm mapped RNA-Seq reads to exon 10 as opposed to exon 9 indicating isoform switching from Pkm1 to Pkm2 ( Figure 4 A). Pkm2 is typically expressed in rapidly proliferating cancer cells and is characteristic of the Warburg effect, with a lower catalytic rate in the dimeric or monomeric form, leading to the accumulation of upstream glycolytic intermediates. QRTPCR confirmed a 9-fold induction of Pkm1 compared with 40-fold for Pkm2 ( Figs. 4 B and C). Protein measurements found a 13.4-fold increase in total PKM, with a 1.6-fold increase in PKM1, and a 3-fold increase in PKM2 using 3 different antibodies ( Figs. 4 D and F). Although Pkm2 was only induced 2.4-fold compared with 1.7-fold for Pkm1 at 72 h, it was induced 13.8-fold compared with 1.6-fold for Pkm1 after treatment every 4 days for 92 days ( Figs. 5 A and B). Other AhR agonists, including 3,3′,4,4′,5-pentachlorobiphenyl (PCB126), 2,3,7,8-TCDF, and βNF also induced Pkm2 , while 2,2′,4,4′,5,5′-hexachlorobiphenyl (PCB153), which does not bind to the AhR, had no effect on Pkm2 mRNA levels ( Figs. 5 C and D). ChIP-Seq analysis at 2 h after treatment with 30 µg/kg TCDD identified a 7.5-fold increase in AhR enrichment within the Pkm coding region that contained a pDRE suggesting DRE-dependent AhR-mediated regulation ( Figure 4 A). This same peak was previously identified by ChIP-chip and confirmed by ChIP-PCR ( Dere et al. , 2011b ).
Reduction of glycolytic flux and the inhibition of β-oxidation not only reduce hepatic acetyl-CoA levels, but also impact the TCA cycle and reducing equivalent production. Less pyruvate from glycolysis, inhibition of pyruvate dehydrogenase by pyruvate dehydrogenase kinase ( Pdk3 and 4 induced 1.9 and 1.7-fold, respectively), and the inhibition of β-oxidation reduce acetyl-CoA availability for condensation with OAA. In cancer cells and the Warburg effect, central carbon metabolism reprogramming also involves compensatory changes in glutamine metabolism (Gln) to maintain oxidative phosphorylation and ATP production ( van den Heuvel et al. , 2012 ). TCDD increased glutamate (Glu) levels in the serum, altered its metabolism ( Figure 1 ) and induced Slc1a5 (5.1-fold), Slc6a19 (3.6-fold), and Slc38a1 (1.5-fold) transporters while repressing the expression of antiporters Slc38a3 , 4 and 7 (repressed 1.8 -, 1.7 -, and 1.3-fold respectively; Figure 3 C). A concomitant 1.4-fold increase in glutaminase ( Gls1 ), which catalyzes the transformation of Gln to Glu and ammonia, and a 1.9-fold decrease in glutaminase 2 ( Gls2 ) are consistent with increased glutaminolysis, as reported in some cancer cell lines ( Wise et al. , 2008 ). Glutaminolysis may also provide Glu to support GSH biosynthesis.
Coincidentally, mRNA for the cystine (Cyss) transporters, xCT ( Slc7a11 induced 1.8-fold), Slc3a2 (induced 1.9-fold), and Slc7a6 (induced 1.8-fold) are increased. Slc7a11 is regulated by Nfe2l2 (aka Nrf2 ) ( Conrad and Sato, 2012 ) which was also induced 2.4-fold by TCDD. Once inside the cell, Cyss is reduced to Cys, the rate limiting precursor for de novo GSH biosynthesis ( Yin et al. , 2015 ). Although Cys may be synthesized from Met by the transsulfuration pathway, the primary source of intracellular Cys is more likely transport as hepatic Met levels were decreased 5.9-fold in the liver.
Two splice variants associated with tumor environments exist for Gls1 ; kidney-type glutaminase (KGA) consisting of exons 1–14 and 16–19, and the GAC isoform with an alternate carboxy-terminus consisting of exon 15. The GAC:KGA mRNA and protein ratio is increased in tumors ( van den Heuvel et al. , 2012 ), and also increased by TCDD, with GAC and KGA mRNA induced 2.3 - and 1.6-fold, respectively ( Figs. 3 D and E). Accordingly, the GAC:KGA protein ratio was increased 1.7-fold ( Figure 3 F).
Higher hepatic Glu levels (2.5-fold) would be expected to be converted to α-ketoglutarate (αKG) by glutamate dehydrogenase ( Glud1 ) to replenish the TCA cycle ( Figure 3 C), and may explain the 11-fold increase in OAA in the treated liver ( Figure 3 C). The 2-fold repression of PEP carboxykinase ( Pck1 ), the rate-limiting step of gluconeogenesis, is also consistent with gluconeogenesis inhibition ( Ahmed et al. , 2015 ; Diani-Moore et al. , 2010 ; Weber et al. , 1995 ), further contributing to OAA accumulation. However, Glud1 is repressed 1.8-fold as are succinyl-CoA synthetase (1.6 - and 1.5-fold repression of Sucla2 and Suclg2 1.5-fold), succinate dehydrogenase (1.3-fold repression of Sdhd ), and fumarase (1.4-fold repression of Fh1 ) which may reduce glutaminolysis flux towards OAA.
Alternatively, Glu could be used in the malate-aspartate shuttle to produce pyruvate, which is reported not to contribute to lactate production ( Son et al. , 2013 ). In this scenario, mitochondrial OAA and Glu are converted to αKG and aspartate (Asp) by glutamate oxaloacetate transaminase ( Got2 ) which is modestly repressed 1.3-fold ( Figure 3 F). The 1.9-fold repression of cytosolic Got1 would limit the reconversion of cytosolic Asp, and may be the source of the 4.1-fold increase in hepatic Asp levels.
Redirecting Central Carbon and Amino Acid Metabolism to NADPH Production
Increased ROS production and lipid peroxidation create an oxidative environment requiring NADPH-dependent counter-measures for cell survival. In cancer cells, NADPH demands are satisfied by the PPP with smaller contributions from NADP + -dependent malic enzyme ( Me1 ) and NADP + -dependent isocitrate dehydrogenase ( Idh ) ( Jiang et al. , 2014 ). With increased PKM2 expression, accumulating G6P levels are redirected to the PPP ( Figure 6 A) consistent with decreases in hepatic G6P (1.9-fold), F6P (1.9-fold), 6-phospho-D-gluconate (6PG; 4.5-fold), erythrose-4-phosphate (E4P; 2.1-fold), and 5-phosphoribose 1-diphosphate (1.6-fold). Increases in hepatic ribulose 5-phosphate (Ru5P; 1.5-fold), and ribose 5-phosphate (R5P; 2-fold) are consistent with increased PPP flux. In mammalian cells, oxidative stress is reported to induce glucose-6-phosphate dehydrogenase ( G6pdh ) activity so that this enzyme is not rate limiting ( Cosentino et al. , 2011 ; Wang et al. , 2014 ). Phosphoglucomutase ( Pgm1 ), involved in the catalysis of R5P to ribose 1-phosphate (R1P) was also induced 1.6-fold further supported suggestions of increased flux through the PPP that would produce NADPH.
Recently, folate-dependent reactions have been estimated to account for 10–40% of NADPH production, compared with approximately 30% from the PPP ( Fan et al. , 2014 ). PKM2 reduced glycolytic flux is consistent with the shunting of 3-phosphoglycerate (3PG) towards Ser biosynthesis and folate-dependent NADPH production ( Figure 6 B). 3PG is transformed to 3-phosphohydroxypyruvate by phosphoglycerate dehydrogenase ( Phgdh , 5.8-fold) and subsequently to Ser by phosphoserine aminotransferase ( Psat1 , 3-fold), which simultaneously converts Glu to αKG, another potential source of OAA accumulation ( Figure 3 C). Ser and tetrahydrofolate (THF) can then be converted to Gly and 5,10-methylene-tetrahydrofolate (5,10-CH 2 -THF) by cytosolic and mitochondrial serine hydroxymethyl transferase 1 and 2 ( Shmt1 / 2 ). TCDD did not affect Shmt2 expression but repressed Shmt1 2.3-fold. 5,10-CH 2 -THF is then used by methylenetetrahydrofolate dehydrogenase ( Mthfd1/2 ) to produce NADPH. Mitochondrial Mthfd2 was induced 2.4-fold while cytosolic Mthfd1 was modestly repressed suggesting greater NADPH production to support mitochondrial oxidative defenses. Besides de novo synthesis, Ser is available from protein and Ser-containing phospholipid catabolism, as well as transporter-mediated uptake. The 4.2-fold increase in serum Ser levels may be made available via the 3.4 - and 5.1-fold induction of Slc1a4 and 1a5 , respectively, 2 high affinity transporters. In addition to producing NADPH, this pathway also provides Gly, another precursor required for GSH biosynthesis, thus serving dual antioxidant roles.
Cell Cycle Arrest by TCDD
PKM2 lies at the intersection of cell proliferation and cell survival ( Cairns et al. , 2011 ; Gorrini et al. , 2013 ; Harris et al. , 2012 ). The AhR has also been linked to anti-oxidant defenses and cell cycle regulation ( Jackson et al. , 2014 ; Puga et al. , 2002 ). Jackson et al. (2014) report that p21 gene expression regulation ( Cdkn1a ) is critical to TCDD-induced G 1 cell cycle arrest. Similarly, we confirm Cdkn1a was induced 7.6-fold at 28 d with 3 AhR enriched regions, but was repressed 1.2-fold following treatment with TCDD every 4 days for 92 days ( Figure 7 ), while Pkm2 was induced 13.8-fold. Other cell cycle regulators also demonstrated altered expression by TCDD, although their role in TCDD-induced in cell cycle arrest is not as well understood including the 8.6-fold induction of Fhit , the 1.7-fold repression of Gadd45α and the 3.1-fold induction of Gadd45β .
Results from previous TCDD ‘omic’ studies in various in vivo and in vitro models have used different designs, doses, treatment protocols, and models that limit extrapolations between species, and confound integration and interpretation. Consequently, we examined the metabolomic and transcriptomic effects of TCDD in the same animals, and used KEGG pathways to integrate and manually curate the data while incorporating complementary pDRE and AhR-binding (ChIP-Seq) data to further investigate metabolic reprogramming in the mouse liver. Oral treatment every 4 days for 28 days simulated continuous exposure that increased the severity of the hepatic effects but was well tolerated ( Nault et al. , 2015a ). Our results demonstrate that TCDD elicited a dose-dependent reorganization of hepatic metabolism characterized by the reprogramming of central carbon and amino acid metabolism, reminiscent of the Warburg effect, which supported NADPH demand created by increased anti-oxidant defense activities.
TCDD induced hepatic steatosis by altering gene expression consistent with lipid accumulation. This also included the inhibition of de novo fatty acid synthesis and β-oxidation ( Boverhof et al. , 2005 ; Fernandez-Salguero et al. , 1996 ; Lee et al. , 2010 ; Sato et al. , 2008 ; Tanos et al. , 2012a ). Although other studies have suggested that AhR-mediated inhibition of β-oxidation may involve suppression of PPARα activity ( Shaban et al. , 2004 ; Wang et al. , 2011a ), the presence of AhR enrichment within the loci of several de novo fatty acid synthesis and β-oxidation genes suggest AhR regulation involving tethering mechanism ( Beischlag et al. , 2008 ; Dere et al. , 2011b ). Hepatic fat accumulation and the induction of xenobiotic (eg, cytochrome P450 activity) and purine metabolism (ie, Xdh /XO) would not only increase ROS levels but also lipid peroxidation ( Bansal et al. , 2014 ; Sugihara et al. , 2001 ), creating a sustained oxidative burden requiring counter-measures to minimize cell damage. In response, AhR and Nrf2 cooperate to induce antioxidant counter-measures ( Lu et al. , 2011 ; Wang et al. , 2011b ; Yeager et al. , 2009 ). Furthermore, redirection of glycolytic flux and the inhibition of β-oxidation would also compromise NADH and ATP production further compounding cell damage ( Lu et al. , 2011 ).
Notably, TCDD repressed Pklr expression while inducing Pkm2, the isoform typically associated with the Warburg effect and cell proliferation ( Lunt et al. , 2015 ). This is the first report of PKM isoform switching elicited by an exogenous chemical in normal differentiated tissue, although switching has been shown in PTEN -null mice fatty liver tissue ( Panasyuk et al. , 2012 ). PKM2 sits at the intersection of cell growth and cell survival ( Cairns et al. , 2011 ; Gorrini et al. , 2013 ; Harris et al. , 2012 ), and its induction with subsequent metabolic changes highlights PKM2’s putative role in TCDD toxicity. Although we did not determine whether Pkm2 was present in monomeric, dimeric, and tetrameric forms, under oxidative stress conditions such as those elicited by TCDD, PKM2 is oxidized at Cys358 which further reduces its activity resulting in the accumulation of upstream intermediates and shunting towards NADPH producing pathways ( Anastasiou et al. , 2011 ). Interestingly, AhR enrichment has been repeatedly demonstrated within the Pkm loci indicating direct regulation. Furthermore, PCB126, TCDF, and βNF also induced Pkm2 while PCB153, which does not bind to the AhR, had no effect on expression. Collectively, these results indicate AhR-mediated regulation of PKM2 expression is a sustained response as long as TCDD exposure is maintained.
In cancer cells, PKM2 is tightly regulated, both transcriptionally and post-translationally ( Cairns et al. , 2011 ; Gorrini et al. , 2013 ; Harris et al. , 2012 ). In response to TCDD, our data are consistent with the shunting of accumulating upstream intermediates towards anti-oxidant pathways as opposed to cell proliferation, especially given the induction of the cell cycle inhibitor p21 ( Jackson et al. , 2014 ). More specifically, G6P and 3PG are redirected to the PPP and Ser/folate biosynthesis, respectively, to satisfy increasing NADPH demand in order to maintain redox balance that include GSH biosynthesis and GSSG reduction. Similar metabolic reprogramming is reported in cancer cells ( Amelio et al. , 2014 ; Cairns et al. , 2011 ; Gorrini et al. , 2013 ; Wise et al. , 2008 ), lipopolysaccharide-activated macrophages ( Palsson-McDermott et al. , 2015 ), and normal cells experiencing oxidative stress ( Kuehne et al. , 2015 ).
Metabolic reprogramming, including increased intermediate transport, also explains increased hepatic Gly, Cys, OAA and Asp levels. 3PG accumulation due to PKM2 expression is shunted to Ser biosynthesis indicated by Phgdh and Psat1 induction. Ser is then converted to Gly to produce 5,10-CH 2 -THF which is oxidized to produce NADPH. Ser/folate metabolism is estimated to account for 10–40% of NADPH production, compared with approximately 30% from the PPP ( Fan et al. , 2014 ). Moreover, Ser biosynthesis involves the conversion of Glu to αKG, another potential source of OAA, in addition to glutaminolysis and the malate-aspartate shuttle. Decreases in acetyl-CoA, and the inhibition of β-oxidation and gluconeogenesis further contribute to OAA accumulation. Pck1 repression also inhibits OAA utilization in gluconeogenesis and metabolism to lactate. Increased OAA levels likely drive Asp biosynthesis despite Got repression. Nevertheless, more precise tracer studies are required to further characterize TCDD-elicited metabolic reprogramming since metabolic regulation primarily occurs via post-translational mechanisms ( Stincone et al. , 2014 ).
In summary, compromised energy production, increasing lipid peroxidation, and the resulting oxidative environment cause cell damage that triggers not only immune cell infiltration but also the reprogramming of central carbon and amino acid metabolism towards defensive responses ( Figure 8 ). Induction of p21 and its known role in TCDD mediated cell cycle arrest ( Jackson et al. , 2014 ) further supports cell survival by ensuring accumulating intermediates are shunted towards oxidative stress counter-measures and cell survival, as opposed to biomass production and cell proliferation. Resistance to cell death, energy metabolism reprogramming, and the initiation of inflammation are important hallmarks of cancer ( Hanahan and Weinberg, 2011 ). TCDD is a known carcinogen in rodents, but its human carcinogenicity is debated ( Boffetta et al. , 2011 ). The continued induction of PKM2, and return to basal Cdkn1a expression at 92 days suggests possible support of proliferation of surviving cells with accumulated DNA damage that have evaded repair, leading to tumor development. This balance between survival and proliferation regulated by metabolic reprogramming could also have roles in AhR-mediated species- and sex-specific sensitivity and ligand toxicity differences. Consequently, regulating metabolic reprogramming in normal tissue, reminiscent of the Warburg effect, may be an advantageous defensive strategy used by normal cells to promote survival in response to other toxicities that is later co-opted by transformed cells to support proliferation and tumor development.
We thank Dr Ed Dere for his help with computational pDRE identification, and Dr. John Froehlich and Brendan Johnson for ProteinSimple Wes assay assistance.
This work was supported by the National Institute of Environmental Health Sciences Superfund Research Program (NIEHS SRP P42ES04911). T.R.Z. is supported by AgBioResearch at MSU. RN is supported by the MSU Barnett Rosenberg Endowed Assistantship and Integrative Training in the Pharmacological Sciences grant (the National Institutes of Health 5T32GM092715). K.A.F. is supported by the Canadian Institutes of Health Research (CIHR) Doctoral Foreign Study Award (DFS-140386). S.Y.L. is supported by the Department of Defense (W81XWH-12-1-0466). JM is supported by a CIHR operating grant (MOP-125919), an unrestricted research grant from DOW Chemical, and an Early Researcher Award from the Ontario Ministry of Innovation.