Ochratoxin A induced early hepatotoxicity: new mechanistic insights from microRNA, mRNA and proteomic profiling studies

The mycotoxin ochratoxin A (OTA) is found widely in agricultural commodities. OTA can induce various toxicities. In this study, rats were gavaged with OTA for different weeks. Then, the expression of microRNAs, mRNAs and proteins were measured in the rat livers treated with OTA for 13 weeks. Our sequencing data suggests that the medial and the high doses of OTA exert different effects on livers. Five distinctive pathways were induced after OTA treatment as collectively demonstrated at miRNA, mRNA and protein levels. Two (primary bile acid biosynthesis and metabolism of xenobiotics by cytochrome P450) are directly associated with liver damage, whereas the remaining pathways (arginine and proline metabolism, cysteine and methionine metabolism and PPAR signaling pathway) cause metabolic disease. This study reveals OTA-induced early hepatotoxicity for the first time by combining multi-omics methods. The novel metabolic pathways may contribute to the pathogenesis of metabolic diseases later in life.

including induction or maintenance of cell fate in normal, stem and cancerous cells. With these functions, miRNAs have been extensively studied as biomarkers for the diagnosis and therapy in diverse range of diseases, as well as for their impact on gene regulation. The miRNA expression profiling is not only helpful in identifying miRNAs capable of regulating a large range of biological processes, but also can be systematically utilized to study the gene regulation, especially when miRNA measurements are combined with mRNA profiling and other genome-scale data.
Within recent years, the multiple omics profiling technologies are considered powerful to acquire corroborative evidence for hepatotoxicity study. Through a cross-omics method, it is demonstrated that the major metabolic consequences of acute toxicity of aflatoxin B1, one of the well-known hepatotoxic mycotoxin, are associated with gluconeogenesis and lipid metabolism pathways 9 . Moreover, inorganic arsenic methylation in mouse liver may be saturated under chronic exposure to the element based on an integrated-omics approach 10 . The toxicity mechanism of another hepatotoxic compound EMD 335823 was also elucidated by using modern omics technologies, from which the liver pathologies are linked to deregulated PPARa signaling, glucose metabolism, and fatty acid metabolism 11 . Furthermore, omics analyses on several hepatotoxicants have enhanced preclinical safety assessment and helped to identify the nature of injury. In particular, Wnt/b-catenin signaling has been identified as the biomarker of human hepatocellular carcinoma and gastrointestinal disease 12 . Altogether, cross-omics methods indeed have enhanced the understanding of liver metabolic regulation of exogenous chemicals. The exploration of OTA hepatotoxicity also can be better understood by integration of multiple omic platforms (Fig. S1). On the other hand, a whole body model is superior to an in vitro model to better reflect the metabolic and physiological circumstances.
In order to better understand the early hepatotoxic response to OTA exposure in male F344 rats, we herein integrated miRNA, mRNA and proteomic profiling approaches to identify the biological pathways of early OTA-induced hepatotoxicity.

