2022 Volume 47 Issue 8 Pages 337-348
Drug-induced liver injury (DILI) is the main cause of failure in drug development and postapproval withdrawal. Although toxicogenomic techniques provide an unprecedented opportunity for mechanistic assessment and biomarker discovery, they are not suitable for the screening of large numbers of exploratory compounds in early drug discovery. Using a comprehensive analysis of toxicogenomics (TGx) data, we aimed to find DILI-relevant transcription factors (TFs) that could be incorporated into a reporter gene assay system. Gene set enrichment analysis (GSEA) of the Open TG-GATEs dataset highlighted 4 DILI-relevant TFs, including CREB, NRF2, ELK-1, and E2F. Using ten drugs with already assigned idiosyncratic toxicity (IDT) risks, reporter gene assays were conducted in HepG2 cells in the presence of the S9 mix. There were weak correlations between NRF2 activity and IDT risk, whereas strong correlations were observed between CREB activity and IDT risk. In addition, CREB activation associated with 3 Withdrawn/Black box Warning drugs was reversed by pretreatment with a PKA inhibitor. Collectively, we suggest that CREB might be a sensitive biomarker for DILI prediction, and its response to stress induced by high-risk drugs might be primarily regulated by the PKA/CREB signaling pathway.
Drug-induced liver injury (DILI), one of the most commonly recognized causes of acute liver injury, is the main reason for drug failure in the drug development phase and the withdrawal of drugs from the market (David and Hamilton, 2010; Stravitz and Lee, 2019; Temple and Himmel, 2002). Despite efforts for DILI prediction, it remains a worrisome and tough issue because at least 6 mechanisms are involved in DILI, including cell injury/death, oxidative stress, mitochondrial dysfunction, and accumulation of fatty acids, most of which are specific, unpredictable and dose-independent according to the pharmacological effects of drugs (Lee, 2003; Arnesdotter et al., 2021; Katarey and Verma, 2016). Conventional animal studies for the assessment of human hepatotoxicity have obstacles such as cross-species extrapolation, resulting in low concordance and poor prediction of human hepatotoxicity (Olson et al., 2000). Researchers are still looking for reliable biomarkers and cell culture methods to improve the accuracy of DILI prediction through in vitro methods. Therefore, there is a pressing need to establish a high-throughput, sensitive, and robust DILI risk prediction system to improve the detection rate and efficiency of drug hepatotoxicity, thereby reducing the cost of preclinical drug development and improving the safety of clinical medication.
Toxicogenomics (TGx), which aims to apply comprehensive technologies to investigate the adverse effects of environmental and pharmaceutical chemicals on human health, provides an unprecedented opportunity for the clarification of toxicity mechanisms and biomarker discovery in drug risk assessment (Liu et al., 2019). Toxicogenomics databases, such as DrugMatrix (Ganter et al., 2005) in the USA, PredTox (Suter et al., 2011) in Europe, and TG-GATEs in Japan (Igarashi et al., 2015) have been published with global efforts to integrate the use of toxicogenomic data into DILI risk assessment. These databases have facilitated the deeper understanding of DILI mechanisms through the interpretation of archived data. For example, Gao et al. (2010) extracted important genes (AKR7A3, TRIB3, and GSTP1) associated with hepatotoxicity of chemical-induced glutathione depletion from rat data in TG-GATEs. In addition, studies using TGx have made progress in hepatotoxicity prediction. A research team used TGx to establish a model for predicting drug hepatotoxicity in rats from the chemical structure with a correct classification rate of approximately 70% (Low et al., 2011). In addition to these efforts, the US Tox21 in vitro high-throughput transcriptomic screening program (Tox21) combined next-generation sequencing (NGS) technology to build a compound library for faster toxicological characterization (Merrick, 2019). This library contained nearly 10 000 compounds in its production phase (Attene-Ramos et al., 2013). Thus, the constantly evolving TGx can provide researchers with rich and reliable datasets for further data mining to advance the development of toxicity assessments.
Nonetheless, there is a strong need for high-throughput toxicity assessment of large numbers of exploratory compounds in the early stages of drug discovery. Changes in gene expression are fundamentally related to the activity of transcription factors (TFs). We thought that identification of DILI-relevant TFs from the TGx data, would allow their incorporation into a reporter gene assay system. Reporter gene assays are most commonly used for monitoring the dynamics of gene expression and regulation (Schenborn and Groskreutz, 1999). They provide high sensitivity, robustness, and wide linearity, and the ability to rapidly screen a large number of samples on multiwell plates (Fan and Wood, 2007). This idea was supported by Yuan and Kaplowitz, who stated that the mechanisms by which liver cells defend against drug-induced toxicity are related to the expression of various genes (Yuan and Kaplowitz, 2013). Indeed, pregnenolone X receptor (PXR), constitutive androstane receptor (CAR), and steroid and xenobiotic receptor (SXR) have been found to enhance protection against xenotoxicants by regulating the levels of some cytochrome P450 enzymes (CYPs) (Xie et al., 2000a, 2000b). A lot of evidence of the activation of TFs related to DILI has also been reported, among which the most widely reported was the nuclear factor erythroid 2-related factor 2 (NRF2), which responds to oxidative stress (Ma, 2013; Xu et al., 2008). Others include the signal transducer and activator of transcription 4 (STAT4) (Urban et al., 2012), and activator protein-1 (AP-1) (Hasselblatt et al., 2007). In this context, the present study aimed to identify human DILI-specific TFs from a comprehensive analysis of TGx data. We performed gene set enrichment analysis (GSEA) using a database of transcription factor binding sites (TFBSs) and TG-GATEs data. The usefulness of the TFs identified in the reporter gene assay was also validated.
The drugs in the TG-GATEs repository were annotated according to the Liver Toxicity Knowledge Base Benchmark Dataset (LTKB-BD) (https://www.fda.gov/science-research/liver-toxicity-knowledge-base-ltkb/ltkb-benchmark-dataset), which is a list of DILI-concern drugs classified into 3 levels based on FDA-approved prescription drug labels. Microarray data for primary human hepatocytes exposed to drugs registered in both LTKB-BD and TG-GATEs were downloaded from the TG-GATEs website (http://toxico.nibiohn.go.jp/english/) for further analysis.
Correlation of DILI risk and relative expression change for each geneThe raw data for each of the different exposure levels (high, middle) and incubation times (8 hr, 24 hr) were subjected to analysis in R project 4.0.4 (https://www.r-project.org/). First, the downloaded data for each drug were normalized using the robust multiarray average (RMA) algorithm in the “Affy” package (Gautier et al., 2004). Second, the logarithm of the fold change (logFC) in gene expression from untreated controls to treated samples was calculated for each treatment condition. We then divided drugs into 2 groups; less and no-DILI-concern drugs were classified in the low-risk group, whereas most-DILI-concern drugs were classified in the high-risk group. Consecutively, we performed a logistic regression analysis of the modified DILI risks against logFC values of each probe. The correlation for each probe was expressed using Cox-Snell R2 values and annotated using the HGU133Plus2.db package in 20 161 genes. For multiprobe genes, R values were integrated into a single value using Fisher’s Z transformation.
Gene set enrichment analysis (GSEA)For GSEA analysis, genes were first sorted by the R value, which represents the correlation between the DILI risk of a drug and the drug-induced change in gene expression. GSEA analysis was then conducted using the GSEA Pre-ranked tool in the GSEA software 4.1.0 (http://www.broadinstitute.org/gsea/). The TFT_Legacy subset of TFT (transcription factor targets) was selected from the C3 Collection of the Molecular Signatures Database (MSigDB v7.4), which contains 610 transcription factor binding site (TFBS) gene sets. Genes were ranked by R values, and significant TFBSs were identified based on the false discovery rate (FDR) q-value. The permutation test was conducted 1000 times to yield robust q-values.
Chemicals and reagentsAcetaminophen, tert-butylhydroquinone, and tienilic acid were purchased from Sigma-Aldrich (St. Louis, MO, USA). Aminopyrine, sulfamethoxazole, and forskolin were purchased from Fujifilm Wako Pure Chemical Co. (Osaka, Japan). Erythromycin, valsartan, and warfarin were purchased from LKT Laboratories, Inc. (St Paul, MN, USA). Nevirapine, pravastatin, and procainamide were purchased from Tokyo Chemical Industry Co., Ltd. (Tokyo, Japan). Stock solutions of all test drugs were prepared in dimethyl sulfoxide (DMSO) at a concentration of 200 mM. Each solution was adequately diluted with Dulbecco’s modified Eagle’s medium (DMEM) (Nissui Pharmaceutical Co., Ltd., Tokyo, Japan) so that the final concentration of DMSO in the culture medium was 0.5%.
Luciferase reporter gene assayHepG2 human hepatoma cells were obtained from the National Institute of Biomedical Innovation (Osaka, Japan). Cells were maintained in DMEM supplemented with 10% fetal bovine serum (FBS) (Thermo Fisher Scientific Inc., Waltham, MA, USA) and penicillin-streptomycin-glutamine (Thermo Fisher Scientific Inc.) at 37ºC in an atmosphere of 5% CO2 in air.
HepG2 cells were seeded at a density of 5.0 × 104 cells/100 μL in a 96-well plate. After 24 hr in culture, 5 μL of Opti-MEM I reduced serum medium (Thermo Fisher Scientific Inc.) containing 20 mg/mL NRF2 reporter vector (pGL4.37[luc2P/ARE/Hygro]; Promega, Madison, WI, USA) or CREB reporter vector (pGL4.29[luc2P/CRE/Hygro]; Promega), 0.002 mg/mL control vector (pNL1.1.PGK [Nluc/PGK]; Promega) and 9% FuGENE ®HD (Promega) was added to each well for transfection. After 4 hr of transfection, transfection reagents were removed by replacing the medium with fresh one. After a subsequent 20 hr incubation, the medium was replaced with freshly prepared medium containing 1 mg equivalent/mL S9 mix, which consisted of the human S9 fraction (Corning® Gentest™ Pooled human liver S9, Corning, Corning, NY, USA) and cofactor I solution (Oriental Yeast Co., Ltd., Tokyo, Japan). The cofactor reagent was reconstituted with 2.34 mL purified water and filtered through a 0.45 µm Millipore membrane filter (Merck, Darmstadt, Germany). The S9 mix was prepared by mixing S9 and cofactor solutions at a ratio of 3:7. After incubation of cells with the S9 mix and each test drug (100 or 1000 µM) for 8 hr, luciferase activity was measured using the Nano-Glo® Dual-Luciferase® Reporter Assay System (Promega) on a Wallac 1420 ARVOTM MX (PerkinElmer, Waltham, MA, USA). The luciferase activity assay was performed according to the manufacturer’s instructions. Firefly luciferase activity was normalized to NanoLuc luciferase activity. Data are represented as fold induction of the drug-treated group to the control group. The polyserial correlation coefficients between luciferase activity (firefly/NanoLuc luciferase activity) and idiosyncratic drug toxicity (IDT) risks, which range from −1 to 1, were calculated using the “polycor” package in R.
For the protein kinase A (PKA) inhibition study, cells were treated with 50 µM H-89 dihydrochloride hydrate (H-89, Sigma-Aldrich) 1 hr prior to the addition of S9 mix and test drugs.
Statistical analysisResults of the reporter gene assay are expressed as the mean ± standard deviation (SD). Unpaired t-tests were used to analyze the significance between H-89-treated and H-89-untreated groups.
In this study, we investigated DILI-relevant TFs using microarray data of primary human hepatocytes exposed to drugs from the Open TG-GATEs. First, we set up a list of 102 drugs that met the criteria of being registered in both TG-GATEs and LTKB-BD. These drugs were classified into 49 most-, 46 less-, and 7 no-DILI-concern drugs according to the classification by the FDA (Supplementary Table S1). As the number of no-DILI-concern drugs was too small, we combined less- and no-DILI-concern drugs into one low-risk group and compared them to most-DILI-concern drugs, which represented the high-risk group.
The DILI risk of a drug is thus a binary variable, whereas the change in gene expression (i.e., logFC) due to drug treatment is a quantitative variable. Therefore, we expressed the correlation between the 2 variables in terms of Cox-Snell R2 values, with the positive or negative values of their square roots being determined from the slope. Furthermore, we averaged the R values of the probes for each gene using Fisher’s Z-transformation. Figure 1 represents the distribution of the R value with DILI risk for a total of 20 161 genes for each of the 4 treatment conditions (high dose-8 hr, high dose-24 hr, middle dose-8 hr, and middle dose-24 hr). We found that the number of genes with an absolute R > 0.25 in the 4 treatment conditions were 354 (high dose, 8 hr), 199 (high dose, 24 hr), 487 (middle dose, 8 hr), and 108 (middle dose, 24 hr). However, no genes had an absolute R of > 0.25 in all treatment conditions.

Box-whisker plots of R values between logFC and DILI risks for each gene in 4 treatment conditions. A positive R value indicates the upregulation of genes by DILI-concern drugs, whereas a negative R value indicates the opposite.
As is known, TFs regulate the expression of multiple genes. Accordingly, if a group of genes shows the same tendency, it would be meaningful even if the correlation with DILI was weak for individual genes. Therefore, we performed GSEA using the TFBS gene sets to identify TFs associated with DILI. A total of 610 TFBSs were registered in the C3 collection from MSigDB (TFT_Legacy subset of TFT (transcription factor targets)), corresponding to 280 TFs. Table 1 summarizes the number of significant TFBSs with a false discovery rate (FDR) q-value < 0.05. We observed that most normalized enrichment scores (NES) were positive in all 4 treatment conditions, indicating that the majority of statistically significant TFBSs were upregulated by test drugs. Table 2 summarizes 25 TFBSs, that the regulation of which was statistically significant in at least 3 treatment conditions, corresponding to 9 TFs.
| Treatment conditions | NESa) > 0 | NESa) < 0 |
|---|---|---|
| High dose 8 hr | 9 | 1 |
| High dose 24 hr | 45 | 3 |
| Middle dose 8 hr | 42 | 2 |
| Middle dose 24 hr | 21 | 0 |
a) Normalized enrichment score. A positive value indicates that the TFBS is involved in the DILI-related drug-induced upregulation of downstream genes.
| TF | TFBSa) | High dose | Middle dose | ||||||
|---|---|---|---|---|---|---|---|---|---|
| 8 hr | 24 hr | 8 hr | 24 hr | ||||||
| NES | FDRq-value | NES | FDRq-value | NES | FDRq-value | NES | FDRq-value | ||
| NRF2 | GABP_B | 1.808 | 0.009 | 1.833 | 0.004 | 1.544 | 0.015 | 0.855 | 1.000 |
| NRF2_01 | 2.032 | 0.002 | 1.942 | 0.001 | 1.961 | 0.002 | −0.914 | 1.000 | |
| ELK1 | ELK1_02 | 2.215 | 0.002 | 1.651 | 0.011 | 1.983 | 0.002 | 0.930 | 1.000 |
| CREB | CREB_02 | 1.662 | 0.023 | 1.509 | 0.032 | 1.641 | 0.007 | 0.801 | 1.000 |
| CREB_Q3 | 1.731 | 0.014 | 1.627 | 0.014 | 1.526 | 0.017 | 0.855 | 1.000 | |
| E2F | E2F_01 | 0.997 | 0.809 | 1.525 | 0.028 | 1.508 | 0.019 | 1.818 | 0.002 |
| E2F_02 | 1.076 | 0.565 | 1.976 | 0.002 | 1.719 | 0.004 | 2.281 | < 0.001 | |
| E2F_03 | −0.924 | 1.000 | 1.717 | 0.007 | 1.656 | 0.007 | 2.206 | < 0.001 | |
| E2F_Q3 | 0.837 | 1.000 | 1.689 | 0.008 | 1.879 | 0.002 | 1.990 | < 0.001 | |
| E2F_Q3_01 | 0.875 | 1.000 | 1.886 | 0.002 | 1.953 | 0.002 | 2.149 | < 0.001 | |
| E2F_Q4 | 1.053 | 0.655 | 2.072 | 0.001 | 1.837 | 0.003 | 2.348 | < 0.001 | |
| E2F_Q4_01 | 0.945 | 0.941 | 1.780 | 0.004 | 1.952 | 0.002 | 2.194 | < 0.001 | |
| E2F_Q6 | 1.032 | 0.725 | 2.033 | 0.001 | 1.796 | 0.003 | 2.314 | < 0.001 | |
| E2F_Q6_01 | 0.768 | 1.000 | 1.526 | 0.028 | 1.763 | 0.004 | 1.958 | < 0.001 | |
| E2F-1 | E2F1_Q3 | 1.083 | 0.547 | 2.041 | 0.001 | 1.804 | 0.003 | 2.235 | < 0.001 |
| E2F1_Q4 | 0.890 | 1.000 | 1.729 | 0.007 | 1.405 | 0.043 | 1.668 | 0.007 | |
| E2F1_Q4_01 | 0.924 | 0.951 | 1.882 | 0.002 | 1.937 | 0.001 | 2.191 | < 0.001 | |
| E2F1_Q6 | 1.114 | 0.475 | 1.991 | 0.001 | 1.880 | 0.002 | 2.374 | < 0.001 | |
| E2F1_Q6_01 | −0.790 | 1.000 | 1.826 | 0.004 | 1.663 | 0.007 | 2.104 | < 0.001 | |
| E2F-1:DP-1 | E2F1DP1_01 | 1.087 | 0.544 | 1.993 | 0.001 | 1.718 | 0.004 | 2.279 | < 0.001 |
| Rb: E2F-1:DP-1 | E2F1DP1RB_01 | 0.932 | 0.950 | 1.957 | 0.001 | 1.827 | 0.003 | 2.310 | < 0.001 |
| E2F-1:DP-2 | E2F1DP2_01 | 1.075 | 0.565 | 2.006 | 0.002 | 1.743 | 0.004 | 2.310 | < 0.001 |
| SGCGSSAAA_E2F1DP2_01 | 0.991 | 0.818 | 1.975 | 0.001 | 1.852 | 0.003 | 2.401 | < 0.001 | |
| E2F-4:DP-2 | E2F4DP1_01 | 1.007 | 0.780 | 1.933 | 0.001 | 1.755 | 0.004 | 2.257 | < 0.001 |
| E2F4DP2_01 | 1.075 | 0.561 | 1.971 | 0.001 | 1.716 | 0.004 | 2.297 | < 0.001 | |
a) TFBSs with FDR q-values < 0.05 for at least 3 of the 4 treatment conditions.
We investigated the feasibility of using these TFs as markers for DILI detection. Among the 9 TFs, luciferase reporter vectors were commercially available for nuclear factor erythroid 2-related factor 2 (NRF2) and cAMP-response element binding protein 1 (CREB1). We transfected HepG2 cells using the pGL4.37[luc2P/ARE/Hygro] and pGL4.29[luc2P/CRE/Hygro] vectors, as well as a control vector. Fujimoto et al. (2010) performed an LDH assay following exposure of rat primary cultured hepatocytes to 42 approved drugs, which were classified into 3 IDT risk categories (Withdrawn/Black box Warning, Warning, Safe) according to the report by Nakayama et al. (2009). In the present study, we used the results of their toxicity experiments to efficiently perform reporter gene assays on a limited scale. A set of experiments was planned to be completed in a single 96-well plate to minimize variability associated with cell culture, gene transfection, and other processing conditions. We selected 10 drugs out of the 42 drugs and set only two test concentrations, 100 and 1000 μM. Selection criteria for this study were (1) drugs that have been assessed for toxicity at the same concentrations as in the study of Fujimoto et al. (2010), (2) drugs that were not toxic at the given concentrations, and (3) drugs not used in TG-GATEs. Supplementary Fig. S1 shows the cell viability of HepG2 cells treated for 8 hr with the test drugs supplemented with the S9 mix. The lactate dehydrogenase (LDH) assay indicated that HepG2 cell viability remained as high as ~87% after treatment with 100 or 1000 μM of the test drugs. Based on the study of Fujimoto et al. (2010), we decided to use the IDT risk instead of DILI risk, which we consider not a major issue since a majority of DILIs are idiosyncratic (Bell and Chalasani, 2009; Verma and Kaplowitz, 2009).
Figure 2 shows the fold induction of NRF2 activity after 8 hr of treatment with each drug in pGL4.37[luc2P/ARE/Hygro] vector-transfected HepG2 cells. Tert-butylhydroquinone was included as a positive control. We found that in the presence of the S9 mix (Fig. 2a), the “Withdrawn/Black box Warning” drug-treated groups exhibited higher luciferase activity, which was also observed in the “Safe” drug-treated groups (e.g., valsartan and warfarin). However, in the absence of the S9 mix (Fig. 2b), no remarkable change in NRF2 activity was observed. We estimated the correlation between the 3-level IDT risks and change in NRF2 activity using a polyserial correlation coefficient, which was designed for ordinal categorical variables. We accordingly found that the polyserial correlation coefficients between luciferase activity (firefly/NanoLuc luciferase activity) and IDT risks were 0.377 and −0.198 at concentrations of 100 and 1000 µM, respectively.

