TDF was kindly gifted by Guangshengtang Co., Ltd (Fujian, China). The purity of TDF was 99.5%.

The study protocol was approved by the Animal Protection Organization and Ethics Committee of Chongqing Medical University (Chongqing, China). Female C57BL/6J mice (8 weeks) from the Animal Center of Chongqing Medical University were housed in a room with a 12-h light and dark cycle at 22°C with free access to mouse chow and water. After 1 week of acclimatization, mice were divided randomly into two groups. Mice in the TDF group were administered with TDF solution (455 mg of TDF powder + 0.5 g of sodium carboxymethyl cellulose +100 ml of normal saline were mixed thoroughly with a homogenizer until the solution was transparent) at 45.5 mg/kg/day by oral gavage for 4 months. Dose determination of TDF was performed via conversion of human equivalent doses to murine doses based on the body surface area (Reagan‐Shaw et al., 2008; Ng et al., 2015). Mice in the control group received an equivalent volume of vehicle (100 ml of normal saline + 0.5 g of sodium carboxymethyl cellulose were mixed thoroughly with a homogenizer until the solution was transparent) for 4 months. Each mouse was weighed once a week. At study termination, mice were killed after 12 h of fasting. Blood samples were taken by excising the eyeballs. The liver was collected and weighed for subsequent experiments. A flowchart of this study is shown in Figure 1.

Serum was collected after centrifugation at 1,000 rpm for 15 min. Serum ALT levels were measured in the Clinical Laboratory, The Second Affiliated Hospital, Chongqing Medical University.

Liver tissues of mice were fixed in 4% paraformaldehyde, embedded in paraffin, sliced, and stained with hematoxylin and eosin (H&E).

According to manufacturer instructions, TRIzol^® (Invitrogen, United States) was used to extract total RNA from the liver. Then a bioanalyzer (2,200 series; Agilent Technologies, United States) was used to evaluate the concentration, purity, and quality of RNA. Then RNA was sequenced using the DNBseq platform in BGI (Shenzhen, China).

Low-quality reads, adaptor reads, and reads with > 10% unknown bases (poly-N) were removed from raw data to obtain high-quality “clean” reads for subsequent analyses. The Q20 (percentage of bases with a quality value ≥ 20) and Q30 content of clean data were also calculated. Clean reads were mapped to the Mus musculus reference genome (GRCm38.p6) using HISAT2 (http://daehwankimlab.github.io/hisat2/) (Kim et al., 2015). Expression was calculated by RSEM (Li and Dewey, 2011) and represented in fragments per kilobase per million (FPKM) reads. Differential gene expression was identified using the “DEGseq” package with R (R Institute for Statistical Computing, Vienna, Austria) (Wang et al., 2010). The absolute value of fold change (FC) ≥ 2 and adjusted p-value (Q-value) ≤ 0.001 were adopted as criteria for determining the significance of differential expression of a particular gene.

The Gene Ontology (GO) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) database were employed to identical enrichment of function and signaling pathways, respectively, based on Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) (Huang et al., 2009). Following the instructions of the DAVID manual, first, we clicked on the “Start Analysis” on the Internet website. Second, we entered the DEGs list, selected identifiers as “entrez gene ID,” selected list types as “gene lists,” and submitted lists. Third, we chose to limit annotations and background by M. musculus. Finally, the enrichment results of GO and KEGG databases were presented. p < 0.05 and gene count ≥ 2 were considered significant.

Furthermore, additional analyses of enrichment of function signaling pathways were done using Metascape (https://metascape.org/gp/index.html#/main/step1/) (Zhou et al., 2019). First, we pasted the gene list as “entrez gene ID.” Second, we chose to input the species as M. musculus. Third, we clicked on “Express Analysis.” Finally, analyses of enrichment of function and signaling pathways were carried out with the following ontology sources: Biological Process (BP) within the GO database, KEGG Pathway, Reactome Gene Sets, CORUM, TRRUST, PaGenBase, WikiPathways, and PANTHER Pathway. Terms with p < 0.01, minimum count of 3, and enrichment factor >1.5 (the enrichment factor is the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities.

First, Metascape was utilized to construct a PPI network and identify the significant modules. Besides, the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/), a user-friendly online system that provides predicted and experimental interactions of proteins (Szklarczyk et al., 2019), was used to establish a PPI network of DEGs with a confidence score ≥ 0.7 for significant differences. Then the PPI network was visualized using Cytoscape 3.6.1 (www.cytoscape.org) (Shannon, 2003). Molecular Complex Detection (MCODE) 1.5.1 (a plugin of Cytoscape) (Bader and Hogue, 2003) was used to screen and identify the most significant modules in the PPI network. MCODEs were extracted when the node score cutoff was 0.2 and K-core was 2. cytoHubba (a plugin of Cytoscape) was employed to calculate the properties of the network topology for nodes to identify hub genes with a degree ≥ 10. The “degree” indicates the number of edges connected with a specific node. Nodes with a high degree are identified as hub genes (i.e., may contribute to vital biological behaviors).

