Comprehensive RNA-Seq Analysis of Potential Therapeutic Targets of Gan–Dou–Fu–Mu Decoction for Treatment of Wilson Disease Using a Toxic Milk Mouse Model

Background: Gan–Dou–Fu–Mu decoction (GDFMD) improves liver fibrosis in experimental and clinical studies including those on toxic mouse model of Wilson disease (Model). However, the mechanisms underlying the effect of GDFMD have not been characterized. Herein, we deciphered the potential therapeutic targets of GDFMD using transcriptome analysis. Methods: We constructed a tx-j Wilson disease (WD) mouse model, and assessed the effect of GDFMD on the liver of model mice by hematoxylin and eosin, Masson, and immunohistochemical staining. Subsequently, we identified differentially expressed genes (DEGs) that were upregulated in the Model (Model vs. control) and those that were downregulated upon GDFMD treatment (compared to the Model) using RNA-sequencing (RNA-Seq). Biological functions and signaling pathways in which the DEGs were involved were determined by gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses. A protein–protein interaction (PPI) network was constructed using the STRING database, and the modules were identified using MCODE plugin with the Cytoscape software. Several genes identified in the RNA-Seq analysis were validated by real-time quantitative PCR. Results: Total of 2124 DEGs were screened through the Model vs. control and Model vs. GDFMD comparisons, and dozens of GO and KEGG pathway terms modulated by GDFMD were identified. Dozens of pathways involved in metabolism (including metabolic processes for organic acids, carboxylic acids, monocarboxylic acids, lipids, fatty acids, cellular lipids, steroids, alcohols, eicosanoids, long-chain fatty acids), immune and inflammatory response (such as complement and coagulation cascades, cytokine–cytokine receptor interaction, inflammatory mediator regulation of TRP channels, antigen processing and presentation, T-cell receptor signaling pathway), liver fibrosis (such as ECM-receptor interactions), and cell death (PI3K-Akt signaling pathway, apoptosis, TGF-beta signaling pathway, etc.) were identified as potential targets of GDFMD in the Model. Some hub genes and four modules were identified in the PPI network. The results of real-time quantitative PCR analysis were consistent with those of RNA-Seq analysis. Conclusions: We performed gene expression profiling of GDFMD-treated WD model mice using RNA-Seq analysis and found the genes, pathways, and processes effected by the treatment. Our study provides a theoretical basis to prevent liver fibrosis resulting from WD using GDFMD.


