Genome-Scale Characterization of Toxicity-Induced Metabolic Alterations in Primary Hepatocytes.

Context-specific GEnome-scale metabolic Network REconstructions (GENREs) provide a means to understand cellular metabolism at a deeper level of physiological detail. Here, we use transcriptomics data from chemically-exposed rat hepatocytes to constrain a GENRE of rat hepatocyte metabolism and predict biomarkers of liver toxicity using the Transcriptionally Inferred Metabolic Biomarker Response algorithm. We profiled alterations in cellular hepatocyte metabolism following in vitro exposure to four toxicants (acetaminophen, carbon tetrachloride, 2,3,7,8-tetrachlorodibenzodioxin, and trichloroethylene) for six hour. TIMBR predictions were compared with paired fresh and spent media metabolomics data from the same exposure conditions. Agreement between computational model predictions and experimental data led to the identification of specific metabolites and thus metabolic pathways associated with toxicant exposure. Here, we identified changes in the TCA metabolites citrate and alpha-ketoglutarate along with changes in carbohydrate metabolism and interruptions in ATP production and the TCA Cycle. Where predictions and experimental data disagreed, we identified testable hypotheses to reconcile differences between the model predictions and experimental data. The presented pipeline for using paired transcriptomics and metabolomics data provides a framework for interrogating multiple omics datasets to generate mechanistic insight of metabolic changes associated with toxicological responses.

Toxicity is an unintended effect of many compounds, resulting in significant health complications. The liver, kidney, and heart are often subject to adverse, potentially toxic effects because of their role in drug metabolism (Albini ; Awdishu and Mehta, 2017; Chen ). Hepatotoxicity is of particular concern (Church and Watkins, 2017; Rueda-Zárate ; Zimmerman, 1999), highlighting the need to understand how liver metabolism is altered as a result of toxicity. Understanding the metabolic changes to the liver can facilitate understanding the mechanisms associated with toxicity, thereby guiding development of novel strategies to counterbalance any toxic effects. Furthermore, with such mechanistic interrogation of liver metabolism, we can identify potential biomarkers associated with toxicity and potential intervention points involved with toxicological processes.

Genome-scale metabolic network reconstructions (GENREs) have emerged as useful tools for the study of cellular metabolism (Gille ; Karlstädt ; Mardinoglu ; Väremo ). GENREs represent metabolic reactions in a stoichiometric matrix that accounts for the stoichiometric coefficients of chemical transformations and the associated metabolites. GENREs also account for gene-protein-reaction (GPR) rules that map relationships between genes, the proteins they encode, and the reactions they catalyze in the network. With the GPR mappings and stoichiometric matrix to account for associated metabolic reactions, GENREs can be used to predict gene essentiality, changes in metabolites secreted, and the ability of a cell to catabolize particular carbon substrates; because of these characteristics, GENREs are increasingly applied to tackle questions about cellular toxicological responses (Bartell ; Brunk ; Carbonell ; Gatto ; Pannala ).

The incorporation of omics data into GENREs allows for cell-type-specific interrogation of metabolism. Transcriptomics and proteomics data are frequently integrated into GENREs to create cell-type-specific models. Several algorithms to integrate omics data into GENREs have been developed (Shlomi ; Zur ). Often with such methods, the integration of omics data constrains the GENRE by turning “on” and “off” genes and their associated reactions, reflecting gene expression in different conditions. These expression data integration algorithms help to contextualize these omics data and improve predictions of cellular metabolic functions.

Biomarkers are currently used in the diagnosis of cancer, cardiac function, and renal function (Jungbauer ; Lotan ; Pan ; Shlipak ) among other pathologies, often associating the presence or absence of a molecule with a specific diagnosis. For example, alanine aminotransferase is a protein that is used frequently as a biomarker of liver function (Dufour ; Zimmerman, 1999); high levels of this protein indicate that the liver has been damaged. A recently developed computational method for predicting biomarkers called Transcriptionally Inferred Biomarker Response (TIMBR) (Blais ) uses gene expression data contextualized in a GENRE to estimate relative changes in secreted metabolite levels. In a previous study (Blais ), TIMBR predicted changes in extracellular metabolite levels based on gene expression data for cells exposed to various chemical compounds. Predictions of a limited number of metabolite biomarkers for one chemical were validated, but a global evaluation of how well the biomarker predictions matched experimental data was missing. In this study, predictions from TIMBR are compared with paired metabolomics data to observe the differences between computational predictions and experimental data. Agreement between predictions and experimental data can be illustrative of mechanism behind an observed biomarker; disagreements between the computational model and experimental data can facilitate the development of specific testable hypotheses.