Reporter gene assay for assessment of NRF2 activity. (a) Fold induction of NRF2 activity using a reporter gene assay in HepG2 cells in the presence of the S9 mix. (b) Fold induction of NRF2 activity using a reporter gene assay in HepG2 cells in the absence of the S9 mix. HepG2 cells, transfected with the pGL4.37[luc2P/ARE/Hygro] and pNL1.1.PGK [Nluc/PGK] vectors, were incubated with 100 or 1000 µM of 10 test drugs and 100 µM tert-butylhydroquinone (positive control) in the presence/absence of the S9 mix. The final concentration of human S9 was 1 mg/mL. Each value represents the mean ± standard deviation. N = 3. For the negative control, N = 6.
We also assessed the CREB activity after 8 hr of treatment with test drugs in pGL4.29[luc2P/CRE/Hygro] vector-transfected HepG2 cells, using forskolin as a positive control. Among the groups treated with the S9 mix, we observed a higher luciferase activity in the “Withdrawn/Black box Warning” drug-treated groups compared with other groups, except for the warfarin-treated group (Fig. 3). In addition, we noticed that these groups had a higher fold induction of CREB activity than the groups treated without the S9 mix. We identified that the polyserial correlation coefficients against IDT risks were −0.306 and 0.202 under treatment conditions with 100 and 1000 µM drug with the S9 mix, respectively. However, we observed that if the warfarin-treated group was excluded, the polyserial correlation coefficients were as high as 0.858 and 0.841 at each concentration.

