Abstract
Genes involved in drug absorption, distribution, metabolism, and excretion (ADME) are called ADME genes. Currently, 298 genes that encode phase I and II drug metabolizing enzymes, transporters, and modifiers are designated as ADME genes by the PharmaADME Consortium. ADME genes are highly expressed in the liver and their levels can be influenced by liver diseases such as hepatocellular carcinoma (HCC). In this study, we obtained RNA-sequencing and microRNA (miRNA)-sequencing data from 371 HCC patients via The Cancer Genome Atlas liver hepatocellular carcinoma project and performed ADME gene–targeted differential gene expression analysis and expression correlation analysis. Two hundred thirty-three of the 298 ADME genes (78%) were expressed in HCC. Of these genes, almost one-quarter (58 genes) were significantly downregulated, while only 6% (15) were upregulated in HCC relative to healthy liver. Moreover, one-half (14/28) of the core ADME genes (CYP1A2, CYP2A6, CYP2B6, CYP2C8, CYP2C9, CYP2C19, CYP2E1, CYP3A4, NAT1, NAT2, UGT2B7, SLC22A1, SLCO1B1, and SLCO1B3) were downregulated. In addition, about one-half of the core ADME genes were positively correlated with each other and were also positively (AHR, ARNT, HNF4A, PXR, CAR, PPARA, and RXRA) or negatively (PPARD and PPARG) correlated with transcription factors known as ADME modifiers. Finally, we show that most miRNAs known to regulate core ADME genes are upregulated in HCC. Collectively, these data reveal 1) an extensive transcription factor–mediated ADME coexpression network in the liver that efficiently coordinates the metabolism and elimination of endogenous and exogenous compounds; and 2) a widespread deregulation of this network in HCC, most likely due to deregulation of both transcriptional and post-transcriptional (miRNA) pathways.
Introduction
The liver expresses a variety of phase I and II enzymes and transporters and thus it is the primary organ of drug and xenobiotic metabolism and clearance. Phase I enzymes mainly perform oxidative, reductive, or hydrolytic reactions that usually create or reveal a functional group (e.g., hydroxyl, carboxyl, and amino) on a xenobiotic, whereas phase II enzymes generally add a water-soluble group (e.g., glutathione, sulfate, and glucuronic acid) directly to a xenobiotic or to a metabolite of phase I metabolism. The influx and efflux transporters are responsible for the uptake and excretion of various compounds (e.g., drugs) in and out of cells. The group of genes involved in drug absorption, distribution, metabolism, and excretion (ADME) consist mainly of genes encoding phase I and II drug metabolizing enzymes, drug transporters, and modifiers; the latter are a group of factors that either modulate the expression of other ADME genes or influence the biochemistry of ADME enzymes (www.pharmaadme.org). Currently, 298 genes are classified as ADME genes by the pharmaADME Consortium, including 32 core ADME genes and 266 extended ADME genes (www.pharmaadme.org). The combined activities of the ADME genes determine the capacity of the liver to metabolize and clear drugs and hence influence the efficacy of drug treatment.
The expression and activities of ADME genes are influenced by liver diseases, such as viral infection, alcohol liver disease, primary sclerosis cholangitis, nonalcoholic fatty liver disease, and hepatocellular carcinoma (HCC) (Kurzawski et al., 2012; Hardwick et al., 2013). For example, nonalcoholic fatty liver disease significantly alters the expression of ADME genes that encode a variety of phase I and II drug metabolizing enzymes and drug transporters, including cytochrome P450 (P450) enzymes (Fisher et al., 2009), glutathione S-transferases (GSTs) (Hardwick et al., 2010), UDP-glucuronosyltransferases (UGTs) and sulfotransferases (SULTs) (Congiu et al., 2002; Hardwick et al., 2013), and ATP-binding cassette (ABC) transporters (Hardwick et al., 2011).
HCC is the fifth most common cancer and the third highest cause of cancer-related death globally (Lin et al., 2012). Common risk factors for HCC are viral infection, chronic alcohol consumption, aflatoxin exposure, and cirrhosis (Farazi and DePinho, 2006). Surgical resection, liver transplantation, and ablation remain the only potentially curative treatments for early stage HCC, whereas two multikinase inhibitors, sorafenib and regorafenib, are the only approved systemic drugs for treating advanced HCC (Llovet et al., 2008; Bruix et al., 2017). The sequential use of sorafenib (front line) and regorafenib (second line) provides a survival benefit in advanced HCC (Bruix et al., 2017).
Several studies have investigated the differential expression of ADME genes in HCC relative to adjacent noncancerous liver tissue. Many of these studies reported altered expression of a single ADME gene or gene superfamily, such as P450s (Chen et al., 2014; Yan et al., 2015a,b), ABC transporters (Borel et al., 2012), or UGTs (Yan et al., 2015a). Some high-throughput studies using microarray-based gene expression profiling (Okabe et al., 2001; Lee et al., 2004; Park et al., 2006; Jin et al., 2015a), or more recently deep RNA sequencing (RNA-seq) (Ho et al., 2015; Woo et al., 2017), have revealed altered expression of ADME genes from a number of superfamilies (e.g., P450s, UGTs, and GSTs); however, these untargeted studies did not draw specific conclusions about dysregulation of ADME pathways. Based on the limitations of these previous studies, we sought to provide the first targeted, yet comprehensive, dysregulation analysis of all known ADME genes in HCC.
The Cancer Genome Atlas (TCGA) liver hepatocellular carcinoma (LIHC) project contains RNA-seq data from 371 HCC patients including 50 paired HCC and corresponding adjacent noncancerous liver tissues (https://portal.gdc.cancer.gov/). In the present study, we performed differential gene expression analysis of the 50 paired HCC and noncancerous tissues focusing on the defined ADME gene set. Key findings were that one-half of the core ADME gene set and about one-quarter of the extended ADME set were significantly downregulated in HCC, while few genes were upregulated. Furthermore, about one-half of the core ADME genes showed significant correlation with each other and were also significantly correlated with nine transcription factors (TFs) [AHR, ARNT, HNF4A, pregnane X receptor (PXR), constitutive androstane receptor (CAR), PPARA, PPARD, PPARG, and RXRA], indicating coordinated regulation at the transcriptional level. The study also identified specific microRNAs (miRNAs) likely to be involved in ADME gene dysregulation in HCC.
Materials and Methods
Liver Hepatocellular Carcinoma Data Set of TCGA.
The Cancer Genome Atlas is a community resource project that collects and analyzes cancer specimens from a variety of human cancers (https://cancergenome.nih.gov). Data generated from TCGA projects such as RNA-seq and miRNA sequencing (miRNA-seq) for various human cancers are made freely available for browser, download, and use for commercial, scientific, and educational purposes. The results shown in the present study were mainly based on RNA-seq and miRNA-seq data generated by the TCGA Research Network of the LIHC project (https://portal.gdc.cancer.gov/projects/TCGA-LIHC). The TCGA LIHC project contains RNA-seq data of HCC tumors from 371 patients. RNA-seq data from normal liver tissues adjacent to the corresponding tumors are also available for 50 of these patients. The TCGA LIHC project also contains miRNA-seq data for these patients, including paired normal and cancerous specimens from 49 patients. As described in detail subsequently, the RNA-seq and miRNA-seq data from the paired normal and cancerous specimens were used for differential gene and miRNA expression analysis, respectively, and the RNA-seq data from the 371 tumors were used for gene expression correlation analysis.
The detailed clinical data for each of the 371 HCC patients are available from the cBioportal (http://www.cbioportal.org/study?id=lihc_tcga#clinical). This cohort includes 250 male and 121 female patients with the age at diagnosis ranging from 16 to 90 years. The ethnicities of the subjects are available for 361 patients, including Asians (158), Caucasians (184), Black or African Americans (17), and American Indian or Alaska Native (2). The American Joint Committee on Cancer proposes a staging grouping based on T, N, and M to define the three stages of tumors: the primary tumor (T), regional lymph nodes (N), and distant metastases (M) (Kamarajah et al., 2018). The progression stages of the tumors are available for 368 patients. According to this TNM staging system, all tumors were primary hepatocellular carcinomas at stage T, including 181, 94, 80, and 13 tumors at stages T1, T2, T3, and T4, respectively. Tumors at stages N and M were not included in the TCGA LIHC project. The TCGA LIHC clinical data do not provide the details of the pre/postoperative chemotherapeutic regimens of the patients.
ADME Genes Defined by the PharmaADME Consortium.
Thirty-two core and 266 extended ADME genes (a total of 298 genes) are currently defined by the PharmaADME Consortium (http://www.pharmaadme.org). The core ADME gene set represents the most important genes directly involved in drug metabolism/clearance; the extended ADME gene set represents other genes related to drug metabolism/clearance. Four extended ADME genes, either pseudogene (CYP2D7P1), antisense gene (SLC22A18AS), or functionally unknown genes (LOC728667 and LOC731931) were excluded from analysis in this study. The remaining 294 ADME genes (32 core; 262 extended) are listed in Supplemental Table 1.
Download of the TCGA LIHC RNA-Seq and miRNA-Seq Data.
For differential gene expression analysis, TCGA LIHC RNA-seq data from the 50 paired HCC and corresponding adjacent noncancerous liver tissues that have been mapped to the human genome GRCh38 reference assembly were downloaded as high-throughput sequencing (HT-seq) counts from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). The HT-seq counts contained the raw read counts of each gene generated using a HT-seq tool (http://htseq.readthedocs.io/en/release_0.9.1/) to measure gene expression level from the Genomic Data Commons mRNA quantification analysis pipeline (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/). The HT-seq counts of 609 genes that were included in the analysis in this study are given in Supplemental Table 2.
For gene expression correlation analysis, the TCGA LIHC RNA-seq data from 371 HCC tumors were downloaded in the form of HT-seq counts from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). Genes (protein coding and noncoding) with a mean of less than 10 counts were discarded; the counts of the remaining genes were normalized using the upper quantile normalization method. The normalized read counts of the 609 genes that were included in the analysis of this study are given in Supplemental Table 3.
For differential miRNA expression analysis, TCGA LIHC miRNA-seq data of the 49 paired HCC and corresponding adjacent noncancerous liver tissues were downloaded from The Cancer Genome Atlas data portal (https://portal.gdc.cancer.gov/) and normalized using the upper quantile normalization method. The normalized read counts of miRNAs that were included in the analysis of this study are given in Supplemental Table 4.
ADME Genes and Non-ADME Genes Analyzed in this Study.
Most of the 298 ADME genes belong to gene superfamilies such as ABC transporters, P450 enzymes, solute carrier (SLC) transporters, and UGT enzymes. However, many of these gene superfamilies also contain members that are not involved in drug absorption, distribution, metabolism, and excretion. These genes are termed non-ADME genes in this study. The TCGA LIHC RNA-seq data showed HT-seq counts for 743 genes, including 290 ADME genes and 453 non-ADME genes. For quality control, genes with HT-seq counts of <32 (equivalent to <5 of log 2-transformed counts) were excluded from the analysis. Based on this cutoff criterion, a total of 609 genes (Supplemental Table 2), including 233 ADME genes and 376 non-ADME genes, were included in the final analysis. Inclusion of the non-ADME genes in our analysis allowed us to determine whether ADME genes were specifically dysregulated, rather than whole gene superfamilies.
Statistical Analysis of the TCGA LIHC RNA-Seq Data and miRNA-Seq Data.
The downloaded RNA-seq (HT-seq counts) of the 50 paired HCC and corresponding adjacent normal liver tissues as described previously were subjected to differential gene expression analysis using the DESeq2 program (Love et al., 2014). Briefly, DESeq2 takes the HT-seq counts per gene and normalizes gene counts using a generalized linear model. After normalization, it estimates negative binomial distribution and performs variance shrinkage to determine the fold changes of genes between normal and cancerous tissues. DESeq2 uses a Wald test for statistically significant testing. Briefly, the shrunken estimate of log fold change is divided by its S.E., resulting in a z-statistic, which is compared with a standard normal distribution. The Wald test allows testing of individual coefficients, or contrasts of coefficients, without the need to fit a reduced model as with the likelihood ratio test. The Wald test P values from the subset of genes that pass the independent filtering step are adjusted for multiple testings using the Benjamini-Hochberg test. As detailed in Supplemental Table 5, the final report of the DESeq2 analysis included a log 2 fold change in the expression level of each gene between the tumor and adjacent nontumor liver tissues and the associated S.E. estimate for the log 2 fold change estimate and the Wald test result (P value) and the Benjamini-Hochberg test results (adjusted P values). An adjusted P value of <0.05 was considered statistically significant. Differentially expressed genes were defined by a log 2 fold change of >1 (equivalent to an upregulation of >2-fold) or <−1 (equivalent to a downregulation of <2-fold). Supplemental Table 5 gives the final results pertaining to the 609 genes that are included in analysis in this study.
The normalized gene expression levels (read counts) of the 371 HCC tumors as described previously were subjected to Spearman ranking correlation analysis 1) between core ADME genes and 2) between core ADME genes and nine ADME modifier transcription factors and NRF2 and its coregulatory KEAP1 using GraphPad Prism (version 7.03). A P value of <0.05 was considered statistically significant.
The normalized expression levels of miRNAs of the 49 paired HCC and corresponding adjacent normal noncancerous liver tissues as described previously were analyzed using a paired t test via GraphPad Prism (version 7.03 GraphPad, LA Jolla, CA). A P value of <0.05 was considered statistically significant.
Analysis of Microarray-Based Gene Expression Profiles of Hepatocellular Carcinoma.
Most of the results reported in this study were obtained from the analysis of the TCGA LIHC RNA-seq data and miRNA-seq data as described previously. To support the findings from this data set, we further analyzed microarray-based whole genome gene expression profiles from four other studies (Chen et al., 2002; Wurmbach et al., 2007; Roessler et al., 2010) via Oncomine, a publicly accessible database (Rhodes et al., 2007). As stated by Oncomine (www.oncomine.org), Roessier et al. (2010) analyzed 225 HCC and 220 adjacent nontumor liver tissues using the Affymetrix human genome 96 HT HG-U133A 2.0 microarray platform (GSE14520). Using the same U133A 2.0 microarray, Wurbach et al. (2007) studied 35 HCC tumors and 10 heathy liver tissues (GSE6764). Chen et al. (2002) analyzed HCC tumor liver tissues (71–76 specimens) and adjacent nontumor liver tissues (98–104 specimens) using custom-made cDNA microarrays containing 23,075 cDNA probes, representing 17,400 genes. Thurnherr et al. (2016) studied 100 paired HCC tumor and adjacent nontumor liver tissues using Agilent’s 60-mer oligo microarray (Agilent Technologies, Santa Clara, CA) (GSE62043).
Results
Deregulation of ADME and Non-ADME Genes in HCC.
This study examined the expression of 233 ADME genes (Supplemental Table 6) in 50 paired HCC and corresponding noncancerous liver tissues using the TCGA LIHC RNA-seq data set. Three hundred seventy-six non-ADME genes from ADME gene superfamilies that are not involved in drug metabolism and clearance were also analyzed for comparison. Overall, of the 609 genes analyzed, 109 genes (18%) were downregulated (Supplemental Table 7) and 55 genes (9%) were upregulated (Supplemental Table 8) in HCC. Among the 233 ADME genes, 58 genes (24%) were downregulated (Supplemental Table 9) and 15 genes (6%) were upregulated (Supplemental Table 10) in HCC, whereas out of the 376 non-ADME genes, 51 genes (14%) were downregulated (Supplemental Table 11) and 40 genes (11%) were upregulated (Supplemental Table 12) in HCC. Chi-square tests showed that ADME genes were significantly more likely to be downregulated than upregulated (P < 0.0001), and that ADME genes were significantly more likely to be downregulated than non-ADME genes (P < 0.001).
Deregulation of Core ADME Genes in HCC.
The PharmaADME Consortium (www.pharmaadme.org) currently defines 32 core ADME genes (Supplemental Table 13), four (GSTT1, SLC15A2, SLC22A2, and SLC22A6) of which were not expressed in the TCGA HCC specimens. Of the remaining 28 core ADME genes, 14 (50%) were significantly downregulated in HCC tissues compared with paired adjacent noncancerous liver tissues, including eight phase I enzymes (CYP1A2, CYP2A6, CYP2B6, CYP2C8, CYP2C9, CYP2C19, CYP2E1, and CYP3A4), three phase II enzymes (NAT1, NAT2, and UGT2B7), and three transporters (SLC22A1, SLCO1B1, and SLCO1B3) (Fig. 1; Table 1). None of the core ADME genes were upregulated. Chi-square tests showed that core ADME genes were significantly more likely to be downregulated relative to extended ADME genes (P < 0.01) or non-ADME genes (P < 0.0001). Our analyses of four other microarray gene expression profiling studies of HCC (Chen et al., 2002; Wurmbach et al., 2007; Roessler et al., 2010; Thurnherr et al., 2016) showed that each of the 14 downregulated core ADME genes was reported to be downregulated in HCC in at least one of these other studies (Table 2).
Deregulation of ADME Genes Coding for Phase I Drug Metabolism Enzymes in HCC.
One hundred and thirty of the 298 ADME genes code for phase I drug metabolizing enzymes (PharmaADME). Of these genes, 106 (82%) were expressed in the TCGA HCC specimens, including 34 downregulated genes (32%) (Table 1) and five upregulated genes (5%): ALDH3A1, CBR3, CYP17A1, NOS2A, and CYP7A1. Among the downregulated genes, 27 encode members of three enzyme superfamilies, cytochromes P450 enzymes, alcohol dehydrogenases (ADHs), and aldehyde dehydrogenases (ALDHs). The remaining seven genes were AOX1, DHRS1, EPHX2, PLGLB1, PON1, PON3, and XDH.
There are 57 putatively functional CYP genes in the human genome that are grouped into 18 families and 44 subfamilies based on sequence similarity (Nelson et al., 2004). Forty-seven CYP genes are classified as ADME genes (PhamaADME). Of these ADME CYP genes, 37 were expressed in the TCGA HCC specimens: two of these were upregulated (CYP7A1 and CYP17A1) and 17 were downregulated (Fig. 2A; Table 1). The downregulated CYP genes included the eight aforementioned downregulated core ADME CYP genes (Fig. 1) and nine other CYP genes (CYP2C18, CYP3A7, CYP3A43, CYP4A11, CYP4F2, CYP4F12, CYP8B1, CYP26A1, and CYP39A1).
The human genome has seven functional ADH genes (1A, 1B, 1C, 4, 5, 6, and 7) that are all classified as ADME genes based on their roles in drug clearance (Edenberg, 2007). Of these ADH genes, six (1A, 1B, 1C, 4, 5, and 6) were expressed in the TCGA HCC specimens, five of which (1A, 1B, 1C, 4, and 6) were significantly downregulated in HCC, consistent with the findings of recent reports (Wei et al., 2012; Shang et al., 2017).
There are 19 ALDH genes in humans (Marchitti et al., 2008), 15 of which are classified as ADME genes (www.pharmaadme.org). Of the ADME ALDH genes, 14 were expressed in the TCGA HCC specimens, including five genes that were downregulated (ALDH1A3, 1B1, 2, 6A1, and 8A1) and one that was upregulated (ALDH3A1) in HCC (Fig. 2B; Table 1). The downregulation of ALDH1B1 (Yang et al., 2017) and ALDH2 (Jin et al., 2015b) in HCC was consistent with recent reports, whereas the downregulation of ALDH1A3, ALDH6A1, and ALDH8A1 and the upregulation of ALDH3A1 in HCC were novel findings.
Deregulation of ADME Genes Coding for Phase II Drug Metabolism Enzymes in HCC.
Sixty-eight of the 298 ADME genes encode for phase II drug metabolizing enzymes (PharmaADME). Of these ADME genes, 52 (76%) were expressed in the TCGA HCC specimens, including 11 genes (23%) that were downregulated and three genes (6%) that were upregulated in HCC. Seven downregulated genes encode for members of three enzyme superfamilies (GSTs, SULTs, and UGTs), as described in detail subsequently, and the other four genes were carbohydrate (N-acetylglucosamine 6-O) sulfotransferase 4 (CHST4), nicotinamide N-methyltransferase (NNMT), and two arylamine N-acetyltransferases (NAT1 and NAT2). The upregulated genes were CHST10, SULT1C2, and UGT2B11.
The human genome contains 25 GST genes (www.genenames.org), 21 of which are classified as ADME genes (www.pharmaadme.org). Of these ADME GST genes, 17 genes from three GST subfamilies (cytosolic, mitochondrial, and microsomal) were expressed in the TCGA HCC specimens, including two GST genes (GSTM5 and GSTZ1) that were significantly downregulated in HCC (Fig. 3A). GSTZ1 downregulation in HCC was recently reported (Jahn et al., 2016). The downregulation of GSTA1 (Li et al., 2008) and GSTP1 (Zhong et al., 2002; Zhang et al., 2005; Shen et al., 2012) and the upregulation of GST1A4 (Liu et al., 2017) in HCC were also previously reported. There was no overall significant deregulation of these three GST genes in HCC, but we observed reduced GSTA1 and GSTP1 mRNA levels and increased GST1A4 mRNA levels in HCC tumor tissues compared with corresponding adjacent normal liver tissues in nearly 50% of the patients (data not shown).
There are 13 SULT genes in the human genome, 10 of which are classified as ADME genes (Blanchard et al., 2004). Of these SULT genes, six (1A1, 1A2, 1B1, 1C2, 1E1, and 2A1) were expressed in the TCGA HCC specimens, including three genes that were downregulated (1A2, 1E1, and 2A1) and one gene that was upregulated (1C2) in HCC (Fig. 3B). In agreement with these findings, the downregulation of SULT1E1 and SULT2A1 in HCC was recently reported (Xie et al., 2017).
There are 21 functional UGT genes in the human genome that are subdivided into four families: UGT1, UGT2, UGT3, and UGT8 (Hu et al., 2014a). Eighteen UGT genes are classified as ADME genes, including four core ADME genes (1A1, 2B15, 2B17, and 2B7) and 14 extended ADME genes (1A3–10, 2A1, 2B4, 2B10, 2B11, 2B28, and UGT8) (www.pharmaadme.org). Of the ADME UGT genes, 12 (1A1, 1A3, 1A4, 1A6, 1A9, 2A1, 2B4, 2B7, 2B10, 2B11, 2B15, and 2B17) were expressed in the TCGA HCC specimens (Fig. 3C; Table 1). We show here for the first time, UGT2B10 downregulation and UGT2B11 upregulation in HCC (Fig. 3C). We also show UGT2B7 downregulation in HCC, a finding that was recently reported by others (Lu et al., 2015; Yan et al., 2015a). Previous studies have also reported downregulation of three other UGTs (1A1, 1A4, and 1A9) in HCC (Lu et al., 2015; Yan et al., 2015a). There was no overall significant deregulation of these three UGT genes in the TCGA HCC specimens (Fig. 3C) but their mRNA levels were decreased in HCC compared with paired adjacent noncancerous liver tissues in nearly 50% of the patients (Supplemental Fig. 1). As expected, the mRNA levels of the three extrahepatic UGTs (1A7, 1A8, and 1A10) were extremely low in the TCGA HCC specimens (data not shown).
Deregulation of ADME Genes Coding for Drug Transporters in HCC.
Seventy-seven of the 298 ADME genes encode for drug transporters (PharmaADME) that belong to two transporter superfamilies: ABC and SLC transporters. Of these ADME genes, 54 (70%) were expressed in the TCGA HCC specimens, including nine genes (17%) that were downregulated and seven genes (13%) that were upregulated in HCC. The details of the deregulation of ABC and SLC genes are described subsequently.
The human genome contains 48 functional ABC genes (Dean et al., 2001; Wlcek and Stieger, 2014), among which 26 genes are defined as ADME genes (PharmaADME). Of these ADME ABC genes, three (ABCA4, ABCC4, and ABCC10) were upregulated and one (ABCC9) was downregulated in HCC tumor (Fig. 4A). The expression of three core ADME genes (ABCB1, ABCC2, and ABCG2) showed high interindividual variabilities in HCC cancerous and noncancerous liver tissues (Chou et al., 1997; Richart et al., 2002; Zollner et al., 2005), and the degregulation of these ABC genes in HCC has also been reported (Moustafa et al., 2004; Zollner et al., 2005; Borel et al., 2012). Indeed, we observed similar high interindividual expression variabilities for these three genes in the TCGA HCC specimens with upregulation or downregulation in HCC relative to corresponding adjacent noncancerous liver tissue in nearly 30% of the patients (data not shown).
The human genome contains 395 SLC genes (Hediger et al., 2013), among which 51 genes are defined as ADME genes (PharmaADME). Of these ADME SLC genes, 34 were expressed in the TCGA HCC specimens, including four genes that were upregulated (SLC22A4, SLC22A11, SLC22A15, and SLCO2A1) and eight genes that were downregulated (SLC10A1, SLC22A1, SLC22A10, SLC7A8, SLCO1B1, SLCO1B3, SLCO2B1, and SLCO4C1) in HCC (Fig. 4B; Table 1). The downregulation of three SLC genes (SLC22A1, SLCO1B1, and SLCO1B3) in TCGA HCC was consistent with recent reports (Vavricka et al., 2004; Vander Borght et al., 2005; Tsuboyama et al., 2010; Heise et al., 2012). The downregulation of SLC22A3 in HCC was also reported (Heise et al., 2012). This gene was downregulated in nearly 50% of the TCGA HCC specimens. Northern blot analysis showed that SLCO4C1 is only expressed in kidney, and thus is considered a kidney-specific SLC (Mikkaichi et al., 2004). We showed here that this gene was abundantly expressed in noncancerous liver tissues but was significantly downregulated in HCC (Fig. 4B).
Deregulation of ADME Genes Coding for Modifiers in HCC.
Twenty-four of the 298 ADME genes encode for modifiers that either modulate the expression of ADME genes or affect the biochemistry of ADME enzymes (PharmaADME). Of these ADME modifiers, 21 were expressed in TCGA HCC specimens, including nine transcription factors (AHR, ARNT, HNF4A, PXR, CAR, PPARA, PPARD, PPARG, and RXRA) (Fig. 5; Table 1). Among these expressed modifier genes, four (catalase, cytidine deaminase, MAT1A, and PXR) were downregulated in HCC and one (KCNJ11) was upregulated in HCC. The upregulation of KCNJ11 (Zhang et al., 2018) and the downregulation of catalase (Cho et al., 2014), PXR (Chen et al., 2014), or cytidine deaminase (Nwosu et al., 2017) in HCC have been recently reported. KCNJ11 is an HCC oncogene and its overexpression in HCC promotes cell proliferation and tumor progression (Zhang et al., 2018). MAT1A upregulation in HCC was also recently reported (Nwosu et al., 2017); however, the expression of this gene was significantly downregulated in the TCGA HCC specimens. This discrepancy awaits further investigation.
In addition to the nine transcription factors defined as ADME modifiers by the PharmaADME Consortium, many other transcription factors are involved in transcription regulation of ADME genes (Hu et al., 2014a; Zanger and Schwab, 2013). For example, a large number of ADME genes are regulated by the oxidative stress-responsive transcription factor Nrf2 (Hayes and Dinkova-Kostova, 2014). Mutations of NRF2 and its negative regulator KEAP1 in HCC have been reported (Guichard et al., 2012). There was no significant overall deregulation of NRF2 and KEAP1 in TCGA HCC specimens (Fig. 5); however, we show here that NRF2 had reduced mRNA levels in HCC in 42 of the 50 patients (Supplemental Fig. 2).
Correlation Analyses Between Core ADME Genes and Nine ADME Modifiers Coding for Transcription Factor Genes in HCC.
As mentioned previously, nine ADME modifiers (AHR, ARNT, HNF4A, PXR, CAR, PPARA, PPARD, PPARG, and RXRA) are transcription factors that are believed to be able to alter the expression of other ADME genes. Indeed, three of these modifiers (PXR, CAR, and HNF4A) are involved in transcriptional regulation of some ADME genes in liver (Xu et al., 2005; Congiu et al., 2009; Zhong et al., 2017). Analysis of TCGA RNA-seq data from 371 HCC tumor tissues revealed that most of the core ADME genes were significantly positively correlated to five TF modifiers (HNF4A, PXR, CAR, PPARA, and RXRA), and that about one-half of the core ADME genes were significantly positively correlated to two other TF modifiers (AHR and ARNT) (Supplemental Table 14). In particular, PXR was significantly positively correlated with 25 core ADME genes (Fig. 6). By contrast, PPARD was only positively correlated with one core ADME gene (DPYD) but was negatively correlated to 17 core ADME genes (Fig. 7). Similarly, PPARG was positively correlated with only four core ADME genes (ABCC2, CYP2E1, SLCO1B3, and UGT1A1) but negatively correlated with 14 core ADME genes (Supplemental Table 14). Of note, GSTP1 was the only core ADME gene that was significantly negatively correlated with most of the modifier TFs (seven out of nine) (Fig. 8). We further examined whether these TFs were differentially expressed in HCC in individual patients relative to their paired adjacent noncancerous liver tissues. As expected from an overall downregulation of PXR in HCC (Fig. 5), it was downregulated in HCC in almost all patients (Supplemental Fig. 2). There was no significant overall deregulation of AHR, CAR, PPARD, and PPARG (Fig. 5), but we show here that AHR and CAR were downregulated in HCC in nearly 50% of the patients, and that there was upregulation of PPARD and PPARG in HCC in nearly 50% of the patients (Supplemental Fig. 2). Taken together, about one-half of core ADME genes are positively correlated to seven ADME TFs (AHR, ARNT, HNF4A, PXR, CAR, PPARA, and RXRA) and negatively correlated to two ADME TFs (PPARD and PPARG) in HCC. Furthermore, NRF2 was also positively correlated with more than 20 core ADME genes (Supplemental Table 14). Given that six of these 11 transcription factors and their coregulators are deregulated in HCC, it is likely that they contribute to the widespread downregulation of core ADME genes in HCC.
The Expression Levels of Core ADME Genes Are Correlated with Each Other in HCC.
The correlation of the core ADME genes with ADME TFs described previously suggested that core ADME genes are coregulated. To examine this idea, we performed expression correlation analysis of all 28 core ADME genes. Twenty-two core ADME genes (ABCB1, ABCC2, ABCG2, CYP1A1, CYP1A2, CYP2A6, CYP2B6, CYP2C8, CYP2C9, CYP2D6, CYP3A4, CYP3A5, NAT1, NAT2, SLC22A1, SLCO1B1, SULT1A1, TPMT, UGT1A1, UGT2B7, UGT2B15, and UGT2B17) were significantly positively correlated with at least 20 other core ADME genes (Supplemental Table 15); five core ADME genes (CYP2C19, CYP2E1, DPYD, GSTM1, and SLCO1B3) were significantly positively correlated with about one-half of the core ADME genes (Supplemental Table 15). For instance, UGT2B7 was significantly positively correlated with 25 core ADME genes (Fig. 9). GSTP1 was the only core ADME gene that was significantly negatively correlated to more than one-half of the core ADME genes (Supplemental Table 15), consistent with the aforementioned negative correlation of GSTP1 with seven ADME TFs. Taken together, the widespread expression correlation of core ADME genes with each other and with transcription factors as summarized in Supplemental Table 16 support the hypothesis that ADME genes are coordinately regulated by transcription factors in the liver.
Deregulation of miRNAs that Are Known to Regulate Core ADME Genes in HCC.
miRNAs generally repress gene expression through mRNA degradation or translational repression. Therefore, the upregulation of miRNAs that regulates ADME genes might contribute to their downregulation in HCC. To test this hypothesis, we examined the expression levels of 12 miRNAs (miR-132, miR-128, miR-130b, miR-29a-3p, miR-27a, miR-107, miR-25, miR-629, miR-103-5p, miR-378a, miR-126, and miR-142) (Supplemental Table 4) that are reported to target the downregulated core ADME genes (CYP1A2, CYP2C9, CYP2C9, CYP2C19, CYP3A4, CYP2C8, CYP2B6, SLCO1B3, CYP2C8, CYP2E1, CYP2A6, and UGT2B7, respectively) in HCC and paired adjacent noncancerous liver tissues (Mohri et al., 2010; Zhang et al., 2012; Nakano et al., 2015; Rieger et al., 2015; Shi et al., 2015; Yu et al., 2015a,b; Higuchi et al., 2016; Jin et al., 2016; Chen et al., 2017; Papageorgiou and Court, 2017). Strikingly, most (eight out of 12) of these miRNAs were upregulated in HCC tumors (Fig. 10), suggesting that they may play an important role in core ADME gene dysregulation in HCC.
Deregulation of Aldo-Keto Reductases in HCC.
Aldo-keto reductases (AKRs) are involved in the phase I metabolism of numerous endogenous and xenobiotic compounds including steroids, carbohydrates, prostaglandins, aldehydes, ketones, and many therapeutic drugs such as the cytotoxic anthracyclines (e.g., doxorubicin, daunorubicin, and idarubicin) (Novotna et al., 2008; Skarka et al., 2011; Matsunaga et al., 2012; Penning, 2015). There are 15 AKR genes in the human genome (Penning, 2015) but none of these genes are classified as ADME genes by the PharmaADME consortium (pharmadame.org). We show here that 11 human AKR genes (e.g., AKR1A1, AKR1B1, AKR1B10, AKR1C1, AKR1C2, AKR1C3, AKR1C4, AKR1D1, AKR1E2, AKR7A2, and AKR7A3) were expressed in TCGA HCC specimens (Fig. 11). Of these genes, AKR1B10 and AKR1C3 were significantly upregulated, whereas AKR1D1 and AKR7A3 were significantly downregulated in HCC. AKR1B10 upregulation is consistent with the overexpression of AKR1B10 enzyme in HCC as reported recently in an immunohistochemistry study (Matkowskyj et al., 2014). AKR7A3 downregulation in HCC (Chow et al., 2016) and AKR1C3 upregulation in prostate and breast cancer (Lin et al., 2004; Fung et al., 2006) has also been reported.
Discussion
ADME genes are involved in drug absorption, distribution, metabolism, and excretion and are abundantly expressed in the liver, the major organ of drug metabolism and clearance in humans. Although several studies have investigated the deregulation of subsets of ADME genes in HCC (Okabe et al., 2001; Lee et al., 2004; Park et al., 2006; Borel et al., 2012; Chen et al., 2014; Ho et al., 2015; Jin et al., 2015a; Yan et al., 2015a,b; Woo et al., 2017), there has been no systematic, targeted analysis of the entire ADME gene set. In this study, our comprehensive analysis of paired cancerous and noncancerous specimens from the TCGA HCC data set identified for the first time the complete list of ADME genes that are deregulated in HCC. Overall, around one-quarter of all ADME genes were downregulated while only 6% were upregulated in HCC. The deregulated genes included members of all four groups of ADME genes, namely, phase I (CYPs, ADHs, and ALDHs) and phase II (GSTs, SULTs, and UGTs) drug metabolizing enzymes, drug transporters (ABCs and SLCs), and modifiers. Moreover, one-half of the core ADME genes were significantly downregulated in HCC. Among the ADME gene superfamilies, CYPs were the most significantly downregulated. Our analyses of four other microarray gene expression profiling reports of HCC observed similar downregulation of these core ADME genes in HCC (Chen et al., 2002; Wurmbach et al., 2007; Roessler et al., 2010; Thurnherr et al., 2016). The deregulation of many other ADME genes reported by us in the present study is consistent with the findings of previous studies as already discussed comprehensively. However, we report in the present study for the first time on the deregulation of a variety of ADME genes in HCC, such as AKR1C3, AKR1D1, ALDH1A3, ALDH3A1, ALDH6A1, ALDH8A1, UGT2B10, and UGT2B11. Most ADME genes belong to superfamilies that also contain non-ADME members (for example, genes that mediate biosynthetic processes or transport endobiotics). To better understand the specificity of ADME deregulation, we compared the ADME gene set to non-ADME genes from the same superfamilies. Our results showed that the ADME genes were far more likely to be downregulated than non-ADME genes in HCC. We further performed unbiased pathway analysis (Reactome) (https://reactome.org) using the entire ADME and non-ADME gene set (609 genes) and found that the downregulated gene set was significantly more enriched for drug metabolism pathways than the unchanged or upregulated gene sets (data not shown). Together with the data on ADME TFs and miRNAs discussed subsequently, this supports the idea that ADME-specific regulatory pathways are deregulated in HCC.
Decreased expression of genes encoding hepatocyte-specific gene products (e.g., albumin) and detoxification enzymes (CYPs) and associated loss of liver-specific functions are recognized phenomena in HCC and are in part due to de-differentiation of cancer cells (Okabe et al., 2001; Park et al., 2006; Andrisani et al., 2011; Obaidat et al., 2012). The widespread downregulation of ADME genes in HCC shown herein is likely part of this de-differentiation. The molecular mechanisms underlying de-differentiation are not completely understood but could be related to epigenetic (chromatin) modifications or deregulation of transcription factors and miRNAs that regulate liver-specific genes. The deregulation of TFs and miRNAs in HCC, as we showed herein in this study, supports this assumption. Briefly, we found that most of the miRNAs known to regulate the downregulated core ADME genes were upregulated in HCC. Furthermore, transcription factors (PXR, AHR, and CAR) that are positively correlated with core ADME genes were downregulated in HCC, whereas the transcription factors (PPARD and PPARG) that are negatively correlated with core ADME genes were upregulated. The recent reports of DNA hypermethylation at the promoters of three downregulated core ADME genes (CYP1A2, CYP2C9, and SLC22A1) in HCC links epigenetic modifications to their deregulation in this cancer (Schaeffeler et al., 2011; Shen et al., 2012).
There are no curative therapeutic options for advanced HCC (Llovet et al., 2008). Most of the clinical trials evaluating anticancer agents in chemotherapy, hormone therapies, and immunotherapy have failed to show convincing survival benefits for patients with advanced HCC (Llovet et al., 2008; Llovet and Hernandez-Gea, 2014). The failure of these trials might be related to the expression of nearly 80% of ADME genes (233 out of 298) in HCC as we found in the present study. Two multikinase inhibitors, sorafenib and regorafenib, are currently the only two drugs approved by the Food and Drug Administration for treating advanced HCC (Llovet et al., 2008; Bruix et al., 2017). The deregulation of the ADME genes involved in the metabolism and clearance of these two drugs may help to explain their relative efficacy. The key ADME genes involved in the metabolism and clearance of sorafenib and regorafenib in liver include the enzymes CYP3A4 (Ghassabian et al., 2012) and UGT1A9 (Peer et al., 2012; Miners et al., 2017), influx transporters (SLCO1B1, SLCO1B3, and SLC22A1) (Swift et al., 2013; Zimmerman et al., 2013; Ohya et al., 2015), and efflux transporters (ABCB1, ABCC2, and ABCG2) (van Erp et al., 2009; Ohya et al., 2015). While the efflux transporters were not significantly altered in HCC, the expressions of CYP3A4, UGT1A9, SLCO1B1, SLCO1B3, and SLC22A1 were all significantly downregulated in HCC, consistent with previous observations (Vander Borght et al., 2005; Tsuboyama et al., 2010; Heise et al., 2012; Ye et al., 2014). The downregulation of CYP3A4 and UGT1A9 could reduce the metabolism of sorafenib and regorafenib, and hence enhance therapeutic efficacy; however, the downregulation of influx transporters (SLCO1B1, SLCO1B3, and SLC22A1) may decrease the cellular uptake and reduce efficacy. Ultimately, the combined effects of deregulated uptake and metabolism will determine the therapeutic efficacy of sorafenib and regorafenib. The opposing influence of ADME gene deregulation on the therapeutic efficacy of this effective HCC therapeutic regimen warrants further investigation.
Our observation of deregulated UGT expression may also provide some insights into cytotoxic drug efficacy in HCC. The chemotherapy regimen epirubicin, 5-fluoruracil, and cisplatin has been found to be poorly effective as a first line therapy for HCC; however, it did provide survival benefits when used as a second line therapy for advanced sorafenib-resistant HCC (Lee et al., 2014). UGT2B7 is the only UGT known to mediate epirubicin conjugation (Innocenti et al., 2001), and this metabolism is suggested to be responsible for its lower systemic toxicity relative to doxorubicin (Hu et al., 2014b). We observed downregulation of UGT2B7 in HCC as well as upregulation of UGT2B11. In recent studies performed in vitro, we found that UGT2B11 inhibits the ability of UGT2B7 to glucuronidate epirubicin (Meech et al, manuscript in preparation). Thus, the combined downregulation of UGT2B7 and upregulation of UGT2B11 could reduce intratumoural epirubicin conjugation and clearance, and this may contribute to the modest efficacy of this regimen in advanced HCC.
CBR1 and multiple AKR subfamily 1 enzymes (AKR1A1, AKR1B1, AKR1B10, AKR1C3, and AKR1C4) are involved in metabolism of doxorubicin, a common cytotoxic drug for treating HCC (Kassner et al., 2008; Novotna et al., 2008). We show in this study that all of these enzymes were highly expressed in HCC and two of them (AKR1B10 and AKR1C3) were significantly upregulated in TCGA HCC specimens. The abundant expression of these enzymes and their deregulation in HCC may contribute to the development of cancer resistance to doxorubicin-containing HCC chemotherapeutic regimens (Le Grazie et al., 2017).
It has been proposed that the process of absorption, distribution, metabolism, and excretion of xenobiotics (e.g., drugs and dietary constituents) is coordinated in the liver mainly through the coregulation of ADME genes by transcription factors (Congiu et al., 2009) (Fig. 12). The direct evidence to support this hypothesis came from the discoveries that a large number of ADME genes (e.g., CYP, UGT, ABC, and SLC transporters) are induced (rarely repressed) by ligand-activated transcription factors (e.g., PXR, FXR, CAR, PPAR, AhR, and Nrf2) upon exposure to xenobiotics (e.g., drugs) (Scotto, 2003; Chen et al., 2013; Zanger and Schwab, 2013; Hu et al., 2014a). Further indirect evidence to support this hypothesis came from the studies that have shown positive correlation of the mRNA levels between ADME genes and transcription factors such as PXR, CAR, and HNF4A in normal (Slatter et al., 2006; Wortham et al., 2007) and non-HCC diseased (Congiu et al., 2002, 2009; Kurzawski et al., 2012) liver tissues. However, these correlation studies mainly relied on quantitative real-time polymerase chain reaction to measure gene expression levels and thus only analyzed a limited number of ADME genes. In this study, a comprehensive expression correlation analysis of the TCGA HCC RNA-seq data showed for the first time that most core ADME genes were significantly positively correlated with each other. We further showed for the first time that the expression levels of about one-half of the core ADME genes were significantly positively (AHR, ARNT, HNF4A, PXR, CAR, PPARA, and RXRA) or negatively (PPARD and PPARG) correlated with those of the nine transcription factors that are classified as ADME modifiers. Collectively, these observations provide compelling evidence supporting a coordinated hepatic regulation of the broader ADME transcriptome. Moreover, this study suggests that the coordinate regulation model of ADME genes may also extend to post-transcriptional regulation by miRNAs (Fig. 12); however, further work will be required to define the ADME gene/miRNA regulatory network.
In conclusion, this study reveals a widespread deregulation of ADME genes in HCC and widespread expression correlations of ADME genes to each other and with transcription factors. The deregulation of ADME genes represents part of the de-differentiation of liver function in HCC and is likely due to deregulation of epigenetics, transcription factors, and miRNAs. The coregulation of ADME genes mediated by liver-enriched (e.g., HNF4A, PPARA, PPARD, and PPARG) and xenobiotic-activated (e.g., AHR, PXR, CAR, and NRF2) transcription factors ensures that the process of absorption, distribution, metabolism, and excretion is coordinated efficiently upon exposure to endogenous and exogenous compounds such as dietary constituents and therapeutic drugs (Fig. 12).
Authorship Contributions
Participated in research design: Hu, Marri, McKinnon, Mackenzie, Meech.
Conducted experiments: Hu, Marri.
Performed data analysis: Hu, Marri, Meech.
Wrote or contributed to the writing of the manuscript: Hu, Meech, Marri, Mackenzie, McKinnon.
Footnotes
- Received November 12, 2018.
- Accepted December 19, 2018.
This study was supported by the National Health and Medical Research Council of Australia [Grant ID 1143175 (to R.M., R.A.M., P.I.M., and D.G.H.), Grants ID1020931 (to P.I.M.) and ID1085410 (to P.I.M., R.A.M., and R.M.)]. The project was also supported by funding from the Flinders Medical Centre Foundation.
R.A.M. is a Cancer Council/SA Health Beat Cancer Professorial Chair. During the preparation period, P.I.M. was a National Health and Medical Research Council Senior Principal Research Fellow and R.M. was an Australian Research Council Future Fellow.
↵This article has supplemental material available at jpet.aspetjournals.org.
Abbreviations
- ABC
- ATP-binding cassette
- ADH
- alcohol dehydrogenase
- ADME
- absorption, distribution, metabolism, and excretion
- AKR
- aldo-keto reductases
- ALDH
- aldehyde dehydrogenase
- CAR
- constitutive androstane receptor
- GST
- glutathione S-transferase
- HCC
- hepatocellular carcinoma
- HT-seq
- high-throughput sequencing
- LIHC
- liver hepatocellular carcinoma
- miRNA
- microRNA
- miRNA-seq
- microRNA sequencing
- P450
- cytochrome P450
- PXR
- pregnane X receptor
- RNA-seq
- RNA sequencing
- SLC
- solute carrier
- SULT
- sulfotransferase
- TCGA
- The Cancer Genome Atlas
- TF
- transcription factor
- UGT
- UDP-glucuronosyltransferase
- Copyright © 2019 by The American Society for Pharmacology and Experimental Therapeutics