Transcriptome Analyses Revealed Adaptive-Responses in Livers of Mice Treated With Hua-Feng-Dan and Its “Guide Drug” Yaomu

Hua-Feng-Dan is a patent Chinese medicine for stroke recovery and is effective against Parkinson’s disease models with modulatory effects on gut microbiota, but its effects on hepatic gene expression are unknown. This study used RNA-Seq to prole hepatic gene expression by Hua-Feng-Dan and its “Guide Drug” Yaomu. Mice received orally Hua-Feng-Dan 1.2 g/kg, Yaomu 0.1-0.3 g/kg, or vehicle for 7 days. Liver pathology was examined, and total RNA was isolated for RNA-Seq. The bioinformatics, including GO and KEGG pathway enrichment analysis, two-dimensional clustering, Ingenuity Pathways Analysis (IPA), and Illumina BaseSpace Correlation Engine were used to analyze differentially expressed genes (DEGs). qPCR was performed to verify selected genes.


Results
Hua-Feng-Dan and Yaomu did not produce liver toxicity as evidenced by histopathology and serum ALT and AST. GO Enrichment revealed Hua-Feng-Dan affected lipid homeostasis, protein folding and cell adhesion. KEGG showed activated cholesterol metabolism, bile secretion and PPAR signaling pathways.

Conclusion
Hua-Feng-Dan at clinical dose did not produce liver pathological changes but induced metabolic and signaling pathway activations. Low dose of its Guide Drug Yaomu produced similar changes to a lesser extent, but high dose of Yaomu had little effects. The effects of Hua-Feng-Dan on liver transcriptome changes may produce adaptive responses to program the liver to produce bene cial or detrimental (overdosed) pharmacological effects.

Background
Hua-Feng-Dan has over 300 years of history in the treatment of stroke and various diseases in China, and is still used today alone or in combination with other medications [1]. Hua-Feng-Dan contains cinnabar (HgS) and realgar (As 4 S 4 ), and our initial research focused on safety evaluation of Hua-Feng-Dan following acute [2], subacute [3,4], and subchronic [5,6] administrations, and demonstrated that the toxicity of Hua-Feng-Dan is different from environmental mercury (Hg) and arsenic (As) compounds [1].
Yaomu is considered as the "Guide Drug" in Hua-Feng-Dan. "Guide drug" is an important theory in Chinese medicines as exempli ed by Licorice (also called Gan-Cao, Glycyrrhiza uralensis Fisch). Licorice as a "Guide Drug" is essential to enhance e cacy and reduce toxicity when combined with other herbal medicine preparations [10,11]. Yaomu in Hua-Feng-Dan consists of 7 components (Typhonium giganteum Engl; Pinellia ternata (Thunb.) Breit; Arisaema heterophyllum Blume; Aconitum carmichaeli Debx; Curcuma wenyujin Y. H. Chen et C.Ling; Medicated Leaven (Shen-Qu); and cow bile water), and is subjected to fermentation for three months. Fermentation is an important processing technology in traditional medicines under appropriate temperature, humidity, and microorganisms to reduce toxicity, enhance e cacy, and generate active inter-metabolites [12,13]. Limited studies on Yaomu in Hua-Feng-Dan are found in regard to processing technology [14], and fermented Yaomu could reduce the content and toxicity of aconitine, mesaconitine, and hypaconitine [15]. However, little is known about hepatic effects of Yaomu.
The liver is the main organ of the phase I and phase II drug metabolism and plays a leading role in the elimination of oral drugs [16]. The liver is also an important organ of detoxi cation in response to xenobiotics. To "program the liver" is an important concept in that xenobiotics at appropriate doses could evoke adaptive responses to protect against environmental and external damage, and to induce a number of molecular events to exert bene cial effects against various disorders [16]. The use of traditional medicines to modulate innate and adaptive immune functions has been reviewed [17]. Oleanolic acid, a triterpenoid from many herbs/fruits, could reprogram the liver to protect against a variety of hepatotoxicants at the low doses, but is hepatotoxic at higher doses [18]. Thus, learn to program the liver would help us to understand the bene cial or harmful effects of traditional medicines such as Hua-Feng-Dan.
RNA-seq is an advanced technology for high-throughput analysis of mRNA, small RNA, and non-coding RNA, etc. RNA-Seq is also superior to microarray for genomic pro ling [19], and has been recommended to study traditional Chinese medicines [20]. Therefore, this study utilized RNA-Seq to pro le the liver transcriptome after Hua-Feng-Dan and Yaomu treatment. After quality evaluation of RNA-Seq analysis data, differentially expressed genes (DEGs) were identi ed with the DESeq2 method compared to controls. A comprehensive bioinformatics including two-dimensional clustering, GO and KEGG pathway enrichment analysis, ingenuity pathways analysis (IPA), and Illumina BaseSpace Correlation Engine (BSCE) were used to analyze DEGs and to mine adaptive response changes.