A Venn diagram was delineated to identify significant union genes among “Metascape_MCODE,” “Cytoscape_MCODE,” and “Cytoscape_cytoHubba” by Bioinformatics (www.bioinformatics.com.cn/), an online platform for the analyses and visualization of data. Summaries for the basic information of the significant genes were obtained via Mouse Genome Informatics (www.informatics.jax.org/). A hierarchical clustering heatmap of significant genes was plotted by using OriginPro 2021 based on gene expression, and classified by the biological function of genes. Correlation analyses among significant genes were achieved through Pearson’s correlation test.

DisGeNET 7.0 (www.disgenet.org/) is a discovery platform containing one of the largest publicly available collections of genes and variants associated to human diseases (Piñero et al., 2019). The current version of DisGeNET contains 1,134,942 gene–disease associations (GDAs), between 21,671 genes and 30,170 diseases, disorders, traits, and clinical or abnormal human phenotypes. The relationships between the significant genes and nine common CLDs were analyzed via this tool. First, we entered the name of the diseases in the “Search” box. Subsequently, the summary of the GDA score and evidence for GDAs were presented, and the results were downloaded. Finally, the target genes associated with CLDs were identified from the results.

To further explore how TDF improves CLDs, Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST 2, https://www.grnpedia.org/trrust/), a database containing 6552 TF–target interactions for 828 mouse TFs (Han et al., 2018), was employed to predict the TFs of target genes. First, we selected the species as “mouse,” and submitted the list of target genes in the bottom panel of the search area titled “Find key regulators for query genes.” Then the results were downloaded. Meanwhile, prediction of miRNAs of the target genes and TFs identified from TRRUST was done using DIANA-TarBase v8 (www.microrna.gr/tarbase/), a database that curates experimentally verified miRNA targets manually and contains 665,843 unique miRNA–target pairs (Karagkouni et al., 2018). We defined the species as M. musculus, entered the genes and TFs one-by-one, and tabulated the results. Ultimately, pairwise-related genes, TFs, and miRNAs were screened, and Cytoscape was utilized to visualize the TF–miRNA co-regulatory networks of the target genes.

TDF was administered to mice by oral gavage for 4 months to simulate the long-term use of TDF in humans. Before intervention, there was no significant difference in the body weight (BW) of mice (Supplementary Figure S1A) between the control group and TDF group. At study termination, no mice died from TDF administration. Besides, no significant differences were found between the TDF group (n = 9) and control group (n = 11) in BW (Supplementary Figure S1B), liver weight (LW) (Supplementary Figure S1C), liver index (LW/BW) (Supplementary Figure S1D), ALT level (Supplementary Figure S1E), and liver histology (Supplementary Figure S1F). Taken together, these results indicated that the TDF dose we employed was non-toxic, and had no significant effect on liver function of normal mice.

To explore the transcriptional changes in the liver induced by TDF administration, RNA-seq was done using the liver tissues of the mice: 639.99 million raw reads (Supplementary Table S1) were generated. After removing adaptors and low-quality reads, we obtained 616.74 million clean reads, with a high quality of Q30 ≥ 93.36%. Then the trimmed clean reads were mapped onto the M. musculus reference genome, and 76.32–81.35% of clean reads were mapped uniquely to the genome (Supplementary Table S1B). The uniquely mapped reads were used in all subsequent analyses.

A total of 17,199 mRNAs were annotated (Supplementary Table S2). DEGseq was employed to screen for DEGs. A total of 1,341 annotated genes were identified to be differentially expressed when considering exclusively a stringent threshold of Q-value ≤ 0.001 and an absolute value of FC ≥ 2, which is presented as a volcano plot (Figure 2A). After removing DEGs that caused differences between groups due to abnormal expression within the group, we obtained 854 DEGs eventually (217 downregulated DEGs and 637 upregulated DEGs) (Figure 2B and Supplementary Table S3). Hence, TDF could affect gene expression in mouse livers directly. To obtain a global view of these 854 DEGs, hierarchical clustering (Figure 2C) was done with normalized FPKM values, and indicated that our samples were of “good” quality with gene expression of similar proportion in each group.