INTRODUCTION
Wilson disease (WD), also known as hepatolenticular degeneration, is a rare autosomal recessive genetic disease of the nervous system related to disorders of copper metabolism. The clinical manifestations of this disease include liver damage, extrapyramidal symptoms, pigmented corneal rings, and kidney damage. The worldwide incidence rate is about 1 in 30,000 individuals (Bandmann et al., 2015). It is reported that 40-60% of patients with WD have liver disease as the first symptom, and about 50% of patients eventually develop liver cirrhosis (Ferenci et al., 2019). The liver is the central organ where copper is metabolized, and liver fibrosis is the most basic pathological change in the liver caused by disorders of copper metabolism (Bandmann et al., 2015). Liver fibrosis in WD is a pathological process in which the extracellular matrix (ECM) is reversibly deposited in the liver after chronic injury caused by copper ions. It is the main pathological change in almost every WD patient, and is an essential stage in the progression to cirrhosis (European Association for Study of Liver, 2012; Gerosa et al., 2019). Clinically, the main cause of death of patients with WD is liver cirrhosis and complications associated with it. Treatment of such patients with metal chelating agents often needs to be stopped because of severe side effects (Gerosa et al., 2019). This is a clinical problem that needs to be solved urgently and is, therefore, a focus of research globally.
Despite concerns about the need to protect the liver of patients with WD, no major breakthrough has been achieved, and treatment of such patients is still very difficult. It is, therefore, important to find key molecular targets for prevention of liver fibrosis in WD. Both prevention and treatment have great clinical significance. Because the mechanism underlying the pathogenesis of liver fibrosis is complicated, the repair response after liver injury involves the whole body; for this reason, drugs developed for a single target are seldom clinically effective. There are no chemical or biological drugs with prominent curative effects for clinical application. In recent years, research on traditional Chinese medicine (TCM) and its use has revealed therapeutic advantages of these medicines in the prevention and treatment of liver fibrosis (Chen et al., 2016).
In China, Gan-Dou-Fu-Mu decoction (GDFMD), a formula of TCM herbs, has been clinically used to treat WD (Yang et al., 2014). It can improve liver function and liver fibrosis indices in patients with WD exhibiting liver fibrosis. This formula consists of seven commonly used herbs, namely Reynoutria multiflora (Thunb.) Moldenke (Zhiheshouwu in Chinese, ZHSW), Lycium barbarum L (Gouqi in Chinese, GQ), Panax notoginseng (Burkill) F.H.Chen (Sanqi in Chinese, SQ), Curcuma aromatica Salisb (Yujin in Chinese, YJ), Smilax glabra Roxb (Tufuling in Chinese, TFL), Bupleurum chinense DC (Chaihu in Chinese, CH), and Paeonia lactiflora Pall (Baishao in Chinese, BS) . The anti-liver fibrosis effect of GDFMD has been reported to involve regulation of the expression of factors related to the TGF-β1/Smad signaling pathway at gene and protein levels (Yang et al., 2018). However, the specific molecular mechanisms underlying the effects of GDFMD are not completely understood. Therefore, in the present study, we tried to comprehensively elucidate the potential molecular mechanisms underlying the effects of GDFMD on a mouse model of WD using RNA-Seq technology.