Result
Health status of rats gavaged with OTA. After the administration of OTA, body mass of the high-dose rats reduced about 7% (p , 0.05) at 9 and 10 weeks 13 , but gradually recovered thereafter. The liver weights were comparable between OTA-treated and control rats at 2, 4 and 13 weeks (Fig. S2). These results are consistent with a previous toxicity study of orally administering OTA for 13 weeks 5 .
LDH of a and c groups, and AST of c group were lower in the highdose than the control rats, which were consistent with previous reports 14, 15 . In contrast, ALB activity of c group was slightly higher in high-dose group than in control rats 13 . The liver damage in highdose group is more serious than that in medial dose group at 13 weeks. Hence, we focused on the high-dose group in the subsequent study.
Analyses of a complete gross necropsy and microscopic anatomic pathology demonstrated that rat livers treated with OTA for 2, 4 and 13 weeks were comparable to those without OTA treatment. Bile duct hyperplasia appeared in rat liver treated with OTA for 26 weeks (data not shown). We concluded that OTA did not apparently induce pathological damage in the livers of the rats in the 13 weeks, in which point the cross-omics analysis of rat livers was conducted.
OTA differentially affects miRNA expression profiles in the livers. To evaluate miRNA expression profile in livers of rats treated with OTA, an Illumina Hiseq 2000 platform was utilized. A total of 304, 301 and 304 known mature miRNAs were identified in control, medial-dose and high-dose groups, respectively. Of these, 272 mature miRNAs were found to be expressed in all the three groups. We also found that 17, 11 and 11 miRNAs were confined to control, medial and high groups respectively.
Hierarchical clustering results depicted that miRNA expression altogether was similar between control and high-dose groups, but they were differentially expressed compared to medial-dose group (Fig. 1). Our sequencing data indicates that OTA-induced hepatotoxicity is differentially mediated in medial-dose and high-dose groups.
Furthermore, an additional miRNA-target prediction was conducted with at least three different programs 23 and pathway enrichment analysis was carried out. The pathways resulting with ''at least two algorithms'' and ''at least three algorithms'' were comparable in which ''pathways in cancer'' and ''MAPK signaling pathway'' were highly significantly enriched in all OTA-treated groups.
OTA differentially affects mRNA expression profiles in the livers. Illumina HiSeqTM2000 was used for identifying the differentially expressed mRNAs in livers of OTA-treated and control rats. The differentially expressed genes between two groups were screened based on DESeq analysis, and the controlling FDR was set to 5% level of significance. We identified 28 up-regulated and 17 downregulated genes in high-dose compared to control group. Six genes were highly elevated, while three genes were significantly suppressed in medial-dose group compared to control group. Only one gene was differentially expressed in high-dose compared to medial-dose group. These differentially expressed mRNAs (Log 2 fold change . 1, p , 0.05) are shown in Table 1. The KEGG pathway enrichment analysis on these differentially expressed mRNAs resulted in 50 and 24 enriched pathways in high-dose and medial-dose groups, respectively. Furthermore, there were 20 overlapping pathways, among which ''amino acid metabolism'', ''xenobiotics biodegradation and metabolism'', ''energy metabolism'', and ''environmental information processing'' were enriched.
Overlapping pathways between deregulated miRNAs and mRNAs. In order to identify the commonly influenced mechanisms due to OTA treatments, the significantly enriched pathways of miRNAs and mRNAs analyses were compared. By comparing the mRNA and miRNA results, 21 and 13 commonly enriched pathways were identified in high-dose and medial-dose groups, respectively ( Seven pathways including amino acid metabolism, lipid metabolism, signaling molecules and interaction, and xenobiotics biodegradation and metabolism were commonly identified in high-dose and medialdose groups. In the high-dose group, the pathways influenced by OTA were solely associated with circulatory system, digestive system, endocrine system, excretory system, lipid metabolism, and signal transduction. Nevertheless, in the medial-dose group, the effects of OTA were relatively specific, primarily centering on carbohydrate metabolism, amino acid metabolism, metabolism of cofactors and vitamins, and signaling molecules and interaction. These pathway results also suggest that different toxicity mechanisms are provoked by different doses of OTA. Cross-omics analyses and the identification of five common pathways. Some miRNAs induce mRNAs degradation and/or repress mRNA translation without changing mRNA expressions. A comprehensive regulatory network was constructed by integrating hepatic miRNA, mRNA, and protein profiles in response to OTA treatments. In high-dose group, a proteomic approach was adopted by using two-dimensional gel electrophoresis and applying a threshold cutoff of at least 1.5 range ratio and p , 0.05. A total of 61 proteins were significantly different between the high-dose and control group (Fig. 3). A pathway enrichment analysis on these 61 proteins was performed and compared with those pathways obtained from miRNA and mRNAs profiling studies in the high-dose group. Five pathways encompassing cysteine and methionine metabolism, PPAR signaling, primary bile acid biosynthesis, arginine and proline metabolism, and metabolism of xenobiotics by cytochrome P450 were commonly identified among all the three profiles. Differentially regulated miRNAs among these five different pathways are described in Table 3. The mRNAs and their corresponding miRNAs related to these five most interesting pathways are shown in Table 4. Table 5 depicts the possible interactions among the identified proteins and the deregulated miRNAs in the five most relevant pathways.
Interestingly, these miRNAs were found to have a distinct expression pattern in medial-dose and high-dose groups compared to control group. The expression patterns of these miRNAs were found to be dose-related in PPAR signaling, arginine metabolism and proline metabolism pathways. However, in cysteine and methionine metabolism pathway, all of these miRNAs showed highly elevated expressions in medial-dose group compared to high-dose group. In primary bile acid biosynthesis pathway, some but not all of these miRNAs were noticeably dose-related. Although the proteomic analysis was not conducted for the medial-dose group, analyses of the transcriptomic data demonstrated that similar signaling pathways (cysteine and methionine metabolism, arginine and proline meta-  bolism and metabolism of xenobiotics by cytochrome P450) were enriched as compared to control group. This suggests that these three pathways are also modulated by medial-dose of OTA. In addition, the expression of three genes (Sds, Ass1 and Rgd1559459) were significantly different between medial-dose and control group. These genes are the representative members of cysteine and methionine metabolism, arginine and proline metabolism, and metabolism of xenobiotics by cytochrome P450 pathways.
Proteins involved in these five most interesting pathways were listed in Table 6. Carbamoyl-phosphate synthase (Cps1) and S-adenosylmethionine synthase isoform type-1 (Mat1a) were upregulated, while alpha-methylacyl-CoA racemase (Amacr), hydroxymethylglutaryl-CoA synthase (Hmgcs2), glutathione S-transferase Y-b subunit (Gstm2), ornithine carbamoyltransferase (Otc) and cytosol aminopeptidase (Lap3) were down-regulated. Surprisingly, only carbamoyl-phosphate synthase and its encoding gene both were up-regulated, whereas the remaining six genes were not differentially expressed in the high-dose group based on the transcriptomic data. These observations suggest a mode of action of miRNAs in which they only inhibit the translational process but not drive mRNA repression.
PCR validation of miRNAs and mRNAs in these pathways. The miRNAs and mRNAs involved in these five most relevant pathways were validated by qRT-PCR. The expression of rno-miR-195-5p, miR-130a-3p, miR-30a-5p, miR-140-5p, and miR-128-3p were examined. Compared to control livers, the expression of rno-miR-195-5p was significantly increased, while the expression of rno-miR-130a-3p was significantly decreased in high-dose group (p , 0.05). Despite of being statistically insignificant, rno-miR-30a-5p was observed with an up-regulation trend while two other miRNAs: rno-miR-140-5p and rno-miR-128-3p were identified with a down-regulation trend in the high-dose compared to control group (Fig. 4A). Additionally, qRT-PCR experiments were carried out to verify transcriptional results in which the expression patterns of Sds and Cps1 were significantly increased in the high-dose group (p , 0.05). Ass, Cyp8b1, and Cyp2c23 were exhibited an upward trend (p 5 0.083, 0.051, 0.549, respectively) (Fig. 4B). qRT-PCR results of miRNA and mRNA were mostly consistent with the sequencing findings, illustrating an agreement with the sequencing profiles.
Three important enzymes related with miRNA processing. The expression of key miRNA processing regulators (Drosha, Dicer1 and DGCR8) was determined. As shown in Fig. 5, mRNA levels of DGCR8 and Dicer1 were significantly reduced in high-dose group compared with control group. These results are consistent with previous reports on liver carcinogenesis and regeneration studies 24,25 . Also, our proteomic data showed that some of the differentially expressed proteins were closely related with liver regeneration as well (data not shown).
Impact of OTA on the liver. To elucidate the impact of OTA on hepatotoxicity, we further investigated the five most significantly altered pathways (Table 4) induced by OTA. Fig. 6 illustrated the regulatory network of significantly regulated genes, miRNAs and proteins by high-dose OTA in the livers. Three miRNAs: rno-miR-30a-5p, miR-30d-5p, and miR-30e-5p were involved in both cysteine and methionine metabolism, as well as in primary bile acid biosynthesis. rno-miR-130a-3p, rno-miR-130b-3p, and Cyp8b1 were associated with both PPAR signaling pathway and primary bile acid biosynthesis. rno-miR-122-5p was participated in both primary bile acid biosynthesis and arginine and proline metabolism, whereas the rno-miR-128-3p was involved in arginine and proline metabolism and metabolism of xenobiotics by cytochrome P450. Altogether, analyses of these pathways suggest an implication of OTA in hepatotoxicity by affecting liver functions or inducing metabolic diseases.

Discussion
To the best of our knowledge, this is the first study of its own kind in which an integrative approach encompassing miRNA, transcriptomic and proteomic profiling was adopted to explore OTA-induced early hepatotoxicity. We have employed and integrated three omics approaches (miRNA, mRNA and proteomic profiles) and identified five enriched pathways (primary bile acid biosynthesis, metabolism of xenobiotics by cytochrome P450, PPAR signaling pathway, cysteine and methionine metabolism, and arginine and proline metabolism) that are most likely to be unique and critical in the understanding of hepatic response to OTA. In the present study, rats were gavaged with OTA doses (0, 70 or 210 mg/kg body weight (bw)), and miRNA, mRNA and protein pro- filing studies were performed in the livers after thirteen weeks of gavage. The rat livers treated with OTA for 2, 4 and 13 weeks were comparable with controls (without OTA-treatment). However, bile duct hyperplasia appeared in rat liver treated with OTA for 26 weeks (data not shown). Moreover, lactate dehydrogenase (LDH) and aspartate transaminase (AST) at 13 weeks were lower in the highdose compared to control rats, while at 4 weeks, LDH and AST were not significantly different between high-dose group and the control rats. Both LDH and AST could reflect the liver damage to some extent and suggest that the liver damage could be more severe in high-dose group at 13 weeks compared to OTA-treated four-week animals. In addition, the differences of the serum biochemical indexes at 13 weeks were higher in high-dose (compared to control) livers than medium-dose (compared to control) livers. Medial-dose OTA does not have impact on the serum chemical or histopathological indexes, which is consistent with OTA studies 4-6 . Nonetheless, the derivatives found in the liver indicates that OTA must be metabolized in order to act as a carcinogen 26 . Therefore, we selected 13 weeks as a critical time point to find out the potential candidates and their signaling cascades which may be closely associated with early liver damage. No pathological changes were noticed at 13 weeks, but we wanted to elucidate the early changes in livers caused by OTA using multi-omics approach which could identify the early metabolic and toxicological impacts of OTA and provide an insightful resource for the prevention and/or attenuation of OTA toxicity in humans. Some of the OTA pathways resulting from our cross-omics approaches are in line with those reported previously [27][28][29] . Briefly, in a previous study, three pathways (PPAR signaling pathway, cysteine and methionine metabolism, and arginine and proline metabolism) have been found to be closely associated with metabolic syndrome, especially obesity 27 . In another study, the digestive system is shown to be impacted by OTA in the high-dose group 28 . Moreover, the carbohydrate metabolism pathway has been identified enriched in the medial-dose group 29 .
The results of our study provide a comprehensive landscape of OTA-induced hepatotoxicity that has never been revealed before. Furthermore, our sequencing data suggest an OTA hepatotoxicity based on OTA doses. In particular, the findings resulting from hierarchical clustering demonstrate that all miRNA profiles in control tissue are reminiscent of those in the high-dose group but not similar in the medial-dose group. Also, the enriched pathways obtained on miRNAs and mRNAs clearly indicate that, although both high-and medial-dose OTA markedly changed the gene expression, their influenced effects were quite different. The number of differentially expressed genes by high-dose OTA was marginally higher than those observed in medial-dose OTA. These observations suggest that the high-dose OTA exerts a more extensive effect on the rats, whereas, the medial-dose OTA acts in a more concentrated manner.
Mechanistically, OTA at the medial dose mainly influences energy metabolism such as carbohydrate metabolism, and basic substance metabolism (including amino acid metabolism, metabolism of cofactors and vitamins). However, in the high-dose group, the pathways influenced by OTA were solely associated with the different body systems such as circulatory system, digestive system, endocrine system, excretory system, lipid metabolism, and signal transduction. Leire et al. reported that the hepatic gene expression changes induced by OTA. The authors also observed that circadian exercise, cholesterol biosynthesis and lipid metabolism were influenced in the rat liver which is similar with our high dose group 30 . In another study with 300 ml OTA/kg bw, similar dose in our study, carbon metabolism, amino acid metabolism were influenced in the liver 31 . This comparative analysis suggests that a high-dose of OTA influences not only the basic metabolism but also endocrine environment, leading to the imbalance of whole body, while, a low-dose of OTA enables       the imbalance of energy metabolism, but the detoxification of the liver still functions. In contrast, a high-dose OTA treatment exhausts the liver repair system, influencing the overall defense reaction, for instance, circulatory system and excretory system broken down. Although, our sequencing data suggests that there are differential effects between high and medial dose OTA, this could not be confirmed by qRT-PCR. Further studies are needed to confirm our current findings. It should be noted that in the current manuscript, we only studied the toxicological mechanisms between the high-dose and medial-dose OTA at thirteen weeks. However, various toxic mechanisms could be expected at the different time intervals 30 . Moreover, several studies have been carried out to elucidate the impacts of different doses of other substances. For example, a recent study demonstrated that a high-dose cyclophosphamide is solely attributed to its direct cytotoxicity in contrast to the response to a low dose 32 . In another study, the expression of a few critical genes was altered by medial but not with the high doses of diethylnitrosamine, leading to non-linearity effects 33 .
In the present study, we identified dysregulation in critical pathways encompassing primary bile acid biosynthesis and metabolism of xenobiotics by cytochrome P450 in response to OTA, which provide a pathophysiological explanation to OTA-induced live damage. Bile acid is a primary component of bile and plays an important role in the absorption of not only fat and fat-soluble vitamins, but also xenobiotics. OTA is absorbed in the upper portion of the gastrointestinal tract and then transported to the liver. Thereafter, OTA is secreted along with bile acid from the liver to gastrointestinal tract through enterohepatic circulation. An increase in bile acid likely indicates increasing OTA absorption, resulting in the elevation of OTA-induced toxicity. Also, aflatoxins-induced acute liver injuries have been shown to be promoted by elevated serum bile acid levels. The metabolism of xenobiotics by cytochrome P450 is another most relevant pathway for proper liver functioning. Cytochrome P450 system protects liver by metabolizing exogenous substances. We also found a decrease in the expression level of glutathione S-transferase (Gstm2). This gene has previously been reported to be significantly down-regulated during the development of hepatocellular preneoplastic foci in rats 34 . In our study, the down-regulation of Gstm2 was only observed at the protein level, but not during mRNA profiling study. The one of the possible reasons behind this observation could be that miRNAs can regulate the gene expression at translational level instead of decaying mRNAs 35 . Another possibility could be that the post-transcriptional regulation, such as phosphorylation, ubiquitination, glycosylation, and acetylation, can make changes on proteins but not on mRNAs 36,37 .
Deregulated PPAR signaling pathway, cysteine and methionine metabolism, and arginine and proline metabolism are known to cause metabolic syndrome including obesity, dyslipidemia, hypertension, and elevated plasma glucose level. A review of OTA on renal cells has indicated that glomerular capillaries are damaged by hypertension 38 , which is thought to contribute to OTA-induced disorders. PPARs are closely related to energy status and metabolism. An aberration of PPARs leads to abnormal expression of genes in metabolic pathways, thereby contributing to etiology of metabolic syndrome 39 . In addition, PPARs are the major regulators of lipid and fatty acid metabolism. This may explain the metabolic shift of lipid pathways by high dose OTA, as well as changes of protein expression in lipid metabolism based on our proteomic data (data not shown). Moreover, coronary artery disease is associated with abnormal methionine metabolism. An elevated level of total homocysteine in the blood is a strong risk factor for atherosclerotic vascular disease in the coronary, cerebral, and peripheral vessels, and for arterial and venous thromboembolism. Altered methionine metabolism also plays a significant role in the protection against oxidative stress, suggesting that OTA could cause oxidative damage. Further future experiments are required to support the notion whether OTA treatment promotes cardiovascular disease -as arginine and proline metabolism signaling was highly significantly enriched in our study. Interestingly, arginine metabolism has previously been shown to play a major role in cardiovascular physiology and to ameliorate the metabolic syndrome in Zucker diabetic fatty Rats 40 . Furthermore, OTA has been reported to stimulate ceramide formation which participates in the pathophysiology of cardiovascular diseases 27 . Additionally, glutaminase 2 (Gls2) expression is induced in response to DNA damage or oxidative stress 41 , suggesting that OTA may also induces DNA damage or oxidative stress.
Notably, most of the pathways identified in our investigation have previously been explored in OTA toxicity related studies. In addition, Table 6 | Identification of liver proteins related with cysteine and methionine metabolism, PPAR signaling pathway, primary bile acid biosynthesis, arginine and proline metabolism and metabolism of xenobiotics by cytochrome P450 between cH and cC group using MS/MS analysis. Exopeptidase which selectively removes arginine and/or lysine residues from the N-terminus of several peptide substrates.
Spot numbers correspond to those in Fig. 6; Range ratio is the average change in abundance from three independent treatments (cH vs cC); Protein name, matched protein description; Protein ID, accession number from the NCBI database of matched proteins; Function, the funtion of the protein. the role of biotransformation in OTA toxicity has already been addressed in previous studies. Some of these investigators reported a down-regulation in the expression of the representative genes of transporters and xenobiotic metabolism 30,42 , while the remaining reports found an elevation in their expression 43 . Furthermore, energy metabolism was deregulated induced by OTA in vitro and in vivo studies. Also, many genes (for example, Glutathione S-transferase) involved in energy and xenobiotic metabolism of HepG2 liver cell exposed to OTA were down-regulated 7 . Our results also showed a reduction in this gene. Male Tsc2 EKER rats gavaged with OTA also changed their energy provision 43 . In addition, liver lipid metobolism was down-regulated of male F344 rats after OTA treatments in several studies 30,31 . In our result, the protein (such as Hmgcs2) was down-regulated in PPAR pathway, which is the major regulator of lipid and fatty acid metabolism. Also, the amino acids metabolism (such as alanine and aspartame metabolism, tryptophane metabolism) was affected in rats' liver 30 . We compared the miRNA pathways with at least two and three algorithms, and found that the lost pathways with at least three algorithms mainly concentrated on amino acid metabolism, carbohydrate metabolism, glycan biosynthesis, lipid metabolism, and metabolism of cofactors and vitamins. In previous studies on OTA, amino acid metabolism, carbon metabolism pathways were influenced in the liver 31 . In addition, genes related with glycan biosynthesis in HepG2 cells are changed in response to OTA treatment 7 . Also, the lipid metabolism was influenced in the rat liver induced by OTA 30 . Moreover, vitamins are important cofactors of many enzymes and thereby influence metabolic activities, playing a part in the effect of mycotoxins 44 . We therefore think that pathways  www.nature.com/scientificreports resulting from at least two algorithms are more comprehensive. Among the five interested pathways (primary bile acid biosynthesis, and metabolism of xenobiotics by cytochrome P450, arginine and proline metabolism, cysteine and methionine metabolism, and PPAR signaling pathways) in the current study, metabolism of xenobiotics by cytochrome P450, and arginine and proline metabolism were not identified with at least three algorithms. Cytochrome P450s are the critical biotransformation enzymes involved in OTA toxicity 45 . These pathways are lost either due to this highly stringent criterion (at least 3 algorithms) to reduce the number of false positive targets or could be that these pathways may be slightly modulated by miRNAs, but previous studies suggest their potential relationship with OTA. Pathways resulting with at least two algorithms contain more information; hence, the method using at least 2 algorithms is reasonable.
The proteomic and transcriptomic expression results of carbamoyl-phosphate synthase out of the five most interesting pathways were comparable; however, unfortunately, the correlation between mRNA and protein results was relatively low. Recently, the relationship between miRNA, mRNA and proteins has been explored by many researchers. In principle, a perfect or near-perfect Watson-Crick complementarity between the miRNA and its target mRNA results in mRNA cleavage and decay, whereas, an imperfect complementarity directs the translational repression. Such imperfect interactions of miRNAs regulate the gene expression at translational level instead of decaying mRNAs. Most animal miRNAs are imprecisely complementary to their mRNA targets and they inhibit protein synthesis through an unknown mechanism that preserves the stability of the mRNA target. Apart from the repression, several groups have recently reported that miRNAs can also activate rather than repress their mRNA targets under certain conditions, which is in contrast to the known ''dogma'' of miRNA regulation as was known previously. It is probably the reason that rno-miR-195-5p induction by OTA treatment along with its mRNA target (Sds) induction as well (Fig. 4). Therefore, it is ideal and necessary to incorporate results from systematic analyses of miRNA profiling, transcriptomics and proteomics for making a comprehensive and robust conclusion.
The combination of omics technologies allows an improved detailed analysis and identification of potential mechanisms to advance the overall understanding of OTA hepatotoxicity. Our crossomics results demonstrate that different doses of OTA exert a diverse pathophysiological impact and elucidate OTA-induced five most relevant pathways (cysteine and methionine metabolism, PPAR signaling, primary bile acid biosynthesis, arginine and proline metabolism, and metabolism of xenobiotics by cytochrome P450) are induced by OTA. Nonetheless, future investigations are required to study the mechanisms of OTA hepatotoxicity by further exploring these five pathways.

Methods
Ethics statement. Every possible effort was made to minimize animal suffering. All procedures were approved by the the Ethics Committee of China Agricultural www.nature.com/scientificreports University (permission number: 120020). All experiments were performed in accordance with relevant guidelines and regulations.
Animals. Specific pathogen-free male F344 rats were obtained at the age of six weeks and were kept in accordance with institutional guidelines [The Supervision and Testing Center for GMOs food safety, Ministry of Agriculture (Beijing, China) with the license number SYXK (Beijing) 2010-0036]. Rats were housed three per cage in a temperature controlled room (22 6 2uC) with a relative humidity of 40 , 70% and a 12 h light-dark cycle. Feed and sterilized water were consumed ad libitum. Rats (six per group) were randomly assigned to control group (C), medial-dose group (M), high-dose group (H), gavaged with OTA at 0, 70 or 210 mg/kg bw, for 2, 4, 13 and 26 weeks (designated as a, b, c and d) respectively. OTA was dissolved in corn oil (Aladin, Shanghai, China). Body weights were recorded weekly. After the gavage, rats were anesthetized using chloral hydrate (6%, 5 ml/kg bw, ip). Blood samples were collected from the inner canthus before decapitation for determination of blood parameters including alanine (ALT), aspartate transaminase (AST), alkaline phosphatase (ALP), albumin (ALB), creatinine (CREA), glucose (GLU), blood urea nitrogen (BUN), lactate dehydrogenase (LDH), high-density lipoprotein (HDL) and low-density lipoprotein (LDL) using a Hitachi 7020 automatic biochemical analyzer (Hitachi, Tokyo, Japan) 13 .
Preparation of livers. The liver was weighed, and a complete liver necropsy was performed. The livers were fixed in 4% buffered formaldehyde for at least 24 h before histological processing, embedded in paraffin, sectioned to 4-6 mm thick, and then stained with hematoxylin-eosin (H&E). Histopathological examination of tissue sections was conducted by a pathologist at the Experimental Animal Research Center, China Agricultural University. Part of livers were frozen immediately in liquid nitrogen and kept at 280uC for further studies.
miRNA profiling and statistical analysis. Small RNA-containing total RNA was extracted from the thirteen-week livers using mirVana TM miRNA isolation kit (Ambion, Austin, TX, USA). Six samples from one treatment were pooled together to generate one group. RNA quality was assessed and verified by using a BioAnalyzer 2100 (Agilent Technology, Santa Clara, USA). For each sample, small RNA enriched by PEG8000 precipitation from 10 mg total RNA was used for library constructions following the Illumina' protocol. The small RNAs ligated with 59 and 39 adaptors were reverse transcribed and amplified. The libraries were quantified by ECO (Illumina, San Diego, USA) and sequenced by using the Solexa's proprietary sequencing-bysynthesis method. The resulting images were analyzed to generate raw digital-quality data. After masking adaptor sequences and removal of contaminated reads, clean reads were selected for further analysis. These clean reads were aligned with the Rfam database (ftp://selab.janelia.org/pub/Rfam) to match with known rRNA, snRNA, snoRNA and tRNA sequences. After non-coding RNAs were removed, the remaining clean reads were further compared with the pre-miRNAs database in miRBase. The matched reads were used to identify mature miRNAs, and their reads were counted. Chi-square test was applied to calculate the difference between groups with 5% level of significance. The numbers of miRNA reads were normalized by tags per million (TPM) values (TPM 5 (miRNA total reads/total clean reads) 3 10 6 ). Different genes were set by the threshold of false discovery rate (FDR) , 0.05. Hierarchical clustering based on average linkage was performed using J-express. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis on differentially expressed miRNAs was performed using R script by implementing Fishers' exact test with the Benjamini and Hochberg (BH) as multiple testing methods with 5% level of significance. mRNA sequencing and statistical Analysis. Total RNA in the livers of rats gavaged for thirteen weeks were extracted by using mirVana TM miRNA isolation kit (Ambion, Austin, TX, USA). Two rats were randomly choosen as biological replicates in one treatment. RNA quality was assessed by using a BioAnalyzer 2100 (Agilent Technology, Santa Clara, USA). Enriched mRNAs were used for library constructions following Illumina's protocols. The mRNAs were reversely transcribed, end repaired, and then the adaptor addition PCR products were enriched. The library was assessed using a BioAnalyzer 2100 (Agilent Technology, Santa Clara, USA) and the concentration reached 2 nM. mRNA was sequenced on the Illumina HiSeq TM 2000 (Illumina, San Diego,USA). Clean reads were acquired after masking adaptor sequences and removing reads of low quality. Gene expression level was estimated by using Reads Per Kilo bases per Million reads (RPKM) method. Different genes were detected based on DESeq, controlling FDR , 0.05. KEGG analysis on differentially expressed genes was performed by using KOBAS software.
miRNA target predictions and enrichment analysis. An ''in-silico'' screening of miRNA binding site prediction was carried out within deregulated genes by adopting the miRWalk database 15 , including miRanda, miRDB, miRWalk, PITA, TargetScan and RNAhybrid. The targets predicted by at least two different programs were chosen for further analysis to effectively reduce false positive targets 13,[17][18][19][20][21][22] . An in-silico analysis was used for identifying possible interactions between miRNAs and mRNAs. The predicted target genes were subjected to KEGG enrichment analysis with Benjamini and Hochberg (BH) as multiple testing method (background correction) with 5% level of significance.
Relative quantitative protein profiling and mass spectrometry. The livers of rats gavaged for thirteen weeks at th cC and cH were homogenized and re-suspended in a lysis buffer containing 7 M urea, 2 M thiourea, 4% (w/v) CHAPS, 1% (w/v) DL-Dithiothreitol (DTT), 1% (v/v) pH 4-7 IPG buffer (GE Healthcare, NJ, USA) and 0.5% (v/v) protease inhibitor cocktail, sonicated intermittently for 90 s and then centrifuged at 15000 g for 20 min at 4uC. The supernatants were mixed (six rats per group) and carried out in triplicate. Protein concentrations were determined by the Coomassie blue method 46 . For the first dimension, samples containing 550 mg total proteins were mixed in 250 mL of rehydration buffer containing 7 M urea, 2 M thiourea, 2% (w/v) CHAPS, 65 mM DTT, 2% (v/v) pH 4-7 IPG buffer (GE Healthcare, NJ, USA), and 0.002% bromophenol blue. Isoelectric focusing (IEF) was performed in IPG strips (pH 4-7, 13 cm, GE Healthcare) on a Multiphor IIIsystem (GE Healthcare, NJ, USA) using the following program: 500 V for 500 Vh, 1000 V for 800 Vh, 8000 V for 11500 Vh, 8000 V for 7500 Vh, 500 V for 10 h. The second dimension electrophoresis was performed in 12.5% SDS polyacrylamide gels at 2 W/ gel for 45 min, followed by separation at 4 W/gel until the dye front reached the bottom of the gel. The proteins were visualized with Coomassie Brilliant Blue R-250 after a 1 h protein fixation in 50% ethanol, 10% acetic acid and 40% water. Gels were destained with a solution of 30% ethanol, 8% acetic acid and 62% water for 2 h. The resulting digital images were analyzed and quantified using the Image-Master 2D 7.0 software (GE Healthcare, NJ, USA). Three independent experiments were conducted. One-way ANOVA (p , 0.05) was used for selecting the different spots between the two groups (cH vs cC).
For mass spectrometry (MS), spot picking of interest was carried out with preparative gels and subjected to in-gel trypsin digestion according to Sun et al. 47 with minor modifications. Matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) and tandem TOF/TOF MS analyses were carried out on a 4800 Proteomics Analyzer (Applied Biosystems, Foster City, CA, USA) according to Sun et al. 48 . GPS Explorer TM software, version 3.6 (Applied Biosystems, Foster City, CA, USA) was used for creating and searching files with the MASCOT (Matrix Science, http://www.matrixscience.com/) for peptide and protein identification. The highest protein score (top rank) was singled out from the multiprotein family. Precursor error tolerance was set to 60.2 Da, and MS/MS fragment error tolerance was set to 60.3 Da. All of the identified proteins had a protein score .66 with an expected pvalue less than 0.05.
Reverse transcription and quantitative real-time PCR (qRT-PCR) analysis of miRNAs and mRNA. Isolated miRNAs and mRNAs were reverse transcribed by using miRcute microRNA first-strand cDNAsynthesis kit (Tiangen, Beijing, China) and TIANscript RT Kit (Tiangen, Beijing, China), respectively. Quantitative real-time PCR (qRT-PCR) analyses on the miRNAs and the mRNAs were performed by using the RealMasterMix (SYBR green I) (Tiangen, Beijing, China) in a Bio-Rad CFX96 Real-time PCR System (Bio-Rad, Hercules, CA, USA). The thermal cycling program was set as follows: 95uC for 5 min, followed by 40 cycles of 95uC for 30 s, 60uC for 30 s, and 72uC for 30 s. Following the amplification, melting curve analysis was conducted by heating from 65uC to 95uC with increments of 0.5uC for 5 s to determine the specificity of the PCR reactions. The internal controls for miRNA 49 and mRNA 50 normalization were 5 s RNA and b-actin. The primer sequences were shown in Table S1. Relative gene expression was evaluated using the 2 2DDCT method. Six biological replicates were used, and each were performed in triplicates.
Statistics. The false discovery rate (FDR) corrected p-value of 0.05 was considered as statistically significant. KEGG analysis on differentially expressed miRNAs was performed using a customized ''R script'' by implementing Fishers' exact test with the Benjamini and Hochberg (BH) as multiple testing method with 5% level of significance. For RT-PCR experiments, results were statistically evaluated by using ANOVA with Duncan test and were given as mean 6 standard deviation (SD). A pvalue of ,0.05 was considered statistically significant.