First, enrichment analyses using the GO database and KEGG database were undertaken using DAVID. We discovered that 774 out of 854 profiled DEGs were assigned to 394 GO terms: 264 for biological process (BP), 32 for cellular component (CC), and 98 for molecular function (MF). Variations in DEGs related to BP were involved mainly in “immune system process,” “inflammatory response,” “lipid metabolic process,” and “glucose metabolic process” (Figure 3A). With regard to CC, DEGs were significantly enriched in the “extracellular region,” “extracellular space,” “organelle membrane,” and “cell surface” (Figure 3B). Variations in DEGs associated with MF were significantly enriched in “small molecule binding,” “insulin-activated receptor activity,” “iron ion binding,” and “chemokine activity” (Figure 3C). The KEGG database indicated that 356 out of 854 profiled DEGs were assigned to 46 signaling pathways. The significant pathways relevant to DEGs were immune, inflammatory pathways (“cytokine-cytokine receptor interaction,” “NOD-like receptor signaling pathway,” “chemokine signaling pathway,” and “TNF signaling pathway”), and metabolic pathways (“retinol metabolism,” “steroid hormone biosynthesis,” “arachidonic acid metabolism,” and “glutathione metabolism”) (Figure 3D). The data of GO and KEGG classifications are shown in Supplementary Table S4.

To gain further insight into the functions of DEGs, analyses of enrichment of signaling pathways and function were carried out via Metascape. The results (in accordance with the results using DAVID) indicated that the DEGs induced by TDF were significantly enriched in inflammatory processes (“leukocyte migration,” “neutrophil degranulation,” “acute inflammatory response,” “positive regulation of leukocyte migration,” and “regulation of interleukin-1 production”) and metabolic processes (“retinol metabolism” and “monocarboxylic acid metabolism”) (p < 0.05, Figures 3E–G).

First, a PPI network of DEGs was constructed through Metascape (Figure 4A). Fourteen MCODE modules were identified from the PPI network. Notably, MCODE1 with the highest score consisted of 42 genes: Sirpb1c, Gm9733, Sirpb1b, Gm5150, Sirpb1a, Fcgr4, Ticam2, Sucnr1, Siglece, P2ry13, Cd177, Aldh3b1, Cxcl13, Lair1, Atp11a, Tlr2, Sstr2, Sell, Cxcl5, Cxcl2, Ccl6, Ccl5, Ccl4, Pld1, Pirb, Mtnr1a, Clec4d, Mme, Lgals3, Fabp5, Itgb2, Itgam, Cxcl1, Fpr1, Fpr2, Fcer1g, Ccr2, Ccr1, Cxcr4, C5ar1, Adra2a, and Adam8 (Figure 4B).

Simultaneously, construction of the PPI network was also established by STRING with a confidence score of ≥ 0.7 for significant differences. There were 1937 edges and 443 nodes in the PPI network (PPI enrichment p-value <0.001) (Figure 4C). Thirty MCODE modules were identified from the PPI network by the Cytoscape_MCODE. In particular, MCODE1 with the highest score comprised 24 genes: Serpinb6b, Sirpb1c, Pirb, Aldh3b1, Sell, Lilrb4, C5ar1, Fabp5, Pira2, Mme, Tlr2, Siglece, Lgals3, Sirpb1b, H2-Bl, Fcgr4, Sirpb1a, Gp49a, Ticam2, Gm5150, Gm9733, Fpr1, Gm14548, and Serpinb10 (Figure 4D). With degree ≥ 10 considered as the standard of judgment, 10 genes were identified as hub genes with Cytoscape_cytoHubba: Fpr2, Cxcl2, Fpr1, Cxcl1, Tnf, C5ar1, Serpinb6b, Aldh3b1, Tlr2, and Itgam (Figure 4E).

Finally, a VENN diagram was delineated and showed 50 significant union genes among “Metascape_MCODE1,” “Cytoscape_MCODE1,” and “Cytoscape_cytoHubba,” including four common genes: Aldh3b1, Tlr2, Fpr1, and C5ar1 (Figure 5A). The basic information of the 50 significant genes is summarized in Supplementary Table S5. Hierarchical clustering indicated that the significant genes could largely differentiate the TDF group from the control group (Figure 5B). Moreover, we found that these 50 genes were involved mainly in “immunity” (44%, 22/50), “inflammation” (34%, 17/50), and “metabolic” processes (10%, 5/50), and six genes (12%, 6/50) exhibited other or undefined functions (Figure 5B). The Pearson correlation analysis showed a positive correlation among most genes except for Sucnr1, H2-Bl, Mme, and Cxcl5 (Figure 5C).