HPLC-MS/MS Analysis of GDFMD
The proportion of the seven herbs in GDFMD was as follows: ZHSW:GQ:SQ:YJ:TFL:CH:BS 4:20:3:12:12:12:15. All the herbs were obtained from Anqing Huashi Chinese Herbal Medicine Beverage Co., Ltd (Anqing City, China). The steps in the preparation of GDFMD were as follows: First, all the herbal medicines were immersed in water (10-times the weight of medicine for 1 h, and then decocted by boiling for 2 h. After filtration, the residue was decocted with water (8-times the weight of residue) for 1.5 h. The two decoctions were combined and concentrated to 0.2 gmL −1 in vacuum at 50°C to obtain the GDFMD suspension. Thereafter, a 1.0 ml aliquot of GDFMD was added to 9.0 ml ethanol, and the mixture was centrifuged at 8000 × g for 5 min and filtered through a 0.22 µm membrane. Seven reference substances, namely gallic acid, resveratrol, quercetin, curcumin, paeoniflorin, saikosaponin A, and notoginsenoside (the purity of each compound was determined to be more than 98%), were obtained from Mansite Biotech Co., Ltd (Chengdu, China). The reference substances were weighed and dissolved in 90% ethanol to prepare a mixed reference solution. Two microliter aliquots of sample solutions and the mixed reference solution were used for LC-MS/MS analysis.
An Agilent 1290 HPLC system, coupled with an Agilent 6460 series QqQ/MS with electrospray ionization (ESI) source, was employed as the HPLC-MS/MS system for the quality control of GDFMD. The separation of complex components was performed on a BEH C18 column (100 mm × 2.1 mm, 1.7 µm) maintained at 35°C. Formic acid (0.1%) in water (A) and acetonitrile (B) were used as elution solvents. The gradient elution procedure was as follows: 0-30 min, 5-25% B; 30-40 min, 25-45% B; and 40.01-45 min, 5% B. The flow rate was 0.3 mlmin −1 . The mass spectrometry was carried by multiple reaction monitoring (MRM) in the ESI-mode. The optimum values of MS parameters were as follows: nebulizing gas pressure, 25 psi; drying gas (N 2 ) flow rate, 10 Lmin −1 ; capillary temperature, 350°C. All the data were analyzed using the Mass Hunter workstation software (Agilent Technologies, United States).

Animal Experiments and Sample Collection
Male Toxic Milk (TX) mice (verified specific-pathogen-free, 20-30 g) were purchased from the Jackson Laboratory (USA). All the mice were genotyped by Sanger sequencing. The mice were randomly divided into the following four groups: normal, model, GDFMD, and penicillamine groups, with 10 mice per group. The mice were anesthetized by intraperitoneal injection of sodium pentobarbital to minimize pain during surgery. The TX model was prepared as described previously (Reed et al., 2018) and the procedure complied with the National Institutes of Health Animal Care and Welfare Guidelines (National Institutes of Health Publication 80-23). The preparation of GDFMD and its administration (6.96 g/kg, once a day) to animals was done as previously described by us (Tang et al., 2019). Four weeks after the interventions, the mice were sacrificed and their liver was excised and promptly placed in liquid nitrogen for cryopreservation. The experimental protocol was authorized by the Animal Experiment Ethics Committee of the Anhui University of Chinese Medicine (2018AH-08).

Histopathological Examination
The liver was harvested at day 28 after GDFMD administration and fixed with 4% paraformaldehyde. Paraffin-embedded sections were stained with hematoxylin and eosin. Furthermore, we performed Masson staining of the sections and immunohistochemical staining for Col1 and α-SMA.

RNA Extraction, cDNA Library Preparation, Next-Generation Sequencing
A total RNA isolation kit (TR205-200, Tianmo, CN) was used for extracting the total RNA from liver samples, according to the instructions provided by the manufacturer. The quantity and quality of RNA samples were determined using a NanoDrop one (Thermo Fisher Scientific, United States) and Agilent 2100 Bioanalyzer (Agilent Technologies, United States). The RNA library was prepared using the SureSelect chain-specific RNA library preparation kit (Agilent Technologies, United States). Thereafter, Qubit 3.0 Fluorometer and Agilent 2100 Bioanalyzer were used to quantify the purified libraries. The generation of clusters was performed using cBot, and the sequencing was achieved with Illumina NovaSeq 6000 (Illumina, United States). These procedures were outsourced to Origin-Biotech Inc (Ao-Ji Biotech, Shanghai, China).

Analysis of Gene Profile Related to the Regulation of WD by GDFMD
FastQC was used for quality control of experimental samples to ensure the quality of RNA-Seq reads (version. 0.11.3). Fastp was used to trim the discovered Illumina TruSeq adaptor sequences, bad quality reads, and ribosomal RNA reads . The mapping of trimmed reads against the Mus musculus reference genome (mm10) was carried out using Hisat (Kim et al., 2015;Pertea et al., 2016). Thereafter, the quantification of each gene was done with read counts from the trimmed reads using the Stringtie (version: 1.3.0) software (Pertea et al., 2015;Pertea et al., 2016). Trimmed means of M values were normalized with gene counts (Robinson et al., 2010), and a Perl script was used to estimate the fragments per kilobase of transcript per million mapped reads (FPKM) (Mortazavi et al., 2008). To visualize the differences in gene expression in the liver tissue from control, tx-j, and GDFMD mice, we performed principal components analysis (PCA) using RNA-Seq expression data. EdgeR was used to analyze and determine the differentially expressed genes (DEGs) (Robinson et al., 2010;Nikolayeva and Robinson, 2014), and significance was determined by setting the following criteria: p value <0.05 and |log 2 (foldchange)| > 1 (Benjamini et al., 2001). The upregulated (up-DEGs; the intersection of Model and control or GDFMD, respectively) and downregulated (down-DEGs; the intersection of Model and control and GDFMD, respectively) genes were filtered out using Venny.

Functional Classification of Genes Involved in the Regulation of WD by GDFMD
To better understand the functions of DEGs, the R software package was used to analyze the signaling pathways enriched in Gene Ontology (GO) (Ashburner et al., 2000) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000). We used clusterProfiler to analyze the GO terms and KEGG pathways, and top 30 GO terms and pathways are listed (Yu et al., 2012).