Reporter gene assay for assessment of CREB activity. (a) Fold induction of CREB activity using a reporter gene assay in HepG2 cells in the presence of the S9 mix. (b) Fold induction of CREB activity using a reporter gene assay in HepG2 cells in the absence of the S9 mix. HepG2 cells transfected with the pGL4.29[luc2P/CRE/Hygro] and pNL1.1.PGK [Nluc/PGK] vectors, were incubated with 100 or 1000 µM of 10 test drugs and 100 µM forskolin (positive control) in the presence/absence of the S9 mix. The final concentration of human S9 was 1 mg/mL. Each value represents the mean ± standard deviation. N = 3. For the negative control, N = 6.
Since Protein kinase A (PKA) is known to be a predominant kinase in CREB activation, it activates CREB by phosphorylating its Ser133 (De Cesare et al., 1999). To investigate the involvement of PKA in the activation of CREB, we pretreated HepG2 cells with 50 µM H-89, which is an inhibitor of PKA, for 1 hr. Then, we measured the activity of CREB after treatment with Withdrawn/Black box Warning drugs (aminopyrine, nevirapine, and tienilic acid) and the S9 mix. As shown in Fig. 4, compared with the H-89-untreated groups, the activities of CREB augmented by the 3 Withdrawn/Black box Warning drugs were reversed by pretreatment with H-89.