The DisGeNET platform was employed to identify the potential targets of TDF for alleviating CLDs. The nine most common CLDs including non-alcoholic fatty liver disease (NAFLD), cholestasis, primary biliary cholangitis (PBC), primary sclerosing cholangitis (PSC), autoimmune hepatitis (AIH), non-alcoholic steatohepatitis (NASH), liver fibrosis (LF), cirrhosis, and hepatocellular carcinoma (HCC) were analyzed. The gene–disease association score (GDAs) between the potential target genes and these CLDs are portrayed in Figure 6A. Nineteen genes were associated with these nine CLDs. The number of genes involved in “inflammatory,” “immune,” and “metabolic” processes was eight (42%), seven (37%), and three (16%), respectively. From the vertical perspective, 10 genes were related to NAFLD: Tnf, Ccr2, Tlr2, Lgals3, Ccl5, Cxcl, Sell, Lilrb4a, Lilr4b, and Fabp5. The four genes relevant to cholestasis were Tnf, Ccr2, Cxcl2, and Mme. Five genes (Tnf, Tlr2, Lgals3, Ccl5, and Itgb2) were connected with PBC. In addition, Tnf and Ccr2 were associated with PSC. Tnf along with Tlr2 were related to AIH. Moreover, seven genes were involved in NASH: Tnf, Ccr2, Tlr2, Lgals3, Cxcl5, Cxcr4, and Adra2a. Eight genes (Tnf, Tlr2, Lgals3, Ccl5, Cxcl2, Cxcr4, Ccr1, and Sucnr1) were relevant to LF. Besides, there were nine genes associated with cirrhosis: Tnf, Ccr2, Tlr2, Lgals3, Ccl5, Cxcl5, Adra2a, Itgam, and Cxcl1. Tnf together with Ccr2 were related to HCC. In addition, according to the GDA score, the gene most associated with PBC was chemokine C-C motif ligand 5 (Ccl5), and the gene most closely related to the other eight liver diseases was Tnf. From the lateral perspective, we found that tumor necrosis factor (Tnf) was related to all of the nine CLDs, with cirrhosis (GDA score: 0.4) and cholestasis (GDA score: 0.34) showing the closest associations. Chemokine C-C motif receptor 2 (Ccr2) and toll-like receptor 2 (Tlr2) were associated with six CLDs, and galectin 3 (Lgals3) was related to five CLDs. The log₂FC of genes indicated that 19 genes consisted of 17 upregulated DEGs and two downregulated DEGs, and the most significantly altered genes were adrenergic receptor alpha 2a (Adra2a), chemokine C-X-C motif ligand 1 (Cxcl1), and integrin alpha M (Itgam), with log₂FC values of 2.84, 2.83, and 2.76, respectively (Figure 6B). The basic information of these genes can be found in Supplementary Table S5.

First, TF–gene analyses were undertaken with TRRUST, and 13 upstream TFs targeting 12 target genes (seven genes had no results) were identified (Table 1). Quite specifically, NF-κB family (Nfkb1, Rela, and Rel) were the most significantly enriched TFs and targeted seven genes (Ccl5, Cxcl1, Cxcl2, Itgam, Tlr2, Tnf, and Cxcr4). Additionally, the RNA-seq data showed that TDF could also affect expression of Nfkb1 (log₂FC = 0.22, Q-value <0.001), Rela (log₂FC = 0.33, Q-value <0.001), and Rel (log₂FC = 0.23, Q-value = 0.29). Overall, these results indicated that NF-κB might be the most critical TF for the genes targeted by TDF to relieve CLDs.

Subsequently, DIANA-TarBase v8 was used to predict the upstream miRNAs of the 19 target genes, and 80 unique miRNAs were identified (Supplementary Table S6A). The most pivotal miRNAs are shown in Table 2. These results suggested that mmu-miR-155-5p was the most important miRNA and regulated 15 target genes.

Meanwhile, TF–miRNA interactions were predicted for the 13 identified TFs by using DIANA-TarBase, and 156 unique miRNAs were identified (Table S6B). Table 3 shows the most pivotal miRNAs. The results indicated that mmu-miR-155-5p and mmu-miR-124-3p, interacting with eight TFs, might be the most important miRNAs.

Finally, nine genes, 12 TFs, and 34 miRNAs that were pairwise-related were screened out to construct the TF–miRNA co-regulatory networks of the target genes, which consisted of 247 edges and 55 nodes (Figure 7A). These results indicated that mmu-miR-155-5p (involved in 15 edges) might be the most important miRNA. The most important TF might be Nfkb1 (19 edges). Their target genes were Tnf, Ccl5, Tlr2, Cxcl1, and Cxcl2 (Figure 7B).

Taken together, these results suggested that co-regulatory networks consisting of mmu-miR-155-5p and Nfkb1 might be (at least in part) the underlying mechanisms by which TDF improved CLDs.