Construction of PPI Network and Module Analysis
STRING, a database of information on protein interactions including the data on predicted and experimentally determined interactions (Szklarczyk et al., 2017), was adopted for determining the interactions of proteins encoded by genes involved in the regulation of WD by GDFMD. Based on interactions with combined scores ≥0.7, PPIs for the DEGs were generated using the STRING tool. Subsequently, visualization of the PPI network was done with Cytoscape (Shannon et al., 2003), and disconnnected nodes have been hidden. Using default conditions for the functional enrichment analysis module, the PPI network was used for screening the module based on the MCOD plugin in Cytoscape (Bader and Hogue, 2003). The clusterProfiler was used to analyze the DEGs enriched in GO and KEGG in the four groups (Yu et al., 2012).

Validation of Genes Involved in the Regulation of WD by GDFMD Using qRT-PCR
For verifying the reproducibility and reliability of the results of RNA-Seq, qRT-PCR was performed, using GAPDH as an internal control. Six GDFMD-regulated genes, also known to participate in liver injury, namely Timp1, Fbn1, Gas6, Alb, Apob, and Apoa2, were selected for qRT-PCR analysis. The expression levels in the liver samples from the control, model, and GDFMD groups were analyzed. The relative expression of mRNAs was determined using the 2 −ΔΔCT method (Livak and Schmittgen, 2001).

Statistical Analysis
Raw data are expressed as means ± standard deviation (SD). Statistical analysis of the data was performed using the GraphPad Prism (GraphPad Software, United States). A p-value less than 0.05 indicated statistical significance.

Identification and Determination of the Main Chemical Components of GDFMD
To ensure the quality and stability of the GDFMD, HPLC-MS/MS was used to detect and identify the chemical components of GDFMD.  Figure S1.

Effects of GDFMD on Mouse Model of WD
First, all the mice were genotyped. The mice with diploid mutation were tx-j mice (Model group), whereas the diploid wild type mice were used as control mice (control group) ( Figure 1A). To observe the protective effects of GDFMD and penicillamine against WD, TX mice were treated with GDFMD and penicillamine, administered intragastrically. The results indicated that WD induced liver injury in the mouse model. The results of HE staining ( Figure 1B) indicate that the liver was healthy in the control group, whereas diffuse lesions in the liver were observed in the WD model. A larger number of hepatic parenchyma cells were observed to be necrotic in the tx-j group compared with that in the control. GDFMD and penicillamine could effectively improve these lesions.
Collagen was stained blue in liver sections with Masson staining as well as with IHC staining for Col1a. The TX mice showed obvious collagen accumulation compared with those in the control group ( Figure 2B), and administration of GDFMD could reverse this. α-SMA is a marker of hepatic stellate cell activation, which has been widely used in the study of liver diseases. We, therefore, performed IHC staining for α-SMA. As shown in Figure 1B, the expression of α-SMA was upregulated in the tx-j model group compared with that in the control, and was downregulated in the penicillamine and GDFMD groups. This evidence proves that GDFMD can alleviate the death of liver cells, inflammatory cell infiltration, and liver fibrosis, just like penicillamine, in WD mice.

Identification of DEGs
To study the mechanism of action of GDFMD in the mice model, we performed RNA-Seq analysis and compared the transcriptomic profiles of the control, Model, and GDFMD groups. PCA showed that control, tx-j, and GDFMD sample groups could segregate into different areas of the 2D plots ( Figure 2A). As shown in Figure 2B, 611 down-DEGs that were common among the 1351 downregulated mRNAs in the Model vs. control comparison and 752 downregulated mRNAs in the Model vs. GDFMD comparison were filtered. Similarly, 1513 up-DEGs were identified to be common among the 5346