Here, we exposed primary rat hepatocytes to four chemical compounds and characterized their acute metabolic response (Figure 1). After exposure, transcriptomics and metabolomics data were collected from the same sample. We characterized the response of the hepatocytes to the compounds through changes in gene expression and metabolite levels, and evaluated similarities and differences between the cell’s responses across all conditions. The transcriptomics data were integrated into a GENRE of rat metabolism via iMAT (Zur ) to create a hepatocyte-specific network model, then the TIMBR algorithm was used to predict changes in the secreted metabolite profile. We compared these predictions with the coupled metabolomics data. With this methodology, we present a comprehensive strategy to characterize the toxicological response of hepatocytes to compounds of interest, and provide a framework to identify further areas of study in hepatocyte drug and toxicity metabolism.

Figure 1.

Schematic of the experimental set up. A, Primary rat hepatocytes were plated in 12-well format and exposed to acetaminophen, carbon tetrachloride, 2,3,7,8-tetrachlorodibenzodioxin (TCDD), or trichloroethylene for six h. After compound exposure, supernatants were collected and sent for metabolomics analysis. Hepatocytes were lysed and RNA was collected for sequencing. B, Cellular RNA was isolated and sequenced by Genewiz. With the raw sequencing reads as an input, the program kallisto was used to align sequencing reads to a reference transcriptome. The R packages TxImport and DESeq2 were used to summarize transcript counts to the gene level and to perform differential gene analysis respectively. Spent media from the hepatocytes were collected and sent for GC-MS, LC-MS, and HILIC-QTOF metabolomics at West Coast Metabolomics. After receiving metabolite peak intensities, the data were processed in R to generate a list of differentially abundant metabolites in each condition.

MATERIALS AND METHODS

Hepatocyte growth conditions

Frozen, primary rat hepatocytes (male, Sprague-Dawley) were purchased from ThermoFisher Scientific and cultured according to the manufacturer’s directions. Briefly, cells were rapidly thawed in a water bath (37°C), resuspended in plating media (William’s E media base supplemented with FBS, dexamethasone, penicillin/streptomycin, insulin, GlutaMAX, and HEPES; Gibco #CM3000), pelleted (50 ×g, five min), and plated at approximately 85% confluence in 12-well tissue culture plates. After 24 h, plating media was replaced with maintenance media (William’s E media base supplemented with dexamethasone, penicillin/streptomycin, ITS+, GlutaMAX and HEPES; Gibco #CM4000) and cells were incubated at 37°C under 5% CO2 for the remainder of the experiment.

Hepatocyte exposure to compounds

Hepatocytes were exposed to a hepatotoxicant and general toxicants at subtoxic levels. Subtoxic levels were defined as concentrations that resulted in minimal cell death but observed phenotypic changes (e.g., decreases in albumin production, ATP levels, increases in cytochrome p450 activity) (Supplementary Figs 1-3). The compounds were acetaminophen (APAP) at 3 mM, carbon tetrachloride (CCl4) at 10 mM, trichloroethylene (TCE) at 1 mM, and 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) at 1 nM. APAP and CCl4 are known hepatotoxicants, whereas TCDD and TCE are not considered primary hepatotoxicants typically, but do indeed induce hepatotoxicity. APAP, TCDD, and TCE conditions have four replicates, whereas CCl4 and the DMSO controls have three replicates. Solutions were made in WEM containing 0.1% DMSO with 0.1% DMSO as a control. Cells were exposed to the compounds for six hour. Concentrations and the six hour time point were selected based on literature evidence of comparable studies and conditions (Aly and Domènech, 2009; Cai ; Dere ; Forgacs ; Kienhuis ; Mitchell ; Uehara ; Xu ).

RNA isolation, sequencing, and analysis

After supernatants were collected, cells from each condition were treated with TRIzol and then scraped and collected into tubes. Chloroform was added to each tube and after shaking, cells were poured into prespun phase-lock gel tubes (5PRIME). Tubes were then spun in a cold room, the upper phase was collected, and isopropanol and glycogen were added to each tube followed by gentle inversion. Supernatants were again spun in a cold room and the resulting pellet was washed twice with 75% ethanol. The pellet was semidried and then dissolved in nuclease-free H2O. RNA samples were treated with DNA-free DNA removal kit (Ambion/Invitrogen), according to the manufacturer’s instructions, to remove any remaining DNA. RNA was quantified using the Qubit RNA broad range kit and sample integrity assessed using Agilent. RNA samples were subjected to rRNA depletion prior to library construction and sequencing; all services were performed by GENEWIZ. Libraries were sequenced using the Illumina HiSeq2500 platform in a 2× 100 bp pair-end configuration in High Output mode (V4 chemistry). The Unix-based program Kallisto v. 0.43.0 (Bray ) was used to process RNA sequence data in fastq format and quantify transcript abundances. Normalized transcript abundance values (transcripts per million) were calculated by Kallisto, using default settings, and imported to R for differential analysis. To quantify transcript abundances and aggregate toward the gene level, the package tximport in R was used (Soneson ). Differential gene expression was then performed with the standard DESeq2 R package (Love ) to obtain a list of differentially expressed genes (DEGs) with their log2-fold change values.