Effect of H-89 pretreatment on the activation of CREB by Withdrawn/Black box Warning drugs. HepG2 cells were transfected with the pGL4.29[luc2P/CRE/Hygro] and pNL1.1.PGK [Nluc/PGK] vectors, and then cells were preincubated with 50 µM H-89 for 1 hr before treatment with aminopyrine, nevirapine, and tienilic acid, besides the S9 mix. Each value represents the mean ± standard deviation; Unpaired t-test was used to analyze the significance between H-89-treated and H-89-untreated samples; *: p < 0.05, **: p < 0.01 compared with H-89-untreated samples. N = 3.
In the present study, we attempted to identify DILI-associated TFs as new biomarkers and incorporate them into reporter gene assays for exploratory hepatotoxicity prediction. Prior to microarray data analysis, each drug in the TG-GATEs was labeled with DILI risks based on LTKB-BD. LTKB-BD provides a better-balanced categorization of drug DILI potentiality and the most consistent information for DILI research based on FDA-approved drug labels (Chen et al., 2011). It has been used to investigate the characteristics of DILI-concern drugs or validate the adequacy of observations regarding them (Shah et al., 2015; Aleo et al., 2014). However, our results showed that the R value for each gene was not strongly suggestive of an association with DILI. To identify potential TFs related to DILI, we performed GSEA using preranked gene lists based on R values. GSEA is a popular computational method for gaining insights into biological mechanisms by interpreting gene expression data (Subramanian et al., 2005). It focuses on gene sets, which are defined as groups of genes that share common biological functions or regulations based on prior biological knowledge (Subramanian et al., 2005). GSEA might thus help draw new insights into the intrinsic biological processes when the expression of a set of genes is individually modest but mutually associated (Bragde et al., 2020). It can also provide a better performance in a systematic comparison with several other gene set analysis methods (Mathur et al., 2018). In this study, we focused on TFBSs located in the proximal promoter region, considering that transcriptional activation is generally initiated by the specific binding of TFs to the corresponding TFBSs. Association analysis using transcriptome data in TG-GATEs suggested the activation of NRF2, ETS transcription factor ELK1 (ELK1), CREB, and E2F transcription factor (E2F) following exposure to DILI-concern drugs.
Subsequently, we experimentally evaluated the activity of TFs using luciferase reporter vectors to investigate the feasibility of these TFs as markers for idiosyncratic DILI detection. The HepG2 cell line is a well-characterized and widely used human hepatoma cell line for research on drug metabolism and hepatotoxicity in vitro. It offers advantages such as availability, unlimited lifespan, stable phenotype, and ease of handling compared with primary human hepatocytes (Donato et al., 2008). The biggest advantage of using HepG2 cells in this study was the ease of gene transfection, which enabled us to better observe the dynamics of TFs in the presence of drugs with different risk levels. However, HepG2 cells have very limited levels of cytochrome P450 enzymes and some phase II metabolizing activities (Rodríguez-Antona et al., 2002; Westerink and Schoonen, 2007), which play key roles in the metabolic clearance of most drugs as well as in catalyzing the formation of DILI-related reactive metabolites (Feng and He, 2013; Lammert et al., 2010). For this reason, Otto et al. (2008) first reported a method for in vitro drug metabolic hepatotoxicity using HepG2 cells with the addition of S9. They actually confirmed the formation of reactive metabolite in this system. Liver S9 is composed of liver microsome and cytoplasm from the supernatant obtained after centrifugation of liver homogenate (Ames et al., 1975; Brandon et al., 2003). It is a useful tool for drug biotransformation studies because it presents more complete metabolic characteristics with both phase I and phase II activity (Fasinu et al., 2012). This method of combining HepG2 cells with S9 is commonly applied in in vitro drug cytotoxicity research (Boehme et al., 2010; Bhattacharya et al., 2020). In the present study, luciferase activity in the absence of S9 was significantly lower. (Fig. 2b and Fig. 3b). This suggested that the drug-induced activation of NRF2 and CREB might be metabolism-dependent. Potentially, the S9-mediated reactive metabolites of these drugs might be involved in the activation of these 2 TFs.
Drugs with different IDT risk categories can activate NRF2 and CREB to varying degrees. The NF-E2-related factor 2 (NRF2) is mainly regulated by its cytosolic inhibitor Kelch-like ECH-associated protein 1 (KEAP1) in cells. The NRF2/KEAP1 signaling pathway is a protective and adaptive pathway that regulates a wide range of metabolic responses against oxidative stress (Bellezza et al., 2018). Under endogenous and exogenous harmful stimuli, such as xenobiotics and oxidative stress, NRF2 dissociates from KEAP1 and is located in the nucleus where it binds to antioxidant responsive element (ARE) to activate the expression of genes that control diverse functions, including xenobiotic metabolism (GST and NQO1), antioxidative stress (HO-1), autophagy (p62), and cell growth (BMPR1A, IFG1, and JAG1) (Itoh et al., 1999; Jain et al., 2010; Murakami and Motohashi, 2015; Malhotra et al., 2010). Moreover, NRF2 is also activated by the metabolite of acetaminophen (APAP), which is the most widely investigated drug associated with DILI (Copple et al., 2008). The polyserial correlation coefficients between NRF2 activity and IDT risks at 100 and 1000 µM did not show a correlation between IDT risks and NRF2. Hence, whether NRF2 is a potential biomarker for hepatotoxicity research is a matter of debate (Taguchi and Kensler, 2020). Copple et al. (2019) have shown that the chemical-mediated induction of reactive stress of NRF2 might not be a sensitive signal pathway for DILI prediction. Our results further support theirs.
CREB, a cyclic AMP (cAMP) - responsive element (CRE) binding protein, is a TF that activates the transcription of target genes in response to diverse stimuli, including protein kinase A (PKA), mitogen-activated protein kinases (MAPKs), and Ca2+/calmodulin-dependent protein kinases (CaMKs) (De Cesare et al., 1999; Shaywitz and Greenberg, 1999). All these kinases phosphorylate CREB at a specific residue, serine 133 (Ser133), and induce the translocation of cytoplasmic CREB to the nucleus. CREB is known to promote the expression of genes involved in regeneration, proliferation, liver physiology, mitochondrial biogenesis, and stress response (Servillo et al., 2002; Than et al., 2011; Rui, 2014). The polyserial correlation coefficients of CREB between reporter gene assay and IDT risks at 100 and 1000 µM showed a strong correlation, indicating that CREB plays an important role in drug hepatotoxicity defense.
We further investigated the involvement of PKA in CREB activation by Withdrawn/Black box Warning drugs, because PKA is known to be a predominant kinase in CREB activation. As shown in Fig. 4, pretreatment with H-89, a PKA inhibitor, completely reversed CREB activation by aminopyrine, while partly reversed CREB activation by nevirapine and tienilic acid. Yun et al. (2014) reported that PKA inhibition by H-89 pretreatment prior to oral administration of acetaminophen potentiated acetaminophen-induced hepatotoxicity, as determined by serum alanine transaminase and glutathione depletion. They inferred that PKA plays a positive role in protecting against APAP-induced hepatotoxicity. The PKA/CREB pathway has been shown to reduce oxidative stress and avoid cell apoptosis by affecting GSH synthesis or mitochondrial ROS generation (Zhu et al., 2012; Ferretti et al., 2012). In addition, other pathways might be involved in residual CREB activation by nevirapine and tienilic acid in the presence of H-89. The p38/MAPKs signaling pathway, another upstream of CREB, is known to be activated by oxidative stress during drug-induced liver injury (El-Kholy et al., 2017; Wang et al., 2018). Aguilar Mora et al. (2021) found that the EPAC signaling pathway is primarily responsible for protection against diclofenac-induced injury in hepatocytes. EPAC is also responsible for CREB phosphorylation (Ramos-Alvarez et al., 2019; Chen et al., 2014; Glas et al., 2016). Thus, it is likely that CREB, as the downstream hub of a set of upstream signaling networks, plays a key role in exerting hepatoprotective effects against the stress induced by high-risk drugs.
The present study might have raised hopes for the CREB reporter gene assay as a simplified and convenient method for evaluating IDT-associated gene expression. However, there are several issues that need to be considered and resolved. First of all, verification with a sufficient number of compounds is needed. The present experiment was only a preliminary simplified confirmation of the results of toxicogenomic data analysis. Transient reporter gene assays are labor-intensive and less stable and would not be suitable for large-scale screening. The establishment of stable cell lines is essential for increasing the throughput and large-scale validation. Furthermore, there are problems associated with the addition of the S9 mix. It would be useful as a simple method to generate reactive metabolites and assess cytotoxicity. However, metabolites are produced extracellularly, unlike in nature, hence reducing the creditability of obtained results. It should also be noted that the intrinsic bioactivity of tested drugs may affect the expression of the reporter gene. Warfarin is a well-known anticoagulant, but it is also one of coumarin derivatives with various biological activities (Garg et al., 2020). The activation of CREB by warfarin might be a cross-talk between signaling pathways associated with antioxidant and cytoprotective activities. In addition to the review by Hassanein et al. (2020), which suggested that the coumarin derivatives activate NRF2 by directly inhibiting KEAP1, we also detected activation of NRF2 by warfarin (Fig. 2). Taken together, several findings suggest that NRF2 may positive regulate CREB via PPARγ (Lee, 2017; Chiang et al., 2014; Mäkelä et al., 2016). Thus, it was possible that warfarin could activate CREB as one of its biological activities. It reminded us at the same time that the inability to distinguish whether reporter gene expression is due to cytotoxicity risk or drug intrinsic activity would be a limitation of this assay.
In conclusion, we found that CREB might be a sensitive indicator of idiosyncratic DILI. A system using a CREB vector and reporter gene assay, while requiring extensive confirmation and further refinement, might help reduce the risk of clinically relevant hepatotoxicity as a complement to other prediction tools. However, it should always be kept in mind that reporter gene assays can cause false positives due to drug intrinsic activity.
This research was mainly supported by JSPS KAKENHI Grant Numbers JP16H01861 (MH) and JP18H03531 (FY) from the Japan Society for the Promotion of Science and TaNeDS grant program (FY) from Daiichi Sankyo Co., Ltd. The authors thank Professor Hiroshi Mamitsuka of the Bioinformatics Center, Institute for Chemical Research, Kyoto University, for valuable advice on the analysis of transcriptome data, and Dr. Kazunori Fujimoto and Dr. Takuma Iguchi from Daiichi Sankyo Co., Ltd., for useful comments and discussions on toxicity assessment experiments and selection of test drugs.
Data statementThe data will be available at: https://osf.io/mvk42/
Conflict of interestFumiyoshi Yamashita reports financial support was provided by Daiichi Sankyo Co. Ltd. with respect to the research of this article.