GO Analysis of DEGs Altered by GDFMD
To further identify the 611 down-DEGs and 1513 up-DEGs related to GDFMD treatment in the Model group, GO enrichment analysis was performed using the clusterProfiler. The results of enrichment analysis showed that the 611 down-DEGs were significantly (p < 0.05) enriched in 1553 GO terms, 511 of which were related to biological processes (BP), 43 were related to cellular components (CC), and 152 were related to molecular functions (MF). The top five enriched BP terms were organic acid metabolic process (GO: 0006082), carboxylic acid metabolic process (GO: 0019752), monocarboxylic acid metabolic process (GO:0032787), lipid metabolic process (GO: 0006629), and oxidation-reduction process (GO:0055114). Blood microparticle (GO: 0072562), extracellular space (GO: 0005615), endoplasmic reticulum (GO: 0005783), mitochondrion (GO: 0005739), and extracellular region (GO: 0005576) were the top five enriched CC terms. The top five enriched MF terms were oxidoreductase activity (GO: 0016491), catalytic activity (GO: 0003824), oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (GO: 0016705), monooxygenase activity (GO: 0004497), and tetrapyrrole binding (GO: 0046906). The top 30 GO terms related to down-DEGs, according to enrichment factor, are listed in Figure 3A.
These results suggest that GDFMD may regulate the biology process of metabolic process (eg, metabolic process of organic cyclic compound, carboxylic acid, monocarboxylic acid, lipid, etc), response to stimulus (eg, response to stress, organic substance, external stimulus, hormone, lipid, organic cyclic compound, wounding, cytokine, growth factor, etc), immune and inflammatory response, etc. GDFMD can affect almost all cell components, including intracellular, membrane-bounded organelle, membrane, cytoplasm, cell periphery, vesicle, extracellular vesicle, nucleus, cytoskeleton, etc.

KEGG Pathway Analysis of DEGs Altered by GDFMD
More information pertaining to the function of genes and interactions among them can be obtained through pathway enrichment analysis. ClusterProfiler was used for KEGG pathway analysis of 611 down-DEGs and 1513 up-DEGs, which were related to the GDFMD intervention in the Model group. The 611 down-DEGs were enriched in 230 pathway terms, including 65 significant terms with p < 0.05. In particular, the top10 enriched pathways were metabolic pathways (mmu01100), steroid hormone biosynthesis (mmu00140), retinol metabolism (mmu00830), complement and coagulation cascades (mmu04610), chemical carcinogenesis (mmu05204), drug metabolism -other enzymes (mmu00983), glycine, serine, and threonine metabolism (mmu00260), primary bile acid biosynthesis (mmu00120), metabolism of xenobiotics by cytochrome P450 (mmu00980), and peroxisome (mmu04146). The top 30 KEGG pathway terms related to down-DEGs, according to the FDR-value, are shown in Figure 4A.
Thus, dozens of pathways, especially metabolic regulation, immune regulation, extracellular matrix, cell death, and other important pathways, can be affected by GDFMD.
Using the MCODE plugin with default criteria, we obtained 39 modules, and these modules are presented based on MCODE scores in descending order. Following the MCODE scores, the top4 modules were designated as modules 1 to 4. These four modules were selected for visualization of the module network ( Figure 6). Specifically, module 1 consists of 49 genes and 880 interaction pairs, module 2 consists of 51 genes and 598 interaction pairs, module 3 consists of 36 genes and 222 interaction pairs, and the module 4 consists of 12 genes and 65 interaction pairs.

qRT-PCR Verification of the DEGs
Compared with the CN group, the results of qPCR demonstrated that Timp1, Fbn1, and Gas6 were overexpressed in the Model group (Figure 7). In contrast, the expression levels of Alb, Apob, and Apoa2 were low in the Model group as assessed by the two methods (Figure 7). GDFMD reversed the expression of these DEGs to different extent in the Model group (Figure 7). Thus, the results of RNA-Seq and real-time PCR analyses were consistent, indicating the reliability of the RNA-Seq results. The primers used for qRT-PCR analysis are listed in Table 3.

