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

Abstract

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.

Keywords: constraint-based modeling; metabolism; metabolomics; toxicology; transcriptomics.

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.

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 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.

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), CCl₄ (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.

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.