Material And Methods
Hua-Feng-Dan and its "Guide Drug" Yaomu Hua-Feng-Dan (HFD) and its Guide Drug Yaomu (YM) were provided by Hua-Feng-Dan Pharmaceutical Co. (Guizhou, China). Yaomu is the most important part of the traditional medicine Hua-Feng-Dan [15], and fermented for three months before using as the "Guide Drug" to make Hua-Feng-Dan with other herbs, animal-based products, and minerals as listed in Introduction and reported previously [1,[5][6][7][8].
Animal and drug treatment Adult male C57BL/6J mice (20-22 g) were obtained from the SPF (Beijing) Biotechnology Co.,Ltd After acclimatization for one week, mice were randomly divided into four groups: Control group (n = 5),YM-0.1 (n = 5, Yaomu 0.1 g/kg), YM-0.3 (n = 7, Yaomu 0.3 g/kg), HFD (n = 5, Hua-Feng-Dan 1.2 g/kg). Mice received oral gavage daily for 7 consecutive days; the dose of Hua-Feng-Dan selection was based on recent publications [5,6], and Yaomu was estimated from the recipe. Twenty-four hours after the last administration, mice were anesthetized, blood and livers were collected. Histopathology A portion of livers from the same site was cut and xed with 10% formaldehyde, dehydrate with ethanol gradients, and embed in para n. Liver tissue was cut into 4 µm slices with a tissue slicer (model: RM2245), sections were depara nized with xylene, hydrated with gradient ethanol, stained with conventional HE, and observed under upright optical microscope (Olympus, Tokyo, Japan).

Biochemical Analysis
After the mice were anesthetized, the blood was collected and left standing at room temperature for 2 hours and centrifuged at 3500 rpm for 10 min to collect the serum. The ALT and AST test kits were used to detect the activity of ALT and AST in serum (Nanjing Jiancheng Institute of Bioengineering, Nanjing, China).

RNA extraction
Total RNA was extracted from mouse liver tissues using RNAiso plus (TAKARA, Dalian, China). The concentration of total RNA is measured with ND-2000 nanodrop, with A260 /A280 > 1.8.

RNA-Seq
The EASY RNA-Seq method was used to complete the construction of the sample library, and after quality control testing, the Illumina platform was used for sequencing by Chongqing Weilang Biotechnology Co., Ltd. (Chongqing, China). Brie y, RNA is subjected to reverse transcription reaction with oligo(dT) primer to generate the rst-strand cDNA; through the joint reaction of RNase H enzyme, DNA polymerase and T4 ligase, double-strand cDNA was generated. The Agencourt AMPure XP Beads were used to purify cDNA by magnetic separation. After the double-stranded cDNA is fragmented by Tn5 enzyme, the RD sequence required for sequencing is added to both ends, and then P5 and P7 are connected to the RD sequence according to the TruePrep DNA Library Prep Kit with P5/P7 adapter primers (Vazyme, cat. TD503). Sequencing primers were at both ends, and PCR was ampli ed by VAHTS™ DNA clean beads kit (Vazyme, cat. N411-02). The Illumina HiSeq platform with a 150 bp double-ended model was used for sequencing analysis.
Bioinformatics RNA-Seq generated fragments per kilobase of exon per million fragments mapped (FPKM) were transformed to raw gene counts for quality control analysis with all parameters performed by the sequencing company (Chongqing Weilang Biotechnology Co., Ltd., China). The initial bioinformatics analysis was also performed by including two group comparisons using DESeq2 methods, GO, and KEGG analysis with company service.
Gene ontology (GO) and KEGG pathway enrichment analysis. The differentially expressed genes between the two groups were analyzed by Gene Ontology (GO) enrichment, and the hypergeometric distribution was used to test the signi cantly enriched GO entries. ClusterPro ler software [21] was used to analyze the Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway enrichment of differentially expressed genes between groups. Use bubble charts to visualize the top functional pathways in the analysis results of GO and KEGG enrichment (p-value < 0.05). The size of the dot in the gure represents the number of differential genes that can be annotated in the functional pathway, and the shade of the color represents the signi cant degree of enrichment of the functional pathway. Rich factor is a score used to evaluate the degree of enrichment of the functional pathway.
Ingenuity Pathway Analysis. Ingenuity Pathway Analysis (IPA) server (Qiagen, Redwood City, CA) was used to perform upstream regulators analysis. DEGs (p < 0.05) from HFD, YM-0.1, YM-0.3 treatments vs control, respectively were subjected to IPA Core Analysis, followed by Comparative Analysis to compare the Z-score of the changes among treatment groups.
BaseSpace Correlation Engine analysis. BaseSpace Correlation Engine (BSCE, Illumina, CA) is the RNA sequencing and microarray database curated over 23,000 scienti c studies to get data-driven answers for genes, experiments, drugs and phenotypes for the research. DEGs from various treatment groups vs Controls were uploaded to BSCE and compared with all biosets in the database using the Running Fisher test. This method provides an assessment of the statistical signi cance of the correlation of the overlapping genes between each treatment DEGs and biosets curated in BSCE, with a summary p-value.
The results were exported, and each p-value was converted to a -log(p-value). Biosets with -log(p-value) + 4 or -4 were considered signi cant correlation [22]. A column in the gene expression spreadsheet was populated with the -log(p-value) for each bioset. Biosets that were positively correlated with the DEGs were predicted to produce similar effects, either directly or indirectly; biosets that were negatively correlated were predicted to suppress the biological responses.
Real-time qPCR. The total RNA was reverse transcribed with 5 x PrimeScripe Buffer, PrimeScript RT Enzyme MIX, Oligo dT Primer and Random 6 mers. The iTaq™ Universal SYBR green PCR Super Mix (BIO-RAD) was used for real time qPCR analysis. The reaction parameters were 95 ℃ 3 min, 95 ℃ 10 s, 60 ℃ 45 s, 40 cycles in total. The Ct values were normalized to β-actin of the same sample and calculated by the 2 −△△Ct method and expressed as percentage of Control. The primers were designed by Primer3 (version 4), and synthesized by Shanghai Sangon Biotech as described [6].

Statistical Analysis
The DESeq2 method was used to compare differences between two groups, and the criteria were set at p < 0.05. For real time qPCR analysis, the data presented were mean ± SEM. The Normality Test (Shapiro-Wilk) was rst performed. Kruskal-Wallis one-way AONVA was performed, followed by Dunn's multiple range test (SigmaPlot v14), the signi cance was set at p < 0.05.

Results
The effect of Hua-Feng-Dan and Yaomu on serum enzyme activities and liver pathology Seven days after administration, the body weight of mice in each group increased slightly (The initial weights were 22  The liver lobules of the mice in the Control group and treatment groups were clear and complete, without degeneration and necrosis of hepatocytes, and the liver cords were arranged radially around the central vein. No obvious pathological damage was observed in the liver after administration of different doses of Yaomu (YM) and Hua-Feng-Dan (HFD) to mice (Fig. 1).

Gene Expression Pattern overview
The results (~22500/sample) obtained from RNA-Seq FPKM were transformed to gene raw counts, and then subjected to cluster analysis. The clustering heatmaps of gene expressions from 12 individual samples ( Fig. 2A) and 4 groups (Fig. 2B) are shown. Red represents the up-regulation; blue represents the down-regulation. The color brightness is associated with differences. Fig. 2A shows

GO and KEGG pathway enrichment analysis
The GO database standardizes the description of differentially expressed genes in terms of function, participating biological pathways, and cell location. The KEGG is a database of the metabolic pathways of gene products in cells and the functions of these gene products. Here we focus on the comparison between Hua-Feng-Dan and Control group. Figure 3 summarizes the top 20 enriched by GO and the top 7 enriched by KEGG pathway. GO Enrichment Analysis shows that the biological process of differentially expressed genes was mainly involved in lipid metabolism, sterol homeostasis, cholesterol homeostasis, intestinal absorption, protein folding and cell adhesion ( Figure 3A). KEGG Enrichment Analysis shows that pathways of differentially expressed genes were mainly involved in cholesterol metabolism, bile secretion, PPAR signaling pathway, drug metabolism, fat digestion and absorption, retinol metabolism and the like ( Figure 3B).

Differentially Expressed Gene Analysis
Differentially expressed genes (DEGs) were analyzed via the DESeq2 method compared to Controls.
DEGs were set at p < 0.05. Using complete clustering, the DEGs were clustered and input into Treeview to visualize the differences between groups (red indicates up-regulation and blue indicates downregulation). From Figure 4, compared to Hua-Feng-Dan group (806 DEGs), YM-L group had 235 DEGs and YM-H group had 92 DEGs. Four clusters were selected for annotation (two upregulated: line 42-150, 109 genes increased in HFD and YM-L group only, involved in cellular function and signal regulation; line 151-179, 28 genes involved in cellular function, circadian, and signal transduction were increased in all three groups. Two downregulated cluster: line 494-556, 63 genes decreased in HFD and YM-L group only, while line 557-567, 11 genes decreased in all three groups. These genes are involved in Phase-I, II, and III metabolism and immunomodulation. The clustered DEGs with the annotation for gene names are shown in Supplementary Table 1. qPCR analysis of selected DEGs Selected 10 DEGs were further analyzed via real time qPCR ( Figure 5). The expression of Cyp2a4 was increased 2.2-fold by Hua-Feng-Dan (HFD), but was unaffected by YM-0.1 and YM-0.3; the expression of Cyp4a14 was increased 2.3-fold by HFD, YM-0.1 and 0.3 did not affect its expression. The expression of Fgf21 was increased by HFD (3.7-fold), although not reaching signi cance; the expression of Gm3776 was increased 6.5-fold by HFD; The expression of Mt2 (10-fold) and Gsta1 (4.7-fold) was increased by HFD. The expression of Il1rn (24-fold) and Egr1 (8.9-fold) was increased by HFD. The expression of Trib3 (1.7-fold) and Bcl2a1d (1.7-fold) was tended to increase, but not signi cant. The two doses of YM had little effects on the expression of these genes.

Ingenuity Pathways Analysis of DEGs
The IPA upstream regulator analysis ( Figure 6) was used to identify the upstream regulators that may be responsible for gene expression changes observed in the study. IPA predicts which upstream regulators are activated or inhibited to explain the up-regulated and down-regulated genes observed. The top 25 upstream regulators were selected, including 18 increased by Hua-Feng-Dan with molecules (Tgf beta; NRG1; Vegf; NFκB (complex); JUN; STAT3; CREB1; P38 MAPK; CHUK; EGR1; ERK1/2; and EP300) or with chemicals (Tretinoin; carbon tetrachloride; decitabine; ursodeoxycholic acid; AGN194204 and 4hydroxytamoxifen). The 7 downregulated regulators were molecules (Let-7a-5p, ACOX1, and DICER1) and chemicals (MEK1/2 inhibitor U0126; p38 MAPK inhibitor SB203580; JNK inhibitor SP600125; and PIK3 inhibitor LY294002). All of these upstream regulators point towards activation of MAPK signaling pathways and adaptive responses with the downregulation of inhibitors for MAPK and PIK3 signaling pathways. YM-0.1 produced similar upstream regulator alternations in the similar direction, except for ursodeoxycholic acid; while the high dose of YM (0.3g/kg) produced weak alterations in these regulators, with 3 in opposite directions (different color).

Illumina BaseSpace Correlation Engine analysis of DEGs
All the biosets were imported into BSCE for curated studies, ltered by mice, RNA expression and treatment vs control, and the curated les were exported and the -log (p-values) were calculated and the VLOOKUP function was used for correlation comparison with the GEO database. There were very high correlations between HFD-induced DEGs and 25 GSE databases with -log(p-values) from 12 to 26.9 (Figure 7). The deeper the red color, the higher -log (p-values), and the values >4 is considered signi cant [22]. YM-0.1 moderately correlated with 18 GSE databases in the same direction but to the lesser extent, while YM-0.3 only weakly correlated with 9 databases, and negatively correlated with 3 databases (blue color). Thus, Hua-Feng-Dan at the clinical dose mimics chemical-induced adaptive responses, while YM-0.1 could mimic some effects of HFD, but YM-0.3 has little effect, suggesting that the use of Yaomu at appropriate doses is important to induce adaptive responses.

Discussion
This study demonstrated that Hua-Feng-Dan and its "Guide Drug" Yaomu at clinical doses did not produce damage to the liver but produced gene expression changes. RNA-Seq revealed 806 DEGs in Hua-Feng-Dan-treated livers compared to controls, and low-dose of Yaomu produced 235 DEGs in the same direction. GO and KEGG analysis revealed Hua-Feng-Dan affected metabolism and singling pathways. IPA upstream analysis pointed towards activation of MAPK signaling pathways and adaptive responses by Hua-Feng-Dan. qPCR veri ed 10 selected DEGs. Hua-Feng-Dan induced gene expression changes were highly correlated with the GEO database of chemical-treated livers of mice, and Yaomu-0.1 had similar effects, but to a lesser extent. This RNA-Seq study clearly demonstrated Hua-Feng-Dan and low dose of Yaomu induced adaptive responses to reprogram the liver for biological effects as seen in Chinese medicine formulas, or in Chinese medicine extracts such as oleanolic acid [16][17][18].
We have shown that cinnabar (HgS)-containing Hua-Feng-Dan is different from environmental mercury compounds (HgCl 2 , MeHg) in producing liver injury [2,3]. In the present study, Hua-Feng-Dan and Yaomu at clinical doses also did not produce liver injury as evidenced by serum enzyme activities and histopathology, consistent with prior observations. This is an important pre-requirement for RNA-Seq experiments and results interpretation to help us understand how Hua-Feng-Dan could program the liver to produce adaptive responses as seen with oleanolic acid, an active ingredient of many Chinese medicines [16,18].
GO and KEGG pathway enrichment analysis showed the top 20 of GO enrichment and the top 7 of KEGG pathway enrichment. For example, high-fat, high-energy diets can increase the incidence of lipid metabolism disorders [23], and lipid transport disorders can cause lipid accumulation in liver cells, leading to fatty liver [24]. The effects of Hua-Feng-Dan on lipid metabolism could produce adaptive or detrimental effects. For the KEGG pathways, cholesterol metabolism was at the top. Cholesterol is essential for maintaining membrane uidity and is regulated by nuclear receptors [25] and cholesterol homeostasis is important for biological functions and bile acid synthesis [26,27]. The effects of Hua-Feng-Dan on cholesterol metabolism and bile acid homeostasis could be important adaptive responses. Peroxisome proliferation-activated receptors (PPARs) regulate the lipid metabolism and cellular immune and biological functions [28], and effects of Hua-Feng-Dan on PPAR signaling pathway and drug metabolism genes are important to program the liver. Retinoid metabolism and signaling are necessary to maintain immunity, cell proliferation and differentiation [29], which was affected by Hua-Feng-Dan. All these changes imply the ability of Hua-Feng-Dan to program the liver.
Two-dimensional clustering analysis of DEGs showed that Hua-Feng-Dan upregulated genes of cellular function, signal regulation, circadian, and signal regulation; reduced the expression of genes related to Phase-I, II and III xenobiotic metabolism, endogenous metabolism, disease traits, and immune modulation, all point towards adaptive responses. However, the dose and duration of administration could differentiate the adaptive response from toxicity [16][17][18]. Yaomu as the "Guide Drug" of Hua-Feng-Dan [14,15] produced similar effects, but to a lesser extent. The change pattern of DEGs in low-dose Yaomu (0.1 g/kg) is more similar to that of Hua-Fng-Dan, while high-dose Yaomu (0.3 g/kg) is less effective, suggesting appropriate dose of Yaomu could be important for pharmacological effects of Hua-Feng-Dan.
Based on Two-dimensional clustering and fold changes, 10 DEGs were selected for qPCR veri cation. In general, the qPCR results matched the RNA-Seq results. Hua-Feng-Dan induced Cyp2a4 (2.2-fold), Gm3776 (6.5-fold), Mt2 (10-fold), Gsta1 (4.7-fold), Il1rn (24-fold) and Egr1 (8.9-fold). Cyp2a4 is sensitive to circadian regulation and plays roles in liver detoxi cation [30]; Overexpression of Cyp4a14 induced by hyperoxia contributes to increased resistance to oxidative stress [31]; Fgf21 is a metabolic regulator with multiple physiologic functions and is a novel biomarker for non-alcoholic fatty liver disease (NAFLD) in humans and limits hepatotoxicity in mice [32]; Gm3776 is a CAR-target gene in liver of mice as an adaptive responses to xenobiotics [33]; Metallothioneins (Mt2) are sulfhydryl-rich small proteins for heavy metal detoxi cation and free radical scavenging as a cellular antioxidant [34]; Glutathione S-transferase A1 (Gsta1) is a phase II enzymes playing adaptive against toxic stimuli [35]; Interleukin-1 receptor antagonist (Il1rn) is a member of the interleukin 1 cytokine family, inhibits the activities of interleukin 1α (IL-1α) and interleukin 1β (IL1-β), and prevent IL-1β overexpression during in ammatory responses [36]; Early growth response 1 (Egr1) is induced during liver injury and is an important acute-phase adaptive protein [37]; Tribbles Pseudokinase 3 (Trib3) plays an important roles in cellular stress and metabolism, and its induction is considered as an adaptive response [38]; Bcl2a1d is a member of the bcl2 protein family and acts as anti-and pro-apoptotic regulators that are involved in a wide variety of cellular activities [39]. Up-regulation of these genes points towards the adaptive machinery is activated by Hua-Feng-Dan and low dose of Yaomu to a less extent, in an attempt to program the liver to produce pharmacological effects.
IPA analysis of upstream regulatory factors helped our further understanding of Hua-Feng-Dan induced DEGs and pathway alterations. The top 25 upstream regulators clearly indicated the activation of the MAPK signaling pathways. For example, Transforming growth factor beta (Tgf beta) regulates cell fate decisions during embryonic development, tissue homeostasis and regeneration [40]. Neuregulin 1 (NRG1) is a cell adhesion molecule that plays an important role in regulating tissue development, growth and homeostasis [41]. Vascular endothelial growth factor (Vegf) is a sub-family of growth factors, and mediates liver and ischemic injury [42]. Nuclear factor kapa B (NF-κB) is a protein complex that controls transcription of DNA, cytokine production and cell survival, and plays an important role in adaptive immune response [43]. Signal transducer and activator of transcription 3 (STAT3) regulates hepatocyte proliferation and homeostasis [44]. P38 mitogen-activated protein kinases are a class of mitogenactivated protein kinases (MAPKs) that are responsive to stress stimuli [45]. ERK1/2 are broadly expressed protein kinases that are involved in proliferation, differentiation, motility, and cell death [46]. These upstream molecules imply that the adaptive machinery is activated. On the other hand, the downregulation of the PI3K inhibitor LY294002 [47], p38 MAPK inhibitor SB203580 [48], and MEK1/2 inhibitor U0126 [49] further facilitated the activation of these adaptive responses. A network pharmacology analysis of anti-stroke effects of Hua-Feng-Dan revealed potential molecular targets [7], and some of the regulated molecules coincide with the adaptive responses induced by Hua-Feng-Dan.
To further identify biological signi cance of Hua-Feng-Dan-induced DEGs, the keyword "treatment vs control" and "mice" was used to curate biosets in the GEO database. The -log(p-value) was used to evaluate the correlation of the database with Hua-Feng-Dan-induced DEGs [22]. The larger the -log(pvalue), the higher degree of similarity. For example, treatment of CD-1 mice for 4 days with 1,4dichlorobenzene (GSE 44783) [50] had a -log(p-value) of 26.9 for Hua-Feng-Dan; Fed mice with soybean oil-rich high fat diet (GSE68360) [51] had a -log(p-value) 26.4 for Hua-Feng-Dan; Fed mice with luteolin could reprogram the liver to ameliorate deleterious effects of diet-induced obesity (GSE54189) [52] had alog(p-value) of 20.7 with Hua-Feng-Dan. Glucose sensing is required to preserve B cell glucose competence (GSE43581) [53] had a -log(p-value) of 17.3 with Hua-Feng-Dan. Effects of unfolded proteins in the endoplasmic reticulum (ER stress) on glucose intolerance in eIF2 (alphaP) transgenic mice gene pro ling (GSE11210) [54] had a -log(p-value) of 14.5 with Hua-Feng-Dan. Fed Clock mutant mice with beza brate, a PPARα activator, affected circadian clock, lipid metabolism and liver gene expression (GSE20512) [55] had a -log(p-value) of 14.2 with Hua-Feng-Dan. Non-genotoxic carcinogen and peroxisome proliferator diethylhexylphthalate (DEHP)-treated mice showed hepatic gene expression changes associated with multiple biological functions, including complement activation, hemostasis, the endoplasmic reticulum overload response, and circadian rhythm (GSE55733) [56] had a -log(p-value) of 12.9 with Hua-Feng-Dan. Thus, Hua-Feng-Dan at clinical dose evoked hepatic gene expression changes towards adaptive responses, which was mimicked by the low dose of Yaomu to a lesser extent, but high dose of Yaomu was ineffective. These changes could be bene cial or detrimental, similar to the ability of oleanolic acid to program the liver to protect against various hepatotoxicants, but is hepatotoxic at high doses [18]. Caution should be taken when using Hua-Feng-Dan at the high dose and long-term administration.

Conclusions
This study pro led Hua-Feng-Dan and its "Guide Drug" Yaomu on liver gene expression by RNA-Seq. Hua-Feng-Dan and Yaomu at clinical dose did not cause liver toxicity, but induced differential gene expression re ecting adaptive responses to program the liver. Low dose of Yaomu produced similar, but less changes compared to Hua-Feng-Dan, while the high dose of Yaomu had little effect, suggesting that appropriate dose of Yaomu is important. Effects of Hua-Feng-Dan on the liver transcriptome changes were signi cantly correlated with chemical-induced adaptive responses in the livers of mice from the GEO database. Thus, Hua-Feng-Dan could induce adaptive responses to reprogram the liver to produce biological effects, which could be bene cial at appropriate doses, or harmful with high doses and longer time of administration.       Ingenuity Pathways Analysis (IPA) of upstream regulators based on HFD in Z-score. Red indicates the upregulation, and blue indicates the downregulation.