2020 Volume 45 Issue 8 Pages 449-473
Although peroxisome proliferator-activated receptor α (PPARα) agonists are obviously hepatocarcinogenic in rodents, they have been widely used for dyslipidemia and proven to be safe for clinical use without respect to the species difference. It is established that PPARα acts as a part of the transcription factor complex, but its precise mechanism is still unknown. Using the data of Toxicogenomics Database, reliable genes responsive to PPARα agonists, clofibrate, fenofibrate and WY-14,643, in rat liver, were extracted from both in vivo and in vitro data, and sorted by their fold increase. It was found that there were many genes responding to fibrates exclusively in vivo. Most of the in vivo specific genes appear to be unrelated to lipid metabolism and are not upregulated in the kidney. Fifty-seven genes directly related to cell proliferation were extracted from in vivo data, but they were not induced in vitro at all. Analysis of PPAR-responsive elements could not explain the observed difference in induction. To evaluate possible interaction between neighboring genes in gene expression, the correlation of the fold changes of neighboring genes for 22 drugs with various PPARα agonistic potencies were calculated for the genes showing more than 2.5 fold induction by 3 fibrates in vivo, and their genomic location was compared with that of the human orthologue. In the present study, many candidates of genes other than lipid metabolism were selected, and these could be good starting points to elucidate the mechanism of PPARα agonist-induced rodent-specific toxicity.
PPARα agonists, such as fibric acids, are frequently used for dyslipidemia, especially with a high level of triglycerides (Desvergne and Wahli, 1999; Han et al., 2017). However, it is well known that repeated administration of fibrates to rodents induces hepatic hyperplasia within a few days and hepatic tumor in months (Desvergne and Wahli, 1999). The hyperplasia and the tumor have been clearly attributed to the activation of PPARα in studies using knockout mice (Lee et al., 1995). Based on citations for clinical usage of this type of drug, it could be said that there is no concern about the carcinogenicity of PPARα agonists in humans, but the essence of the obvious species difference is still unclear (Holden and Tugwood, 1999; Tateno et al., 2015; Han et al., 2017). Even though there are drugs that are clearly carcinogenic in rodents but used in the clinical field, the evaluation of carcinogenicity of newly developing drugs in rodents is still required by regulatory agencies. Furthermore, it is thought that some potentially useful chemicals are abandoned due to toxicity testing based on a known carcinogenicity in rodents even if it is due to an unknown rodent-specific mechanism. Therefore, the elucidation of the species difference in the hepatocarcinogenesis of PPARα agonists should be important. As activation of PPARα directly leads to gene expression changes, we expect that comprehensive analysis of transcription is the most effective strategy. However, understanding of this problem appears to be stagnated because knowledge of the mechanism of transcription control is presently immature, and it has been revealed that the RNA world is much wider than previously expected, as it consists of coding and non-coding RNA (Cech and Steitz, 2014).
In 2002, fifteen Japanese pharmaceutical companies and the National Institute of Health Science, the National Institute of Biomedical Innovation, started a project that continued for 10 years. The aim of the project, Toxicogenomics Project Japan (TGPJ), was to create a transcriptome database of the rat liver and kidney using more than 150 chemicals consisting of mostly medicinal compounds, and subsequently to make non-clinical toxicity testing more efficient and rapid (Urushidani, 2008). At present, the transcriptome database is open to the public together with various toxicity data as Open TG-GATEs (Genomics Assisted Toxicity Evaluation System developed by TGPJ). The database contains 4 typical fibrates, and the analysis done in the middle of the project (at the stage that 32 chemicals were stored in the database) revealed that PPARα agonists were clearly separated by hierarchical clustering and subsequently, the database proved to be quite useful for the analysis of these types of drugs (Tamura et al., 2006). Although supervised or unsupervised machine learning on transcriptome data is a powerful tool to predict a certain toxicological phenotype based on classification of the data, there exists a limit especially in elucidation of biological and toxicological mechanisms. For this purpose, confirmatory experiments should be performed focusing on a promising target. In the present study, we extracted PPARα agonist-responsive genes that are worth studying in vivo and in vitro using TG-GATEs, and made an analysis of their genomic location to provide additional useful information for further study.
Transcriptome data are from the database, TG-GATEs (http://toxico.nibiohn.go.jp/). The precise protocols are also available from the web page. In general, for the in vivo protocol, 7-week-old male SD rats were treated with three doses of each test drug (low, middle, high, and vehicle: 5 rats each) and sacrificed 3, 6, 9, and 24 hr after single oral dosing as well as 24 hr after repeated dosing for 3, 7, 14, and 28 days. Toxicology data were obtained from 5 animals while gene expression analysis by Affymetrix Gene Chip (Rat Expression Array 230 2.0) was performed for 3 out of 5. As for the in vitro protocol, isolated hepatocytes from 7-week-old SD rats according to the standard method (Tamura et al., 2006) were cultivated in a collagen-coated six-well plate for 24 hr, and 2, 8, and 24 hr after the addition of the test drug (low, middle, high, and vehicle: 2 wells each) and the cells were harvested and analyzed for gene expression. The doses, the vehicle used, and the abbreviations for the drugs analyzed in the present study are summarized in Table 1.
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/45_449-t1.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
Extraction of PPARα agonist-responsive genes in the liver
The digital image files were processed by Affymetrix Microarray Suite version 5.0 and the intensities were normalized for each chip by setting the mean intensity to 500 (per chip normalization). When the absolute expression levels were compared between in vivo rat liver and in vitro rat hepatocytes, the values were normalized by that of peptidylprolyl isomerase A (Ppia, Affymetrix probe ID: 1398850_at) as an internal standard.
As for the data in vivo (N=3 each), the ratio to control of their means and effect size (difference between the means of control and treated group divided by the pooled standard deviation) were calculated for each probe set. When the “present calls” indicated by the system were less than 3 out of 6, the ratio was considered to be not confident, and obligatorily set to 1 (no change) in order to avoid abnormal values.
In vitro data consist of two measures for one point, where they were 2 wells of cultured cells from an animal, meaning one duplicate value, and thus statistical analysis was impossible. Therefore, ratio to control of their means was calculated and the evaluation factor Ef was defined as follows: Ef = (A–B)/C, where A is the smaller one of the pair that has the larger mean value than the other pair; B is the larger one of the pair that has the smaller mean value than the other pair; C is the difference of the means of the two pairs. By this way, the closer the duplicate and the larger the difference of the means, the larger Ef value.
After several preliminary analyses (see results section), it was decided that the first analysis was focused on the data of 24 hr after single administration of high dose of three fibrates, clofibrate (CFB), fenofibrate (FFB), and WY-14,643 (WY). Filtering of in vivo data was done such that ratio to control > 1.5 and effect size > 0.66. As for in vitro data, ratio to control > 1.3 and Ef > 0.3 were selected. Finally, the intersection, common probe sets among these three fibrates were selected both in vivo and in vitro.
Even by the latest information supplied by Affymetrix (http://www.affymetrix.com/analysis/index.affx), a considerable number of extracted probe sets have no annotations. In order to improve the information about the potential gene candidates, a BLAST search for the target sequence was performed for the probe sets without annotation. The probe sets shown below are our personal annotation. These are shown in parentheses in the tables.
1. 1384474_at: cholinergic receptor, nicotinic, alpha 2 (Chrna2)
The target sequence is 83% identical to mouse Chrna2 (9e-80) and another probe set for Chrna2 (1387574_at) shows almost identical results.
2. 1391500_at: aquaporin 7 (Aqp7)
The target sequence is 82% identical to mouse Aqp7 (4e-55), and another probe set for Aqp7 (1368317_at) shows almost identical results.
3. 1390457_at: mitogen activated protein kinase kinase kinase 13 (Map3k13)
Thirty-five % of the target sequence is 75% identical to Chinese hamster Map3k13 (2e-19). This probe set used to be assigned to a non-coding RNA, LOC102546784, which contained a partial sequence of a non-coding region of rat Map3k13, but is presently deleted from the database. Although it is incomplete at present, we have found a sequence containing both the target sequence of 1390457_at and one of the exon of Map3k13 (unpublished observation).
4. 1373468_at: cyclin F (Ccnf)
The target sequence is about 70% identical to cyclin F of several rodents (1e-17 to 1e-42). Moreover, 9 out of 11 probes match with 1382370_at (A), which is officially annotated as cyclin F.
Extraction of PPARα agonist-responsive genes in the kidneyTGPJ collected transcriptome data of in vivo kidney for 50 out of 150 chemicals. As the selection was mainly based on renal toxicity, only CFB was included in the 50 chemicals as a typical PPARα agonist. Thus the CFB-responsive probe sets were extracted in the same way as liver, except that ratio to control >1.5 and effect size >1.0, because it was impossible to take common probe sets among multiple fibrates.
Search for PPAR-responsive element (PPRE) in the upstream of each extracted genePPARα forms a heterodimer with Retinoid X receptor (RXR) and they bind to the DNA binding element, called PPAR response element (PPRE)/Direct Repeat 1 (DR1), and the consensus sequence of PPRE/DR1 is composed of 2 sequences of “AGGTCA” with unknown redundancy separated by a single nucleotide spacer, where PPARα and RXR bind to the 5’ and 3’ AGGTCA sequences, respectively (Schoonjans et al., 1996). As RNA polymerase works bi-directionally from the promoter region (Lloret-Llinares et al., 2016), its everted sequence, TGACCTNTGACCT in the complementary chain, is also PPRE. More recently, Tzeng et al. (2015) proposed WAWVTRGGBBAHRGKTYA (W=A or T, V=A, G, or C, R=A or G, B=G, T, or C, H=A, C, or T, K=G or T, Y=C or T) as the optimal PPARα/RXR heterodimer binding sequence. After the extraction of PPARα-responsive genes, the above sequences were searched in their 5’ upstream region based on the genomic database of NCBI as follows.
a) AGGTCANAGGTCA: Reviewing the reported PPRE, usually 2 or 3 substitutions are allowed and the distance from the first exon is as long as several thousand bases (e.g., Nagasawa et al., 2009). The incidence of this sequence is ca. once per 4.5 x 105 bases and if one or two nucleotide substitutions are allowed, it increases to once per 453,438 or 26,588 bases, respectively. Although gene promoters often locate far from the transcription start sites (Sanyal et al., 2012), three nucleotide substitutions would be less meaningful, because the incidence increases to once per 2,553 bases. According to previous studies, PPREs are located within 5,000 bases upstream or downstream of the transcription start site in most cases (Fang et al., 2016). We then searched PPRE up to 10,000 bases upstream from the first exon allowing one or two substitutions, and when no such sequences were found, the one with three substitutions up to 2,553 bases upstream was recorded. When all else failed, the complementary sequence, TGACCTNTGACCT, was searched based on the above criteria.
b) WAWVTRGGBBAHRGKTYA: Although the incidence of this sequence is calculated to be ca. once per 1.3 x 107 bases, some substitutions of nucleotides are allowed with different contribution (Tzeng et al., 2015). In this case, we simply recorded a sequence with the least substitutions within 10,000 bases.
Analysis of correlation of induction rates between PPARα agonist-responsive genes and their neighbor genesA gene and its neighboring genes are categorized as shown in Fig. 1, i.e., A: downstream on the same strand, B: head to head on the other strand, and C: tail to tail on the other strand. As for the 249 probe sets, which showed at least 2.5 fold average induction by three fibrates in vivo liver, their neighboring genes were identified in the NCBI database and the distance between them (the distance between the last exon of the first gene and the first exon of the second gene for A, and the distance between the first exon of the two genes for B and C) was recorded. In the next step, it was examined whether induction of a gene affects that of its neighboring gene. To do this, drugs with various PPARα agonistic potencies were selected from TG-GATEs. As fibrates, CFB, WY, GFZ, FFB, and the drugs that were reported to have PPARα-stimulating activity, VPA, AM, BZD (Szalowska et al., 2014), BBr (Kunishima et al., 2007), SST (Seo et al., 2008), RGZ (Edvardsson et al., 1999) were selected. As it was reported that NSAIDs, at least some of them, showed PPARα stimulation (Lehmann et al., 1997), all NSAIDs in the database, i.e., ASA, IBU, MEF, SS, IM, NPX, BDZ, NPAA, DFNa were included. In addition to these, CMP, TBF, AAF were included because they were found to induce Acot1, a typical PPARα agonist-responsive gene in the rat, in a preliminary search. As for the total of 22 drugs, Pearson’s correlation coefficient of induction rates of neighboring gene pairs was calculated where the ratio to control was converted to its natural logarithm. When the values in the control and treated groups are both under the detectable level (“absent” or “marginal” call), the induction rate shows an extreme value that disturbs correlation analysis. To avoid such a case, when the “present calls” indicated by the system were less than 3 out of 6 (3 for each control and treated), the data were excluded from the calculation. The pairs with significant (p < 0.05) correlation or inverse correlation were selected, and the genomic location of their human orthologue was searched for in the NCBI data base. When the gene name was not determined in either of the gene pairs, it was assigned to the one without a gene name in the similar position exclusively when the other one of the pair has the same gene name. Finally, the scattered plot was drawn for the correlation coefficient with the distance between the neighboring pair on the rat genome. These procedures indicate the relative position of a gene pair with correlated expression by PPARα agonists in both rat and human orthologue.
In order to select the genes that are reproducible and largely induced by PPARα agonists, 4 fibrates, i.e., CFB, FFB, GFZ, WY, were selected at first, but we consider that the dose level of GFZ was too low to observe inductions of typical PPARα agonist-responsive genes. For example, Fig. 2A-D shows the effects of these 4 fibrates on carnitine palmitoyltransferase 1b (probe sets ID: 1367742_at) of rat liver in vivo. This gene was reported to be induced by fibrates (Silvestri et al., 2006) and at least 9-fold induction was observed at 24 hr after administration of the high dose of CFB (A), FFB (B), WY (C), whereas the fold increase was only 2.4 by GFZ (D). Although clear induction was observed by GFZ after repeated administration (data not shown), we were obliged to exclude this drug from the first extraction step so as not to overlook important responsive genes.
Focusing on in vivo data of 24 hr after single administration of high dose, the induction rate was calculated and sorted. In order to obtain reliable data, probe sets showing ratio to control > 1.5 and effect size > 0.67 were selected after filtering by PA-calls. Extracted probe sets were 3116, 2767, and 3127 from CFB, FFB, and WY, respectively, and the intersection, common ones among three fibrates, was 532. Although there were many probe sets showing significant decrease by fibrates, they were excluded from the present analysis, considering the fact that PPARα forms transcriptional machinery.
In the next step, PPARα agonist-responsive genes were extracted from in vitro data of 24 hr after treatment with high concentration (300, 30, and 150 μM for CFB, FFB, and WY, respectively). As too many known inducible genes dropped out if the cut-off level was set to the same as in vivo, probe sets with ratio to control > 1.3 and Ef value > 0.3 were selected resulting 1634, 1891, 2834, from CFB, FFB, and WY. The intersection, common ones among three fibrates, was only 155.
Table 2A shows the highest 50 probe sets of in vivo liver by sorting the mean induction ratio of three fibrates. Table 2B shows the highest 50 probe sets of in vitro by sorting the mean induction ratio of three fibrates as well as the difference from the rank in vivo, where the probe sets not selected in vivo were uniformly ranked as 533. In Table 2A and B, those commonly selected between in vivo and in vitro are shaded. As the number of probe sets in vivo is 3 fold larger than that in vitro, it is reasonable that there are many unshaded ones in Table 2A. If the induction of gene expression by fibrates occurred similarly in vivo and vitro but differed only in their extent, the uncommon probe sets should be unevenly distributed to the low ranking part of the in vivo list and the difference in the ranking should be small . However, in Table 2A, which contained the top 10% of the extracted 532 probe sets, only 24 out of 50 are selected in vitro, indicating that there are many genes induced by fibrates exclusively in vivo. Looking at the difference between in vivo/vitro ranking in Table 2B, there is no trend in the difference. In contrast to Table 2A, 40 out of 50 probe sets in Table 2B are selected in vivo, meaning that 80% of the top 50 probe sets in vitro are inducible in vivo.
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/45_449-t2A.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/45_449-t2B.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
The present study does not focus on the in vitro data in detail, while there are a few explainable cases of genes that have high rank in in vitro but not selected in the in vivo list. For example, Fig. 3 shows the effect of WY on the expression of fatty acid binding protein 1 (Fabp1, 1369111_at) normalized by that of peptidylprolyl isomerase A (Ppia, 1398850_at). Although significant changes were noticed at several points, the expression of this gene was not largely mobilized by WY in vivo (Fig. 3A) possibly because it was already maximally expressed without agonist (or possibly with endogenous agonists from food) allowing no further stimulation. In contrast, the basal level drastically decreased with time in vitro such that the stimulatory effect emerged after 24 hr of culture (Fig. 3B).
Reviewing probe sets with high ranking exclusively induced in vivo in Table 2A, most of them appear to be unrelated to lipid metabolism, e.g., “cholinergic receptor, nicotinic, alpha 2 (Chrna2)”,”histidine decarboxylase (Hdc)”, “aquaporin 3 (Aqp3)”, “NIM1 serine/threonine protein kinase (Nim1k)”, “androgen-induced 1 (Aig1)”, “proline rich Gla4 (Prrg4)”, “calcium channel, voltage-dependent, alpha2/delta subunit (Cacna2d3)”, “glutaminyl-peptide cyclotransferase (Qpct)”, “ADAM metallopeptidase domain 22 (Adam22)”, “collagen-like tail subunit (single strand of homotrimer) of asymmetric acetylcholinesterase (Colq)”, and “mitogen activated kinase kinase kinase 13 (Mp3k13)” (ranking 2 and 12, 4, 7,9,11,14,15, 17, 18, 42, respectively).
In order to examine whether the induction of genes unrelated to lipid metabolism is a feature of fibrate itself, or dependent on the target organ, PPARα agonist-responsive genes were extracted from the data of kidney in vivo. Unfortunately, drugs tested for both liver and kidney were restricted to 50 and only CFB was included as a fibric acid. Extraction of the data 24 hr after the administration of high dose of CFB by ratio to control > 1.5 and effect size > 0.67 after filtering by PA-call resulted in 334 probe sets. These contained a considerable number of unreliable ones with respect to dose- and time-dependency, and thus the hurdle of effect size was increased to 1.0 to get 310 probe sets. The top 50 probe sets are shown in Table 2C. Of the 310 probe sets from kidney, only 64 are common with 552 from liver, whereas 23 of the top 50 in Table 2C (shaded) are common with liver. Regarding the remaining 27, most of them appear to be related to lipid metabolism. Most of the genes unrelated to lipid metabolism selected in the liver were not extracted from the kidney. Of the 10 genes listed above, only Chrna2 and Colq were selected in the kidney.
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/45_449-t2C.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
Search for PPAR-responsive element (PPRE)-like sequence in the upstream of each extracted gene
We were interested in whether the difference in the responsiveness of the genes to PPARα agonists could be explained by the PPRE-like sequence for each gene. We surveyed PPRE-like sequence for the genes in Table 2A and evaluated their concordance of the sequence and the distance from the first exon (Table 2A, right columns). Of the 50 probe sets in the Table, five were duplicated and one had no annotation, with the result that 44 genes were analyzed. For the consensus sequence “AGGTCANAGGTCA”, the distance from the first exon by the number of bases, and the mismatched bases in the sequence are indicated by lower-case letters. When a sequence with less than 3 mismatches are absent within 10,000 bases and that with 3 mismatches are absent within 2553 bases, the result of its complementary sequence, TGACCTNTGACCT is shown. All of the genes in Table 2A have PPRE fulfilling the condition except Krtap1-5 (No. 28), which has a very short intron (2054 bases) between upstream gene (Krtap1-1) and only has a sequence with 4 mismatches. The other consensus sequence “WAWVTRGGBBAHRGKTYA” was also surveyed and is shown in Table 2A. For all of the 44 genes, including Krtap1-5, sequence with 4 or less mismatches was found within 10,000 bases. However, neither the concordance nor the distance from the first exon of the PPRE-like sequence appears to be related to the extent of induction by PPARα agonists.
The condition of PPRE for AGGTCANAGGTCA employed in the present study allowed 2 to 3 mismatches within several thousand bases. This was practically necessary so as not to exclude well-known PPARα agonist-inducible genes in Table 2A (e.g., Acot3, Cpt1b, Ehhadh, etc.), but this condition appears to be too lenient for selection. For example, α2 subtype of nicotinic receptor (Chrna2) is extensively induced by PPARα agonists (Table 2A, No.2 and 12) but none of other subtypes are detected in the liver with or without stimulation. However, all Chrna subtypes have PPRE-like sequence, i.e., Chrna1, 2, 3, 4, 5, 6, and 7 has AGcTCAGctGTCA at -79, AGGggAAAGGgCA at -817, AGGTCAAgGGaCc at -972, AGGcCAGAGGgCA at -1305,AGtTCAGAGGcCA at -2129,gGGTgACAGGaCA at -106, gGGTCAGAGacCA at -145, respectively. Moreover, similar sequences could be easily found in many genes. For example, ATPase, H+ transporting, lysosomal accessory protein 1(Atp6ap1, 1367873_at) showed constant expression with or without stimulation, but it has AGGTCAGAaGaCA at -2422. The expression of Apolipoprotein CIII(Apoc3, 1370009_at)was markedly suppressed by PPARα agonists but it has AGGcCAGAGGTCA at -1660. These results clearly indicate that a simple survey of PPRE consensus sequence cannot explain the observed gene induction by PPARα agonists.
Expression of cMyc and other proliferation-related genesAn attractive hypothesis for the carcinogenesis of PPARa agonists in rodents proposes that stimulation of PPARα downregulates let-7c and that increases cMyc, leading to hypertrophy and carcinogenesis (Shah et al., 2007; Gonzalez and Shah, 2008). However, cMyc was contained neither in the in vivo 532 probe sets, nor in the in vitro 155 probe sets. Then we checked the data for cMyc probe set, 1368308_at, with single and repeated doses of CFB, FFB, and WY, shown in Fig. 4. In the case of CFB, no obvious induction of cMyc was seen, whereas both FFB and WY showed dose-dependent, obvious induction after 7 days and longer, but not on or before 3 days. As the gene extraction was done for data 24 hr after high single dose, cMyc was excluded from the gene list. As for the in vitro data, ratio to control 24 hr after treatment of high concentration of CFB, FFB, and WY was 0.87, 0.96, and 0.88, respectively, showing no induction of cMyc.
It is well known that cell proliferation increases before induction of hepatocarcinoma by fibrates (Shah et al., 2007; Yang et al., 2007). It is thus of interest when cell proliferation is activated by the administration of PPARα agonists. As there are no genes directly related to cell proliferation in the top 50 genes (Table 2A, except personally annotated cyclin F) of extracted 532 probe sets in vivo, we expanded the analysis to the 250 probe sets showing more than 2.5 fold mean increase by 3 fibrates in vivo liver, and found that many cell proliferation related genes were included. Table 3 shows 57 probe sets selected by the key words directly related to cell proliferation with the guide of GO biological process terms supplied by Affymetrix. Of these, dose- and time-dependency of the induction of three genes, geminin (Gmnn, 1378056_at), topoisomerase IIα (Topo2a, 1372186_a_at), and marker of proliferation Ki-67 (Mki67, 1374775_at), which are established gene markers of proliferation (Whitfield et al., 2006), are shown in Fig. 5A-I. The time course of the induction of these genes by the three fibrates are similar; the maximal induction occurred at 24 hr after single administration and it decreased with repeated administration. It is clear from the data that cell proliferation in rat liver is stimulated by the fibrates 24 hr or earlier after single administration, in contrast to cMyc induction.
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/45_449-t3.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
In order to see organ specificity, cell proliferation-related the 57 probe sets in Table 3 were compared with the ones upregulated in kidney. Of the 310 probe sets extracted from kidney, only 5 (Ccna2, Top2a, Ska1, Brca1, and Cdca5; their ranking was 128, 59, 186, 310, and 139, respectively) were included in the 57 probe sets. These 57 probe sets were also compared with the ones upregulated in vitro hepatocytes. Surprisingly, none of the 155 probe sets was included in the 57 cell proliferation-related probe sets. Inversely, GO analysis of in vitro 155 probe sets revealed that only 3 (Nrli3, Klf11, and Tp53I13; their ranking was 44, 118, 129, respectively) were directly related to cell proliferation and they were not extracted from in vivo data.
In order to elucidate the reason why the cell proliferation-related genes were not induced by PPARα agonists in vitro hepatocytes, time course of the expression level of these gene was analyzed using peptidylprolyl isomerase A (Ppia) as an internal standard to compare the absolute expression levels between in vivo rat liver and in vitro rat hepatocytes. It was found that the expression pattern of the gene in vitro was classified into three types. The first one is that the expression without agonist gradually increased with culture time and reached its apparent maximum at 24 hr, estimating that the expression with the high dose of agonist in vivo at 24 hr is its maximum level. Chaf1a, Gmnn, and Gtse1 belong to this type. As a representative example, expression of Gtse1 in vivo (Fig. 6A) and in vitro (Fig. 6B) treated with FFB are shown as normalized by Ppia expression. It is clear that no stimulation by FFB was observed in vitro at any time points in contrast to the dose dependent increase at 24 hr in vivo. The second one is that the expression was at the apparent maximum without agonist from the beginning of culture. Kif22, Cdkl2, Kif7, and Supt3h belong to this type, and expression of Kif22 in vivo (Fig. 6C) and in vitro (Fig. 6D) treated with FFB are shown. Dose-dependent increase of expression was noted in vivo at 24 hr, whereas the expression appeared to be at its maximum throughout the culture time with or without FFB. The third one is that the remaining 46 (4 genes are duplicated) showed a sharp rise at 24 hr of culture without agonist. As a representative example, the expression of Kifc1 in vivo (Fig. 6E) and in vitro (Fig. 6F) treated with FFB are shown. In this case, control value at 24 hr appears too high to detect stimulation by FFB, whereas that at 3 and 8 hr in vitro appears low enough, but no stimulation by FFB was observed. It is noticed that no stimulation was observed before 24 hr in vivo, too. From Fig. 6A-F, it appears that expression of the genes in Table 3 reached their maximum in vitro at the time point when the stimulation by PPARα agonist is detectable in vivo. These patterns are very similar in other agonists, CFB and WY (data not shown).
Again it appears to be difficult to explain the difference in the expression of cell proliferation-related genes by the analysis of PPRE-like sequence. For example, Cyclin B1 (1370345_at, averaged ratio to control is 3.74) has gGGaCAGAGGaCA at -1236 bases from the first exon, whereas cyclin M2, which showed constant expression with or without agonists, has AGGTtAGAGGaCA at -1499.
Expression of neighboring genes of PPARα-responsive genesGene expression in eukaryotic cells accompanies modification of histones, and consequently, expression of the neighboring genes is affected (Sémon and Duret, 2006). Recent study revealed that RNA polymerase works bi-directionally from the promoter region and transcription of two genes located as B in Fig. 1 should start at once and is called promoter upstream transcription (PROMPT) (Lloret-Llinares et al., 2016). It is of interest whether the lipid metabolism-unrelated genes were transcribed, as they are located next to the proper targets of PPARα agonists. We chose 249 probe sets showing more than a 2.5 fold mean increase by 3 fibrates in vivo liver and their neighboring genes were identified in the rat genome, and categorized into A, B, and C (Fig. 1). In these 249 probe sets, there were many cases such that the gene corresponding to the probe set was not identified as a true gene, both probe sets of the neighboring gene pair were included in 249 sets, and the probe sets for the neighboring gene were not designated as GeneChip, resulting in that 127 pairs for type A, 107 pairs for type B, and 105 pairs for type C were available for further analysis.
If the expression of the neighboring genes affect each other, the extent of their induction by PPARα agonists is expected to be correlated with their agonistic potency. Thus, 22 drugs with various potencies of PPARα agonist were selected, and the log e of mean induction rate at 24 hr after a single administration of the high dose of each drug was calculated for the neighboring pair, and their correlation coefficient was calculated. As a whole, the correlation coefficient distributed higher than 0.9 to –0.88 in the types A to C, showing no trends. Fig. 7A, B show two extremes of type A. Probe set 1388211_s_at does not discriminate Acot1 from its neighboring gene, Acot2, and subsequently, the induction rates of Acot1 and 2 by the 22 drugs were highly correlated (r = 0.97). Excluding this case, the highest correlation was obtained between Acot4 and 3 (Fig. 7A, r = 0.91). In contrast, Melk and its downstream gene Zcchc7 showed a clear inverse correlation (Fig. 7B, r = –0.64). Fig. 8A, B show two extremes of type B. The highest correlation was observed between Hadha and Hadhb (Fig. 8A, r = 0.94) whereas Pex19 and Ncstn showed inverse correlation (Fig. 8B, r = –0.81). Fig. 9A, B show two extremes of type C. The highest correlation was between Esco2 and Pbk (Fig. 9A, r = 0.96) whereas Aqp3 and Nfx1 showed inverse correlation (Fig. 9B, r = –0.88).
In order to pick out the pairs with potential interaction, the pairs with statistical significance (p < 0.05) were selected, and 22 from Type A, 24 from Type B, and18 from type C were obtained. They are summarized together with the distances on the rat genome, and with the relative location of their human orthologues on the human genome (Table 4). In order to see the relation between their co-expression and genomic distance, correlation coefficients and genomic distances were also shown as scattered plots in Fig. 10A, B, and C.
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/45_449-t4.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
Of the genes belonging to type A, Acot1, 2, 3, and 4 appear to form an obvious gene cluster, and they were all induced with high correlation, suggesting that a sequence of genes with similar biological function is read through at once. However, this cannot be generalized because Plin5 (rank 19 in Table 2A; mean ratio to control is 9.5) and Fbp2 (rank 20 in Table 2A; mean ratio to control is 9.3) both have a similar gene directly downstream of each, Plin4 and Fbp1, but the latter were not induced at all (data not shown).
In type B, the positive correlation coefficient tends to be larger as the distance of the gene pair is closer. However, there are several pairs showing negative correlation coefficient with varied distance. Except for the Cyp4a1/Cyp4a3 pair, which exists in a huge Cyp4a gene cluster, pairs showing high correlation coefficient mainly consist of a pair with at least one gene unrelated to lipid metabolism. Of these, the following are the pairs of a typical PPARα agonist-responsive gene and a lipid metabolism-unrelated, but physiologically interesting gene, i.e., Crat/Ptpa, Chnra2/Ephx2, and Map3k13/Ehadh.
In type C, no relationship was seen between the distance and the correlation coefficient. The highest correlation was observed between Esco2 and Pbk, both of which are unrelated to lipid metabolism, i.e., chromosome segregation and the MAPK cascade, respectively.
If these pairs contain genes that are involved in rodent-specific hepatotoxicity, it is conceivable that their relative position is different from that in the human genome. As shown in Table 4, more than two thirds of pairs in rats have their orthologue with the same relative positions in the human genome (indicated as “=” in the column). For three pairs of Acot and a pair of Cyp, they could be considered to be the equivalent gene position in both species. For the pair Ccna2/RGD1307100, it appears that LOC102724158 exists in the human genome corresponding to RGD1307100, but a similarity between these two genes could not be found. For 16 other pairs, there are one or more genes that exist between the human genome. These 17 pairs (16 and Ccna2/RGD1307100) are indicated as open circles in Fig. 10. It is noticeable that some pairs are exceptionally distant from each other (Vwa8/Mtrf1 and Melk/Zcchc7 in type A, Ccna2/RGD1307100 and Aldh1a1/Anxa1 in type B, Hsdl2/Inip and Aqp3/Nfx1 in type C), and they are all categorized in this group. It is conceivable that there could be gene(s) between these pairs in rat genome also. There is no case that human orthologues exist on different chromosomes, or extremely far from each other.
The above-mentioned correlations were observed at one time point, 24 hr after the single high dose. There was not necessarily any similarity in the dose- or time-dependency between the pair. For example, the time- and dose-dependencies of Chrna2/Ephx2 pair, and those of Map3k13/Ehhadh pair, stimulated by WY, are shown in Fig. 11A-D. In the case of the former pair (Fig. 11A, B), the induction of Chrna2 increased with time by single administration and it further increased with repeated administration, whereas the induction of Ephx2 reached its peak at 9 hr after single administration and no further increase was observed in repeated administration. In the case of the latter pair (Fig. 11C,D), the induction of Map3k13 reached its maximum at the low dose on and before 9 hr, and it became dose-dependent on 24 hr after single or repeated administration, whereas the induction of Ehhadh reached its maximum by the low dose at all time points.
It is already established that PPARα agonists induce hepatocyte proliferation leading to carcinoma in the liver of rodents but there is no such risk in humans based on much clinical experience, without any rational explanations (Desvergne and Wahli, 1999). It has been stated that the number or affinity of human PPARα is lower than that of rodents (Holden and Tugwood, 1999; Tateno et al., 2015). However, if it is an issue of dose effect, there should be a low but certain risk. For all-or-none response, there should be gene(s) that are induced and carcinogenic in rodents but not in humans, insofar as PPARα works as a transcription factor. According to Gonzalez and Shah (2008), cMyc is constantly destabilized by an miRNA; let-7C, and WY suppresses the expression of the miRNA in a rodent-specific PPARα-dependent manner to induce overexpression of cMyc, leading to hepatic hypertrophy and carcinogenesis. This indicates that the search for the causal genes should be focused on cMyc-related genes. However, our analyses have failed to support this. In studies to identify non-genotoxic hepatocarcinogens, including PPARα agonists using support vector machine, cMyc was not selected in feature genes in any cases (Uehara et al., 2008, 2011; Yamada et al., 2013). In the present study, cMyc was not selected as a PPARα agonist-responsive gene either. Reviewing the time- and dose-responses of cMyc to the three agonists, only FFB and WY showed significant induction 7 days or later with repeated administration and no change after single administration. Of course, cMyc plays a main role in cell cycle regulation and is essential for cell proliferation and carcinogenesis (Qu et al., 2014). In fact, our previous analysis showed that cMyc itself was not included in the feature genes of the support vector machine, but instead is part of the genes connected to Myc by pathway analysis (Uehara et al., 2011). Therefore, it is reasonable to postulate that there would be species-specific changes in candidate genes other than cMyc at an earlier stage, i.e., 24 hr after single administration. Considering the fact that many cell-proliferation-related genes (including the established marker genes) were upregulated within 24 hr after a single administration of fibrates, the list in the present study could be useful for further study.
In the analysis of in vitro data, the responsiveness of cultured hepatocytes is generally much lower than in vivo liver, and a much smaller number of reliable genes could be extracted. This is understood from the well-known observation that hepatocytes rapidly dedifferentiate during conventional culture (Berry et al., 1997). In a few cases, induction was more marked in vitro than in vivo, where the gene was extensively expressed or changed with a large circadian rhythm without agonists in vivo and the constitutive expression decreased with time of culture such that stimulation by agonists emerged.
The most striking difference between in vivo and in vitro was that the high rank part of the list of PPARα agonist-responsive genes in vivo contains items that appear to be unrelated to lipid metabolism and are mostly not inducible in vitro. It is particularly worthy of note that all of the cell proliferation-related genes extracted from in vivo liver were all absent in the list obtained from in vitro hepatocytes. Similar results were obtained previously (Parzefall et al., 2001; Tamura et al., 2006), and the reason has been attributed to TNFα released from Kupffer cells (Parzefall et al., 2001), or to undefined responses to PPARα in non-parenchymal cells (Qu et al., 2010), which are excluded in the primary culture system. According to the present results, however, the apparent lack of stimulatory effect of fibrates in vitro is not due to loss of responsiveness but to the increase in cell proliferation during culture, at least in part. This could be the result of the cell separation process, inevitable for in vitro experiments.
As the hypolipidemic activity of PPARα agonists is common to many species, rodent-specific toxicity would be due to activity unrelated to lipid metabolism. An attractive hypothesis is that the causal genes are included in these cases unrelated to lipid metabolism. For this purpose, one would propose a strategy to see if the PPARα responsive element exists over the whole genome. However, searching for promoter sequences in the genome is not a successful strategy in general, because the length of promoter sequences is so short that they are most frequently found by chance. The PPARα responsive element, AGGTCA-n-AGGTCA, is relatively long, but because of the ambiguity in the sequence, the usual way is to search upstream of the gene after the target gene is filtered (Tompa et al., 2005). It has been revealed that the transcription controlling region is not restricted to the vicinity of the target gene but often exists far from the target (Sanyal et al., 2012). Moreover, from the fact that expression of genes with similar PPRE in the same organ are differently stimulated, and that a gene with PPRE is differently stimulated in different organ, it is clear that there should be more important, unknown factors in the control of gene expression. Under the present technology, the best way is to drill down the candidate genes to start with.
There are many examples where expression of a gene affects a neighboring gene, such as Hox genes that form a gene cluster, and a group of neighboring genes is immediately regulated (Kondo and Duboule, 1999). Since stimulation of PPARα involves acetylation of histone (Asano et al., 2010), the expression of not only the target gene but also the neighboring genes could be affected. Furthermore, recent studies revealed a phenomenon called “promoter upstream transcription” (PROMPT) (Sigova et al., 2013). Contrary to previous common sense, the direction of transcription is not unidirectional along with the target gene, but randomly bidirectional, i.e., RNA polymerase runs either direction with 50% probability from the promoter region. In the present study, the pairs of Fig. 1B should be transcribed on both sides of the starting point. In general, either or both of these RNAs would be degraded before or after the maturation to mRNA, sequestered as mRNA from translation, or translated to protein. In addition, it is postulated that these RNAs include non-coding RNAs that have some biological roles (Sigova et al., 2013; Almada et al., 2013; Lloret-Llinares et al., 2016). As expected, the induction of a gene by PPARα agonists appeared variously affecting its neighboring genes. In the protocol for GeneChip, oligo dT primer is used for the first reverse transcription step, and subsequently mRNAs were chosen to be measured. According to PROMPT, both genes in the pair belonging to B-type should be transcribed to the same extent, and consequently, in the pairs showing no or inverse correlation, the mRNA of the partner gene is expected to be degraded. This would be also true in the pair of genes showing high but various correlations, where the fate of mRNA is considered to be varied case to case. At present, comprehensive understanding of PROMPT is immature, with only a few elucidated examples. Gene expression profiles in rat liver by PPARα agonists would be a good model to investigate the interaction of gene expression, because the number and the induction rates of genes are exceptionally high among many drugs tested (Tamura et al., 2006; Uehara et al., 2008).
Table 4 summarizes the gene pairs affected by PPARα agonists and candidate genes involved in the species-specific hepatotoxicity that we expect are included. With naive speculation, a particular gene is induced by a PPARα agonist together with its true target leading to cell proliferation or tumorigenesis, but not in human because there are other gene(s) in between or their location is in different chromosomes that interfere with their co-expression. In this case, the candidate genes are the partner of each pair other than that with “=” in Table 4, but there are some problems. Reviewing the candidate genes whose human orthologue pair has gene(s) in between, many of the pairs are separated by a relatively long distance. This does not necessarily exclude the possibility that the gene structure is different between rat and human, but there could be one or more genes to be discovered in the wide blank space in the rat genome, because deciphering rat genome is quite undeveloped compared with human or mouse.
There is another technical limit. Although more than 30,000 probe sets, which are more than the number of total genes, are designed for GeneChip, they are duplicated and important genes are still lacking. In the analysis of neighboring genes, no probe sets are designed for more than a half of the selected genes. These facts suggest that there could be still many important genes overlooked.
Under the present situation, the best way would be to consider all the genes in Table 2A and Table 3 as candidates, since the apparent interaction of the expression of paired genes is not completely dependent on their location. The most attractive idea is to focus on the genes apparently unrelated to lipid metabolism. Focusing on the induction rate, the physiologically interesting ones are Aqp3 and 7, Chrna2, Colq, Hdc, and Cacna2d3. Focusing on the mechanism of hepatocarcinogenicity, genes related to intracellular signal transduction, such as Nim1k, Pbk, Ccdc34, and Map3k13 in Table 2A and the 57 genes directly related to cell proliferation in Table 3 are worth investigating. For cell-cell interaction or cytoskeleton related ones such as Cep128, Adam22, Krt23, Krtap1-5, Rab30, Rhof would be also attractive subjects for future study to examine details.
In the present study, only transcription was analyzed, and the important issue is to see if they are translated and working as protein. Furthermore, it has been revealed that the amount of un-translated RNA is much larger than that of translated RNA, and they have certain physiological roles (Cech and Steitz, 2014; Lloret-Llinares et al., 2016), suggesting that study for the role of the RNA itself is also important. From another point of view, it is reported that PPARα agonists mediate nongenomic signaling, such as activation of mitogen-activated protein kinase pathways (Gardner et al., 2005). Further investigation is clearly necessary, and the genes extracted in the present study would be a clue for a future study to elucidate the mechanism of gene expression as well as the toxicity of PPARα agonists.
We thank Prof. Nobuyosi Akimitsu (Isotope Center, The University of Tokyo) and Prof. Atsushi Ono (Okayama University) for fruitful discussion. We thank Ms. Karin Ozaki for data analysis and Mr. John H. Jennings for editing manuscript.
Conflict of interestThe authors declare that there is no conflict of interest.