Metabolomics

After hepatocytes were exposed to the different compounds, supernatants were collected and stored at −20°C. Supernatants were then shipped to West Coast Metabolomics (http://metabolomics.ucdavis.edu/) at the University of California, Davis and untargeted analysis of primary metabolites, complex lipids, and biogenic amines was conducted on each sample, DMSO controls, and on blank media. An extraction solvent of 3:3:2 acetonitrile/isopropanol/water was prepared to use with the collected samples for Gas Chromatography-Mass Spectrometry (GC-MS) to analyze primary metabolites. External and internal standards for quality control were also prepared along with the samples. Raw results were reported as peak heights for quantification ion at the specific retention index. A full description of the protocol was outlined previously (Fiehn, 2016). Lipidomics analysis was performed by preparing samples with methanol, methyl tert-butyl ether, and water before running Liquid Chromatography-Mass Spectrometry (LC-MS). LipidBlast was used to identify and annotate lipids, and peak heights were reported according to the published protocol (Cajka and Fiehn, 2017).

Biogenic amine peak heights were quantified using Hydrophilic Interaction Chromatography Quadrupole Time of Flight (HILIC-QTOF) Mass Spectrometry, and peak heights were calculated following methods previously described (Meissen ). Samples were processed and analyzed according to West Coast Metabolomics protocols. Proteins and small polar hydrophilic small molecules were separated from lipids according to the protocol published by Matyash . Data were acquired using the following chromatographic parameters. Ultrapure water with 10 mM ammonium formate and 0.125% formic acid (pH 3) for mobile phase A, and 95:5 (vol/vol) acetonitrile: ultrapure water with 10 mM ammonium formate with 0.125% formic acid (pH 3) for mobile phase B. A column temperature of 40°C, with the flow rate of 0.4 ml/min and injection volume of 3 µl for ESI (+) and temperature of 4°C was used. The ESI Capillary voltage was +4.5 kV for ESI (+), the scan range was m/z 60−1200 Da, and the mass resolution was 10 000 for ESI (+) on an Agilent 6530 QTOF MS. After raw peaks were obtained, they were processed by mzMine 2.0 software to find peaks in up to 300 chromatograms. Relative peak intensities of both identified and unidentified metabolites were generated and used for further analyses.

Data analysis

Before differential expression analysis, genes with no counts were removed from analysis to avoid skewing the results. A gene was considered significantly differentially expressed if the false discovery rate (FDR) corrected p-value was <.1. Standard Euclidean hierarchical clustering was performed on all the gene expression data and clustering was done by each individual gene. For the metabolomics dataset, primary metabolites, lipids, and biogenic amines were read in and combined into one data frame to analyze the data similarly. Experimental replicates were averaged together and fold changes were calculated from the cell samples and the fresh media samples. Significance of metabolite differences was determined with a p-value <.05 using the Mann-Whitney U test. All statistical analyses were performed using R version 3.4.0.

Gene enrichment analysis

To perform gene enrichment analysis, the Database for Annotation, Visualization, and Integrated Discovery (DAVID) Bioinformatics Resource was used with a list of DEGs for each compound (Huang ,b). The Functional Annotation Tool was used to determine which Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were overrepresented, or enriched. Entrez gene IDs were submitted to the DAVID Bioinformatics Resource website and the Rattus norvegicus species was selected. The category “KEGG pathways” and functional annotation clustering were selected. KEGG pathway terms were considered significantly enriched if the FDR corrected p-value was <.1.

Flux balance analysis and the creation of tissue-specific models

The stoichiometric matrix (S matrix) was analyzed using the COBRA toolbox v. 2.0.6 (Schellenberger ). The iRno reconstruction of rat metabolism, which accounts for the function of 5620 metabolites, 2324 genes, and 8268 reactions, was used to make computational predictions (Blais ). iRno has been curated to perform liver-specific metabolic tasks, making it appropriate as a base model of liver metabolism. Flux balance analysis was performed using the optimizeCBmodel function in the COBRA toolbox in MATLAB v. R2016b. Condition-specific models were then created using the iMAT algorithm in the COBRA toolbox. The createTissueSpecificModel function in the COBRA toolbox was used, with iMAT set as the method for expression data integration, using reactions associated with DEGs and exchange reactions as high confidence reactions to include in the model. Log-fold changes for DEGs were supplied as inputs along with a model with genes created for exchange reactions, whereas the hepatocyte-specific model was provided as an output.

TIMBR algorithm

The TIMBR algorithm combines the transcriptomics data with the iRno network reconstruction to determine production scores for each exchangeable metabolite relative to a control as previously described (Blais ). The transcriptomics data was used to generate weights for a control case and a treatment case on each reaction in the reconstruction. Next, for each metabolite, the weighted flux through each reaction was minimized whereas maintaining positive flux through that metabolite’s exchange reaction for the control and treatment conditions. Production scores are normalized using the previously described formula (Blais ) to determine whether a metabolite has increased or decreased production relative to the control and used for further downstream analysis. The scripts used to generate each of the datasets can be found on the github site (www.github.com/csbl) published with the TIMBR algorithm.

RESULTS

Transcriptomics Data Reveal Compound-Specific Responses of Hepatocytes

Hepatocytes were exposed to APAP, CCl4, TCDD, or TCE for six hours to characterize the differential toxicity-induced metabolic response. Figure 1 shows the experimental layout; after hepatocytes were exposed to each compound, supernatants were collected for metabolomics analyses and RNA was isolated for transcriptomics analysis. DMSO was used as a nondrug control. The number of DEGs for each condition and time point were determined (Table 1) and a list of genes from the differential gene analysis was produced (Supplementary Data 1). APAP induced the most DEGs in the hepatocytes, whereas TCE induced the least number of DEGs. To further analyze the genes that were differentially expressed, we used the DAVID Bioinformatics platform to identify enriched KEGG pathways for each compound. Figure 2A shows the enrichment results of the DEGs for APAP, CCl4, TCDD, and TCE (Complete enrichment results are shown in Supplementary Data 2). APAP at six hours showed an enrichment for metabolic pathways, whereas CCl4, TCDD, and TCE at six hours did not (Figure 2), suggesting that the hepatocyte’s metabolism was more altered globally in response to APAP compared with the other three compounds. As evidenced in the enrichment analysis, APAP exposure induced a wide variety of gene expression changes, whereas gene expression changes after CCl4, TCDD, and TCE exposure appeared focused toward RNA and protein processing. After investigating the broad effects of the compounds, we then focused on metabolic genes to evaluate how each compound perturbed hepatocyte metabolism.

Table 1.

Comparison of the Number of Differentially Expressed Genes in Response to Each Chemical Compound at the six hour Time Point, Compared With Their Respective Controls

Chemical Compound	Number of Differentially Expressed Genes (FDR <0.1)	Number of Differentially Expressed Genes (FDR <0.1) in the iRno Model
APAP – 6 h	7370	1009
CCl4 – 6 h	824	131
TCDD – 6 h	2493	304
TCE – 6 h	907	151

Figure 2.

Gene enrichment and metabolic gene expression data. A, DAVID enrichment of KEGG Pathways for six hours in APAP-, CCl4-, TCDD-, and TCE-induced toxicity conditions. The heat map above shows the log2 fold changes of the metabolic genes from sequencing (B). Each condition is listed on the x-axis, and the individual genes are listed on the y-axis. Genes that are upregulated are shown in red, whereas downregulated genes are shown in blue. Genes on the x-axis are clustered by Euclidean distance, using complete linkage.

Figure 2B shows a heat map of the log2-fold changes of all the metabolic DEGs with a Benjamini-Hochberg adjusted p-value of <.1 in at least one condition. This heatmap shows that CCl4 and TCE elicit similar changes in gene expression. We observed changes in expression for the Cyp450 family of genes, often associated with metabolizing drugs (Guengerich, 2008). We saw a decrease in Cyp3a4 in APAP (Supplementary Figure 3) but no changes in the other compounds, likely because other Cyp450 genes play a role in rat metabolism of compounds (Tran ; Zuber ). Specifically, Cyp2e1 is induced in hepatotoxicity (Jaeschke ; McGill ). We saw upregulation of Cyp2e1 in APAP-induced toxicity, but not for TCDD- and TCE-induced toxicity; however, there were other genes in the Cyp450 family that were differentially expressed in these other conditions. In APAP- and TCDD-induced toxicity, the Cyp450 gene Cyp2d4, also associated with the metabolism of drugs (Mizuno ), was upregulated. TCDD-induced toxicity resulted in upregulation of most other Cyp450-related genes, whereas TCE-induced toxicity resulted in downregulation for many of the same genes. This result highlights that even though there are common pathways of toxicity associated with the liver, these compounds ultimately result in different specific effects on the hepatocytes. In an effort to identify potential biomarkers specific to each compound, we next interrogated the metabolomics data to identify differential effects of each compound (Table 2).

Table 2.

Comparison of the Number of Differentially Changed Metabolites for Each Chemical Compound at the 6-h Time Point, Compared With Their Respective Controls

Chemical Compound	Statistically Significant Metabolites Changed	Subset of Statistically Changed Metabolites Identified
APAP – 6 h	82	8
CCl4 – 6 h	11	3
TCDD – 6 h	102	6
TCE – 6 h	84	6

Metabolomic Data Discriminates the Response of the Primary Hepatocytes Specific to Each Treatment

APAP produces the most distinct signature of the three compounds, whereas TCDD and TCE display a similar profile. The metabolomics data are illustrated in scatter plots for APAP (Figure 3A), CCl4 (Figure 3B), TCDD (Figure 3C), and TCE (Figure 3D) exposure conditions and fold changes with respect to the DMSO control is described in Supplementary Data 3. The scatter plots show each metabolite, with the fold change of average relative metabolite peak intensity compared with blank medium on the x-axis, and compared with the DMSO controls on the y-axis. With this arrangement, metabolites are classified as having increased or decreased production if the fold change relative to blank is positive, or increased or decreased consumption if the fold change relative to blank is negative. Metabolites are also color coordinated, to help distinguish metabolites that were increased or decreased in their production or consumption. Only metabolites that were significantly changed in either the treated versus control, or treated versus blank cases are displayed. From these data, we see that APAP induces the greatest number of metabolites with an increase in production, whereas the other compounds induced a decrease in production of most measured metabolites. There is a trend for metabolites to either be increased in production (upper right) or decreased in production (lower right). This trend is clear in each condition, as these were the two categories with the most metabolites, although many of these metabolites are not yet identified. In APAP-induced toxicity, there were several amino acids that decreased in production compared with the control case (Supplementary Figure 4A, bottom left and right). This result indicates that hepatocytes consumed more amino acids after being exposed to APAP. TCDD and TCE both caused hepatocytes to decrease production of fatty acids (Supplementary Figs. 4C and 4D, bottom right), whereas APAP triggered an increased production of fatty acids (Supplementary Figure 4A, top right). The results from the metabolomics data suggest a clear metabolic difference in the hepatocytes treated with different compounds, and that the mechanism of action or off target effects of the toxicants may be the likely cause of this shift.

Figure 3.

Overview of the metabolomics data. The scatter plots show the distribution of metabolites that are significantly (p<.05) changed when compared with either the control media or blank media, and colored according to their levels when compared with both sets of media. Metabolites in the top left corner have decreased overall consumption, metabolites in the bottom left corner have increased overall consumption, the bottom right corner indicates decreased overall production, whereas the top right corner shows increased overall production, all with respect to the control media. Plots are shown for APAP- (A), CCl4- (B), TCDD- (C), and TCE- (D) induced toxicity conditions at six hours. The heat map above shows the log2 fold changes for metabolites compared with their respective controls (E). Each condition is listed on the x-axis, and the metabolites are listed on the y-axis. Metabolites on the x-axis are clustered by Euclidean distance, using complete linkage.

We next decided to interrogate the total metabolic response of the hepatocytes to further discriminate treatment conditions. Figure 3E shows a heatmap of the individual metabolite levels, and whether or not the amount of the metabolites increased or decreased with respect to the control condition. Again, we noticed that TCE and TCDD showed a similar but distinct pattern of changes in metabolite levels. Valine and leucine were uniquely increased in TCE, whereas tryptophan, serine, and glutamate were uniquely decreased in TCDD. Between both compounds, nicotinate, glucose-1-phosphate, and aminomalonate all decreased. There were only seven metabolites that increased for both TCDD and TCE, 1,3-diheptadecanoyl-2-(10Z-heptadecenoyl)-glycerol d5 and six unidentified metabolites. There were 80 metabolites that decreased between both compounds including both identified and unidentified metabolites. In the heatmap in Figure 3E there are a few prominent clusters of metabolites. There was a small cluster of unidentified metabolites in the TCDD and TCE condition whose levels were decreased when compared with the control condition. APAP did not follow this trend, as a number of those same metabolites were increased. Within this large cluster the only identified metabolite was nicotinate. Of the 559 metabolites we were able to detect, only 115 could be identified. Of the identified metabolites, we then looked at the unique metabolites altered by each condition to compare and contrast each compound’s effect on the hepatocytes. Common metabolites that consistently decreased across all conditions were l-lactate, glycerate, and alpha-ketoglutarate (AKG), which have been shown to decrease in other toxicity studies (Kim and Moon, 2012). Other studies have shown decreases in citrate and AKG (Ishihara ), which have been attributed to disruptions of the TCA cycle. Finally, there were increased lipid levels in TCE and TCDD compared with their controls, suggesting a strong alteration in lipid metabolism in response to these compounds.

TIMBR Predictions Suggest Unique Responses to Each Toxicant

To make predictions on metabolite production levels relative to control from the gene expression data, we created a hepatocyte-specific metabolic model from the unconstrained iRno GENRE using iMAT (Zur ) along with the gene expression data described earlier. The Transcriptionally Inferred Metabolic Biomarker Response (TIMBR) response algorithm (Blais ) was used to create normalized production scores for each metabolite that could be secreted by the model and we compared these values with the fold changes calculated from the metabolomics data above (Supplementary Data 4). Figure 4A shows a distribution of normalized TIMBR production scores by compound, with the median indicated by the notches and black line and the mean represented by the white diamond in the middle of the box plot. In Figure 4A, we see that each group has its mean at about zero, however, the median for each group is different. The analysis of the APAP and TCE conditions show that more scores have positive TIMBR scores while CCl4 has slightly more negative TIMBR scores. This result suggests that hepatocytes are predicted to produce more metabolites in response to APAP and TCE exposure compared with other toxicant conditions.

Figure 4.

Summary and Distribution of TIMBR production scores. The distribution of TIMBR production scores are shown (A) indicating that the ranges are similar, but scores have a slight skew according to their condition. The APAP condition results in more negative production scores, whereas TCE results in more positive production scores. Lines mark y = 1 and y = −1. Venn diagrams compare all positive (B) and negative (C) production scores for each compound, and the overlap between the 3 conditions. TIMBR scores that are common across all conditions (D), and unique to APAP (E), CCl4 (F), TCDD (G), TCE (H) are illustrated. Here, metabolites are classified into categories taken from the subclass names from the Human Metabolome DataBase (HMDB) if available. Metabolite in a category that increases were given a light color, whereas metabolites in a category that decrease were given a darker color.

We then compared common metabolites that were predicted to increase or decrease after all treatments, which indicate common metabolic shifts in response to drug treatment. Figures 4B and 4C shows Venn diagrams of the metabolites that were predicted to commonly increase or decrease in production, or uniquely increase or decrease in response to APAP, CCl4, TCDD, or TCE at six hours, respectively. Similar to the trend noted in all of the TIMBR production scores (Figure 4A), CCl4 exposure was predicted to decrease a higher number of metabolite production scores (Figure 4C) that do not also decrease in other conditions. However, APAP exposure was predicted to result in the increase of more metabolite production scores that were not increased in other conditions (Figure 4B), which is consistent with the prediction of more positive production scores. We then classified metabolites according to their Human Metabolome Database (HMDB) subclassification (Supplementary Data 5) that were changing in each condition from the results shown in Figures 4B and 4C. The bar charts in Figures 4D−H indicate the number of metabolites uniquely predicted to increase or decrease production after toxicant exposure according to their subclassification. For shared metabolites across all conditions (Figure 4D), a small number of amino acids are predicted to decrease, whereas bile acids are predicted to increase. APAP exposure (Figure 4E) resulted in the highest number of fatty acids predicted to increase in production, followed by amino acids. The increase in bile acids and amino acids suggests alterations in these pathways in response to liver injury, and has been observed in literature (Kumar ; Sun ). In CCl4 exposure (Figure 4F), carbohydrate compounds are predicted to decrease in production, whereas these same metabolites were predicted to increase in the other three conditions. With TCDD exposure (Figure 4G), amino acids are predicted to decrease in production while with TCE exposure bile acids are predicted to decrease (Figure 4H) which is similar to CCl4. Overall the TIMBR predictions illustrate that the response of the hepatocytes to each compound is primarily due to carbohydrate and amino acid metabolism, which could represent a generic response toward toxic compounds. However, predictions from APAP exposure indicate a distinct response in fatty acid metabolism, with CCl4 and TCE eliciting more of a change in bile acid metabolism, suggesting that we can predict compound-specific effects on the hepatocytes.

Comparing TIMBR Predictions and Metabolomics Data

We next wanted to quantify the similarity and dissimilarity between the predictions and metabolomics data, to determine how indicative gene expression changes were in predicting metabolite levels. In addition to the fold changes calculated from the metabolomics data, the Mann-Whitney U test was used to determine statistical significance at the p < .05 level. A change in metabolite levels with p > .05 when compared with the control condition, was classified as “no change,” and represented with a fold change of zero. We then took the subset of secreted, identified metabolites and compared this list with our TIMBR predictions which resulted in 20 metabolites we could validate for each experimental condition. For the TIMBR predictions, metabolite production scores were ranked, and metabolites in the middle 50% of the list were classified as no change and given a value of zero for comparing with the metabolomics data. Figure 5A shows a heatmap of the metabolomics data and the production scores for metabolites on the y-axis with each condition on the x-axis. Of the 20 metabolites in each condition, we predicted five correctly in the APAP condition, 10 correctly in the CCl4 condition, 9 correctly in the TCDD condition, and 9 correctly in the TCE condition. From the list of successful predictions nicotinate, and glycine were correct in three of the conditions, whereas 9 metabolites were correct in two of the conditions. There were six amino acids in the set of 20 that we could make predictions for, and of those six, we correctly predicted 2 in the APAP condition, three in the CCl4 condition, whereas only one prediction was correct in the TCDD condition and five in the TCE condition. Although we were able to predict broad changes in carbohydrate and energy metabolism from the TIMBR predictions as described above, the data were too limited to draw the same conclusions from the subset of experimental data that we were able to validate.

Figure 5.

Validation of TIMBR production scores using metabolomics data. The heat map (A) shows the results from the metabolomics data, and the TIMBR production scores for each metabolite we were able to make a prediction for and validate. Each condition is listed on the x-axis, and the metabolites are listed on the y-axis. Metabolomics data are shown in the upper left triangle, and TIMBR production scores are shown in the bottom right triangle. The bar chart (B) shows the categories a prediction can fall into on the y-axis ranging from increase, decrease, or no change for both the experimental data and the TIMBR predictions. The x-axis contains the number of predictions that fall into the category on the y-axis. Predictions that agree with the experimental data have dark colored bars, whereas disagreement between the data show light colored bars.

Figure 5B quantifies our validation results, and shows exactly where predictions were right and where predictions were wrong. The bulk of the correct predictions came from identifying no change in both the experimental condition and the computational prediction. Overall, our accuracy for our predictions was 41%. We made no correct predictions on metabolites that were measured as increased or predicted to increase. Thirty-seven of the predictions were incorrect from detecting a change and predicting there was none, or vice versa. Our sensitivity for detecting no change was 42% and lower for predicting an increase (0%) or decrease (40%). Our specificity for no change or decrease was high at 75%, but was lower (65%) for the no change condition.

DISCUSSION

There is limited information on biomarkers of toxicity; therefore, novel approaches to elucidate and validate relevant biomarkers are needed. A promise of metabolomics as an approach for identifying biomarkers is its connection to cell phenotype as a change in metabolite levels may represent changes in the functional state of the cells. Here, we present the first use of paired transcriptomics and metabolomics with GENREs to study hepatocytes exposed to different compounds and to integrate these data with metabolic network models to provide insight into the changes that are occurring. At the concentrations and timepoint we selected for exposure, we used standard measures of hepatocyte function (albumin production, Cyp450 activity, etc.) to ensure we were not killing the cells (Supplementary Figs. 1–3). Although we observed minimal changes in traditional measurements of toxicity after six hour exposure (Supplementary Figs. 1 − 3), we observed changes in metabolism as indicated by the transcriptomics and metabolomics data. Additionally, connecting transcriptomic changes to secreted metabolites even at low-toxic compound concentration can be useful in clinical settings, as these secreted metabolites can be measured to gain an early indication of hepatic injury. Secreted metabolites can then be connected back to transcriptional changes using metabolic network models, which allows us to generate mechanistic insight into observed changes in transcript or metabolite levels.

From this study, we observed a number of transcriptional changes in metabolic genes of hepatocytes following exposure to compounds. Analyses of transcriptional changes highlighted that APAP produced the largest change in hepatocyte gene expression, as expected (Ben-Shachar ; McGill and Jaeschke, 2013; Sjogren ; Taguchi ). There were 31 differentially expressed metabolic genes that were shared across all four compounds. There were five genes that were upregulated in each of the four treatment conditions. Among the group of upregulated genes includes two glutathione S-transferase genes, indicative of detoxification mechanisms given glutathione is used to conjugate toxic metabolites (Guengerich, 2008; Monks ). Furthermore, glutathione S-transferase is responsible for the detoxification of NAPQI, a toxic metabolite that is generated from metabolizing APAP (Henderson ). In the metabolomics data, glutathione production was decreased in APAP (albeit p = .34). For the TIMBR predictions in the APAP condition, we did predict decreased production of glutathione, which is attributed toward glutathione detoxifying NAPQI in the hepatocytes.

Twelve differentially expressed metabolic genes that were shared across all compounds were downregulated. Among this group was isocitrate dehydrogenase 3 (Idh3a), which is responsible for the NAD+ dependent conversion of isocitrate to AKG. In the metabolomics data, we observed decreased production of AKG in response to APAP, TCDD, and TCE compared with their respective controls. We also computationally predicted this decrease in AKG in the APAP condition (Figure 5). These examples suggest that there are some transcriptional changes that are indicative of downstream metabolite changes.

Glycolysis and the TCA cycle were disrupted as a result of compound exposure. In the metabolomics data, we observed that glycerate was decreased in response to exposure to APAP, TCDD and TCE, and glucose-1-phosphate was decreased after treatment with TCDD and TCE. Both glucose-1-phosphate and glycerate can feed into glycolysis and then progress to the TCA cycle. Decreases in these metabolites indicate that the hepatocytes are inefficiently producing ATP via the TCA cycle. This observation is also further supported by the measured decrease in AKG in most of the conditions as well as the decrease in citrate in response to APAP. Carbohydrates also feed into glycolysis, and a decrease in carbohydrates can also decrease TCA activity. In Figure 4, we observe that carbohydrates were predicted to decrease in production after exposure to APAP and CCl4, which are both hepatotoxicants. Although we were not able to correctly predict changes in glycerate production in every condition (Figure 5), we were able to predict this shift in metabolism via the carbohydrates, which is supported by the metabolomics data. Thus, TIMBR predictions can be useful for suggesting pathway level differences of a treatment that can be experimentally validated.

We compared our in vitro and computational results with other in vivo toxicity studies that have been done. Across the different studies, lipid metabolism, amino acid metabolism, and energy metabolism (TCA Cycle) were all affected by exposure to different compounds. One study that focused on TCDD-induced transcriptomic changes identified several genes associated with these pathways that were both upregulated and down regulated (Boverhof ). From a metabolomics perspective, TCA cycle intermediates were down regulated in response to APAP (Sun ), which agreed with our data. These same pathways came up in common with our TIMBR predictions (Figure 4) that are based on our measured transcriptional changes. One study noted that in response to APAP-induced toxicity, metabolite levels for glycerol and kynurenine were increased, whereas threonine, serine, ornithine, lysine, glycerate, and glutathione were reduced (Pannala ). The authors also observed enrichment in the glycine, serine, and threonine metabolic pathway (Pannala ). Although we did not observe the decrease in glutathione levels, we did note enrichment in the glycine, serine, and threonine pathway in the APAP condition (Figure 2A). The decrease in glutathione in APAP was shown to occur at later time points, due to increasing progression of liver injury as noted by the authors (Pannala ). Lastly, CCl4 is known to cause hepatocytes to increase urinary bile acid levels (Yang ). We observed that there were a few bile acids predicted to increase (Figure 4D), but unique to CCl4 was the observation that most of the bile acids were predicted to decrease (Figure 4F). Because in vitro conditions do not fully capture in vivo conditions due to differences in time-scales, actual exposure concentrations, among other variables, there is not complete agreement between the in vitro and in vivo results as expected. However, our in vitro experiment provides a means to study changes in hepatocyte metabolism without the variability of an in vivo experiment. Although results may not fully match, general trends in metabolic changes do agree, as indicated by the shift in fatty acid metabolism from TCA cycle and amino acid breakdown noted earlier, which highlights the utility of in vitro systems for interrogating toxicological responses.

One limitation of this study that affected the ability to make predictions was the lack of overlap between the metabolomics data and metabolites for which we were able to make TIMBR predictions. For the primary metabolites in the metabolomics dataset, only 115 out of 559 were identifiable. Of these 115, there were only 21 metabolites in the subset that were secreted and that were accounted for in our current network reconstruction, as shown in Figure 5. Although the number of correct predictions was limited, we were still able to make predictions on glycolysis, the TCA cycle, and amino acid metabolism that were supported by either the metabolomics data or literature from other toxicity studies (Beger ; Kumar ). There are opportunities for further curation of the network reconstruction to account for more metabolites and metabolic reactions, as well as further curation of the metabolomics data.

This study used transcriptomics data paired with metabolomics data to provide insight into the changes induced by these toxicants on hepatocytes. Protein fold changes could be used in place of gene expression data and ultimately could have been used for TIMBR predictions because we can map such data to the metabolic reactions accounted for in the metabolic network reconstruction. As large data sets are made accessible or easy to collect, the use of multi-omic datasets to predict and validate modeling results becomes critical in interrogating specific phenotypes of interest for a chosen system. Our paired experimental and computational approach is one step toward characterizing the cellular response to a compound and identifying potential biomarkers indicative of cell state.

DATA AVAILABILITY

The sequencing and processed data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (Edgar ) with the accession number GSE129814. All supplementary files are stored on the Dryad Digital Repository website (Rawls ), accessible at https://doi.org/10.5061/dryad.04vk390. Source data needed to reproduce figures can be obtained via the code available at github.com/csbl/hepatocyte_omicsdata. All other data supporting the findings of this study are available within the article and its supplementary information files.

SUPPLEMENTARY DATA

Supplementary data are available at Toxicological Sciences online.

DECLARATION OF CONFLICTING INTERESTS

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Click here for additional data file.