DISCUSSION
Several studies have shown that GDFMD treatment can prevent WD in model mice and has clinical application (Yang et al., 2014;Tang et al., 2018;Yang et al., 2018;Tang et al., 2019). In the present study also, we found that GDFMD treatment attenuated the liver injury caused by WD. However, the molecular mechanisms are unclear and require further studies. In the present study, the results of RNA-Seq analysis show that the GDFMD decoction could improve the phenotypic characteristics of WD and alter the expression of several genes. Furthermore, using GO, KEGG pathway, and PPI network analyses, we noticed that most of the genes altered by GDFMD were related to metabolism (including metabolic processes for organic acids, carboxylic acids, monocarboxylic acids, lipids, fatty acids, cellular lipids, steroids, alcohols, eicosanoids, and longchain fatty acids), immune and inflammatory response (such as complement and coagulation cascades, cytokinecytokine receptor interaction, inflammatory mediator regulation of TRP channels, antigen processing and presentation, and T-cell receptor signaling pathway), liver fibrosis (such as ECM-receptor interactions), and cell death (PI3K-Akt signaling, apoptosis, and TGF-beta signaling pathways).  The GO enrichment analysis revealed dozens of metabolic GO terms that were significantly enriched; these included organic acid metabolic process, carboxylic acid metabolic process, monocarboxylic acid metabolic process, lipid  metabolic process, fatty acid metabolic process, cellular lipid metabolic process, steroid metabolic process, alcohol metabolic process, eicosanoid metabolic process, and long-chain fatty acid metabolic process. With KEGG pathway analysis, we also found dozens of metabolic pathways that were significantly enriched, including steroid hormone biosynthesis, metabolic pathways, retinol metabolism, arachidonic acid metabolism, glycine, serine, and threonine metabolism, primary bile acid biosynthesis, drug metabolism-other enzymes, metabolism of xenobiotics by cytochrome P450, linoleic acid metabolism, steroid biosynthesis, and starch and sucrose metabolism. Notably, the PPAR signaling pathway, which is involved in the metabolism of fatty acids, bile acids, lipids, glycerophospholipids, glucose, and other substances, was also significantly enriched (Ceni et al., 2014;Li et al., 2015;Liu et al., 2015;Pawlak et al., 2015). Specifically, we found that GDFMD could effectively reverse the abnormal expression of 20 genes in WD; among these the expression of 18 genes (Apoa2, Acaa1b, Acsl1, Acox1, Ppara, Fads2, Scd2, Slc27a2, Pck1, Hmgcs2, Cyp7a1, Cyp4a31, Cyp4a14, Slc27a5, Slc27a1, Scd1, Cyp8b1, and Fabp1) was low and that of 2 genes (Cyp4a32 and Cyp4a10) was high in the Model group. A large body of evidence suggests that PPAR signaling pathway is involved in different liver diseases (Panebianco et al., 2017;Zhang et al., 2020), especially liver fibrosis (Chen et al., 2015;Pawlak et al., 2015;Chhimwal et al., 2020), as was also confirmed in the present study. Our results show that GDFMD can alleviate liver damage in WD by reversing the changes in the PPAR signaling pathway (Figure 8). GO enrichment analysis revealed that more than 327 DEGs are involved in the immune system process (p 2.156e−14) and more than 120 DEGs are involved in the inflammatory response (p 1.169e−12). It appears that both adaptive (p 2.397e−06) and innate (p 4.022e−06) immune responses are abnormal in WD. Several immune and inflammatory response pathways, such as complement and coagulation cascades, cytokine-cytokine receptor interaction, inflammatory mediator regulation of TRP channels, antigen processing and presentation, and T cell receptor signaling pathway, were enriched. We found that hundreds of immune factors, including 133 cytokines and cytokine receptors, 13 interleukins and interleukin receptors, 8 TNF family members and receptors, 10 TGF-β family members and receptors, and 34 chemokines and receptors, are dysregulated in WD, and that GDFMD could revert these dysregulations. Therefore, we hypothesize that WD can lead to the recruitment and induction of monocytes, T cells, eosinophils, and lymphocytes, among other immune cells. Previous studies have shown that the complement system is involved in a variety of chronic liver diseases (Bykov et al., 2007;Schmitt et al., 2012;Himoto et al., 2019). The activation of the complement and coagulation cascade leads to an inflammatory response, including degranulation, chemotaxis, and phagocytosis. Our study confirms these findings. We found that GDFMD inhibits the changes in the expression of C3ar1, C4R, and C5ar1 caused by WD. We, therefore, speculate that liver injury in WD can be inhibited by GDFMD through the regulation of the immune and inflammatory response, particularly of the complement and coagulation cascades (Figure 9). During liver injury, the expression of and changes in various ECM components in fibrous tissues may convey important pathophysiological information, which is valuable for the development of new biomarkers and antifibrosis interventions. We identified that more than 31 ECMreceptors, except for SDC4, VTN, and FN1, were upregulated in WD ( Figure 10). In particular, 11 collagens (COL1A1, COL1A2, COL4A3, COL4A4, COL4A5, COL4A6, COL6A1, COL6A2, COL6A3, COL6A5, and COL9A3), 5 integrins (ITGA11, ITGA3, ITGA8, ITGB6, and ITGB7), 6 laminins (LAMA2, LAMA4, LAMA5, LAMB1, LAMB2, and LAMC2), and 3 thrombospondins (THBS1, THBS2, and THBS3) were upregulated in WD, and GDFMD could effectively reverse their abnormal expression levels. Collagen has been used as an important marker of liver fibrosis, and a large number of drugs have been developed based on the regulation of collagen, which can play an effective role in inhibiting liver fibrosis (Cho et al., 2000;Williams et al., 2001;Galli et al., 2002;Fontana et al., 2008). We also found that integrins α3, α8, α11, β6, and β7 were upregulated in WD. Popov and Sullivan demonstrated the role of β6 in liver fibrosis (Popov et al., 2008;Sullivan et al., 2010). In general, our results indicate that there is a systemic abnormal expression of ECMreceptors during the pathogenesis of WD. Reversal of these abnormalities could definitely be one of the mechanisms for the treatment of liver injury in WD.
Besides, we also identified many cell death-related genes and pathways that could be affected by GDFMD, such as the PI3K-Akt signaling pathway, apoptosis, and TGF-beta signaling pathway. Sixteen genes related to TGF-beta signaling pathway, namely Ltbp1, Id3, Tgfb3, Thbs1, Dcn, Id4, Bambi, Tnf, Inhbc, Smad7, Bmp5, Inhbb, Gdf7, Bmp6, Id1, and Gdf6, could be enriched, consistent with our findings in a previous study (Yang et al., 2018). In conclusion, GDFMD may attenuate tx-j WD by altering metabolism, immune and tissue remodeling, and cell deathrelated pathways. This is a bioinformatics study to explore the unknown mechanism underlying the effect of a TCM. The results provide valuable information for further research and can provide a new model for development of novel drugs. Although most of our data are consistent with those reported in other studies, some limitations of this study are worth noting. First, the progression of WD may be related to age. In this study, DEGs were examined in adult male mice. Further investigation is needed to verify the candidate targets in different age groups. Second, in view of the regulation of protein translation and posttranslational modification, we do not know the extent to which altered gene expression affects protein function. The functions of proteins encoded by these genes need to be further studied to provide new ideas for the treatment of WD.

CONCLUSION
Through genetic engineering mouse model and RNA-seq, we preliminarily verified the efficacy and mechanism of GDFMD in the treatment of WD. The pathogenesis of WD and the treatment mechanism of GDFMD were interpreted from the whole genome level. However, further experimental is needed to verify for identified genes and pathways.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://bigd.big.ac. cn/bioproject/browse/PRJCA004081.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Experiment Ethics Committee of the Anhui University of Chinese Medicine (2018AH-08).

AUTHOR CONTRIBUTIONS
TW and WH contributed equally to this work. TW and WY conceived and designed the study. TW, LT, and XD performed the HPLC-MS/ MS analysis. TW, LT, WH, YY, and NQ performed the animal experiments. TW and WH performed histopathological examination and qPCR. TW, HW, SH, and JL analyzed the data and pre-viewed the manuscript. All authors read and approved the final manuscript.