Relevant microarray datasets associated with NAFLD were retrieved from the GEO microarray database (Clough et al., 2024) (GEO, https://www.ncbi.nlm.nih.gov/geo/) using the keywords “non-alcoholic fatty liver disease” and “NAFLD”. After comparison and screening, three datasets (GSE260666, GSE126848, and GSE135251) were ultimately selected for the study. Only human-liver RNA-seq studies providing raw counts, ≥10 samples and full NAFLD histology without drug-treatment confounders met the eligibility criteria. To minimise cross-study batch effects, differential expression was calculated inside each dataset with GEO2R; subsequent analyses used the intersection of significant genes across the three datasets, a strategy that bypasses direct data merging and reduces technical bias. Expression data from NAFLD liver tissues and normal liver tissues were extracted from the samples for comparative analysis. Differential expression was determined using a combined criterion of |log₂FC| ≥ 1 and an adjusted P value <0.05; genes meeting these thresholds were classified as significantly up- or downregulated accordingly. Finally, the intersection of differentially expressed genes from the three datasets was identified, and the shared genes were designated as the key differentially expressed genes (DEGs) for this study.

The differentially expressed genes were entered into the STRING database (Franceschini et al., 2013) (https://string-db.org/) to construct a protein-protein interaction (PPI) network, with an interaction threshold set at 0.400, which corresponds to STRING’s medium-confidence score and is widely used to achieve a balance between capturing biologically meaningful interactions and minimizing false positives in large-scale analyses, to identify proteins with significant interactions. Subsequently, the PPI network was visualized using Cytoscape software (version 3.9.1) (Shannon et al., 2003), and highly interconnected protein interaction modules were identified using the MCODE plugin to create corresponding visual network diagrams. The specific parameters were set as follows: degree = 2, node score cutoff = 0.2, k-core = 2, and maximum depth = 100.

The Omicshare online analysis tool (http://www.omicshare.com/tools/index.php/) was utilized to perform Gene Ontology (GO) secondary classification enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on differentially expressed genes. Enrichment significance was assessed with the hypergeometric test, and the resulting P values were corrected for multiple comparisons using the Benjamini–Hochberg false-discovery rate method; pathways with an adjusted P < 0.05 were considered significant. The species selected for this analysis was Homo sapiens, and the results included biological processes, molecular functions, and cellular components from the GO enrichment analysis, as well as related biological pathway enrichment outcomes.

The identified differentially expressed genes associated with NAFLD were entered into the DGIdb database (https://dgidb.org) (Cannon et al., 2024) to identify potential drugs that may target these core genes. Drugs predicted to exhibit known interaction types with the core genes of NAFLD are considered potential candidates for treating this condition. Subsequently, the interactions between the potential drugs and their corresponding target genes were visualized using Cytoscape software.

The differentially expressed genes were entered into the HERB database (Fang et al., 2021) (http://herb.ac.cn/) to filter out traditional Chinese medicine components that have a statistically significant mapping relationship (P < 0.05) with these genes. A statistical analysis was conducted on these traditional Chinese medicine components to examine the distribution of their meridian affinities and efficacy categories.

The three-dimensional (3D) structure of the target compound was downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Target protein structures were then selected from the Protein Data Bank based on the following criteria: a resolution below 2.5 Å and the presence of co-crystallized small-molecule ligands. Subsequently, the selected compounds and target proteins were submitted to the receptor-ligand docking platform DockThor (https://dockthor.lncc.br/v2/) for molecular docking analysis. The top-ranked docking results were visualized using PyMOL 2.4.0 software. The binding stability between the molecules was assessed by analyzing the docking affinity scores.

The gene expression data selected for this study includes: ① A study published in 2004 by Xiangya Second Hospital of Central South University, which analyzed the differential gene expression in liver tissues of patients with NAFLD compared to a control group. The dataset is identified as GSE260666 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE260666), comprising a total of 16 samples, including 10 NAFLD samples and 6 normal liver tissue samples; ② A transcriptomic analysis across the spectrum of non-alcoholic fatty liver disease published in 2020 by Newcastle University, with the dataset identified as GSE135251 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135251), comprising a total of 216 samples, including 206 samples at various stages of fibrosis and 10 normal liver tissue samples; ③ A study published in 2019 by Gubler in Denmark that compared the hepatic transcriptomic characteristics of patients with varying degrees of non-alcoholic fatty liver disease to those of healthy individuals with normal weight. The dataset is identified as GSE126848 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126848), comprising a total of 57 samples, including 14 normal healthy liver tissue samples, 12 normal healthy obese liver tissue samples, 15 samples of simple steatosis, and 16 samples of non-alcoholic steatohepatitis with or without fibrosis.

Using P < 0.05 and |log2FC| ≥ 1 as the criteria for selection, a total of 167 genes were identified in the dataset GSE260666, including 100 upregulated genes and 67 downregulated genes. In the dataset GSE135251, a total of 2,166 genes were identified, including 871 upregulated genes and 1,295 downregulated genes. In the dataset GSE126848, a total of 1,241 genes were identified, including 526 upregulated genes and 715 downregulated genes. By constructing a Venn diagram to analyze the intersection of upregulated and downregulated genes from the three datasets, a total of 34 differentially expressed genes were identified, including 16 upregulated and 18 downregulated genes (Figure 1).

A total of 34 differentially expressed genes were entered into the STRING database to construct the PPI network (Figure 2A), resulting in 22 interactions. The network was then visualized and analyzed for modularity using Cytoscape software (version 3.9.1) along with the MCODE plugin, identifying the most significant module, which is highlighted in red (Figure 2B). Further analysis revealed that CXCL2, CDKN1A, TNFRSF12A, and HGFAC are the key regulatory hub genes within this module of differentially expressed genes.

This study employed a significance threshold of P < 0.05 for selection, resulting in a total of 231 entries from the GO enrichment analysis. Among these, 161 entries pertain to biological processes, including the negative regulation of DNA biosynthesis, RNA splicing, and modulation of fibroblast proliferation. Forty-four entries relate to molecular functions, encompassing interleukin-33 binding, interleukin-1 receptor activity, and cyclin-dependent protein kinase activating kinase activity. Additionally, 26 entries correspond to cellular components, which include integrin α8-β1 complexes, growth factor complexes, and the extracellular environment. The top 20 results for biological processes, molecular functions, and cellular components are illustrated in bubble charts (Figures 3A–C).

The KEGG enrichment analysis identified a total of 81 pathways, of which 17 exhibited a significance level of P < 0.05, as detailed in Table 1. The results indicate that these pathways are enriched in several tumor-related processes, including cytokine interactions, the PI3K-Akt signaling pathway, and butyrate metabolism, with the melanoma pathway showing the strongest correlation. The results of the KEGG enrichment analysis are presented in a circular diagram (Figure 3D).

To explore potential therapeutic strategies, this study utilized the DGIdb database to predict drugs that target the aforementioned core genes (Figure 4). The selection criteria included drugs approved by the U.S. Food and Drug Administration, resulting in a list of chemical compounds (Table 2). A total of 22 drugs or compounds targeting the core genes CDKN1A and CXCL2 were identified. However, no drugs were predicted to potentially regulate the HGFAC gene. The U.S. Food and Drug Administration has not yet approved any drugs or compounds targeting TNFRSF12A. To validate its causal role in NAFLD pathogenesis, an initial step could involve pooled CRISPR-Cas9 knockout screens in human hepatocyte-like cells subjected to lipotoxic stress. Subsequent validation would include arrayed siRNA knock-down assays in hepatic organoids derived from NAFLD patients. Key convergent phenotypic readouts such as neutral lipid accumulation, NF-κB signaling activity, and pyroptotic cytokine release could then be assessed to confirm mechanistic involvement and identify potential assay formats suitable for future small-molecule discovery. Additionally, cyclosporine shows potential activity against CDKN1A and CXCL2, making it a particularly promising novel candidate compound for NAFLD treatment.

After inputting the core genes CXCL2 and CDKN1A, the search results from the HERB database identified five components: curcumin, resveratrol, ursolic acid, berberine, and tetrandrine. Screening through the TCMSP database (https://www.tcmsp-e.com) revealed that these components can be derived from various TCMs, including Alpinia officinarum (Gaoliangjiang), Curcuma longa (Jianghuang), Polygonum cuspidatum (Huzhang), Smilax glabra (Tufuling), Coptis chinensis (Huanglian), Stephania tetrandra (Fangji), and Cornus officinalis (Shanzhuyu).

This study employs molecular docking techniques to investigate the binding modes and characteristics of active components from traditional Chinese medicine with the core genes CXCL2 and CDKN1A. The three-dimensional structures of the proteins corresponding to the core genes were obtained from the Protein Data Bank database. As indicated in Table 3, the binding energies of all active components from traditional Chinese medicine with the target proteins are below −5 kcal/mol, suggesting a strong binding affinity between these components and the core targets (Wong et al., 2022). The visualization results of the molecular docking (Figures 5A–E) further reveal that these small molecules from traditional Chinese medicine can interact with key amino acid residues of the target proteins through hydrogen bonds, resulting in stable complex conformations. The findings of this study not only validate the potential mechanisms of action of these active components from traditional Chinese medicine at the molecular level concerning targets related to NAFLD but also provide a theoretical basis for their application in NAFLD treatment. Guided by the principles of traditional Chinese medicine, the aforementioned active components can be incorporated into modified prescriptions to support the precise treatment of NAFLD, thereby enhancing clinical efficacy.

This study utilized gene expression data from the GEO database to identify differentially expressed genes associated with NAFLD. The analysis revealed a total of 16 upregulated and 18 downregulated genes, with CXCL2, CDKN1A, TNFRSF12A, and HGFAC exhibiting the highest composite scores.

CXCL2, also known as C-X-C motif chemokine ligand 2, is a chemokine that plays a crucial role in hepatic inflammatory responses. Studies indicate that the expression level of CXCL2 is significantly elevated in the liver tissues of NAFLD patients, primarily synthesized by intrahepatic macrophages (Han et al., 2022). Ke-Qi Han et al. (Han et al., 2015) found that inhibiting the expression of CXCL1 can slow the growth of liver tumors in nude mice and suppress the expression of CXCL2, CXCL3, and interleukin-1β. Cyclin-dependent kinase inhibitor 1A (CDKN1A), a member of the Cip/Kip family, is positively correlated with the inhibition of genes related to cell cycle progression and the induction of senescence-associated genes (d'Adda di Fagagna, 2008; Kreis et al., 2019). Lamas-Paz et al. (2025) discovered through data analysis from various patient cohorts and mouse models that the expression level of CDKN1A is significantly associated with non-alcoholic steatohepatitis (NASH), liver fibrosis, and more severe liver diseases such as cirrhosis and hepatocellular carcinoma, with its overexpression exacerbating lipid metabolic disorders and inflammatory responses. Mice with a knockout of CDKN1A exhibit enhanced tolerance to metabolic liver injury, characterized by attenuated liver damage, reduced cell death, inhibited fibrotic progression, and improved lipid metabolism.

The member 12A of the tumor necrosis factor receptor family (TNFRSF12A) plays a critical role in cholestatic liver disease. Experiments have shown that hepatocyte apoptosis is closely associated with the expression of TNFRSF12A, and inhibiting TNFRSF12A expression can alleviate liver damage (Liao et al., 2023). Hepatocyte growth factor activator (HGFAC) is secreted by the liver, circulating in plasma and activated in injured tissues and tumors (Fukushima et al., 2018). The expression level of HGFAC is negatively correlated with the methylation of its promoter region, suggesting that HGFAC expression may be regulated by DNA methylation. Furthermore, the reduction in HGFAC expression is significantly associated with a shortened survival period in liver cancer patients (Yin et al., 2019).

This study suggests that Cyclosporine may represent a highly promising novel candidate compound for the treatment of NAFLD. Cyclosporine, as a calcineurin inhibitor, primarily functions to suppress inflammatory responses and alleviate symptoms related to irritation (Zhong et al., 2020). Clinically, this drug is widely utilized to prevent graft-versus-host disease (GVHD) post-transplantation and to treat various autoimmune disorders (Xie et al., 2024). In recent years, studies have indicated that the potential therapeutic applications of Cyclosporine in liver diseases have garnered extensive attention. Studies have reported (Halestrap and Davidson, 1990) that Cyclosporine can inhibit significant Ca²⁺-induced mitochondrial swelling in the liver and heart, a mechanism likely associated with its regulatory effects on intracellular calcium homeostasis. Calcium ions play a pivotal role in cellular signal transduction and metabolic regulation, and their abnormal elevation is closely associated with the pathogenesis of various liver diseases. Kobayashi et al. (2022) reported that the immunosuppressive agent Cyclosporine A (CsA) confers notable neuroprotection in models of ischemic and traumatic brain injury. This observation supports CsA’s established role in immune modulation and implies its potential hepatoprotective capacity, possibly through the regulation of apoptotic and autophagic pathways. Additionally, there is evidence to suggest that CsA could inhibit the progression of NAFLD by reducing oxidative stress within hepatocytes. Despite its documented hepatotoxicity, CsA emerged as the most promising repositioning candidate because it simultaneously targets the two hub genes CDKN1A and CXCL2 and exhibited one of the lowest docking free energies in our screen, indicating a high likelihood of direct, multi-pathway activity against NAFLD. The drug-vacant status of HGFAC, an extracellular serine protease that zymogen-activates hepatocyte-growth factor, marks this node as an attractive de novo target (Sargsyan et al., 2023). Its solvent-exposed activation loop and trypsin-fold catalytic cleft provide a tractable pocket for either covalent inhibitors or antibody blockade, opening a route for first-in-class therapy. Moreover, decades of pharmacokinetic data and recent advances in liver-targeted nano-formulations provide realistic avenues to administer CsA at micro-doses or with tissue-specific delivery, thereby minimizing systemic exposure and liver injury while preserving therapeutic efficacy. Nonetheless, the known hepatotoxic potential of CsA necessitates careful dose management in future studies. Recent metabolomic studies indicate a clear dose-dependent hepatotoxic effect of CsA: oral micro-doses (<5 mg kg^-1 day^-1) administered for 4 weeks resulted in mild, reversible biochemical changes in rats, whereas higher doses (≥10 mg kg⁻¹ day⁻¹) significantly elevated circulating bile acids and induced hepatic lipid accumulation (Yen et al., 2024). Notably, hepatic lipid accumulation exacerbates NAFLD pathology, thus conflicting with CsA’s therapeutic objectives (Losurdo et al., 2018). Consequently, future NAFLD studies employing CsA should restrict dosing to sub-immunosuppressive micro-doses, accompanied by rigorous therapeutic drug monitoring and routine assessments of bile acids, ALT, and AST, thereby maximizing metabolic benefits while minimizing hepatotoxic risk.

The clinical manifestations of NAFLD primarily include nausea, abdominal distension, loss of appetite, and fatigue. These symptoms can be categorized in traditional Chinese medicine as “intercostal pain,” “phlegm syndrome,” and “accumulation.” The pathophysiological mechanisms of this disease are primarily characterized by liver dysfunction and impaired spleen transportation. This condition often exhibits a mixed pattern of deficiency and excess, accompanied by qi and blood deficiency, with the primary lesions located in the liver, while also involving the spleen and kidneys, internal damp-heat accumulation, and obstruction of the liver meridians. Therefore, the therapeutic principles should focus on promoting liver function, strengthening the spleen, and clearing heat while transforming dampness. At the molecular level, the chemokine CXCL2 drives neutrophil infiltration and local inflammatory “heat,” mirroring the TCM concept of damp-heat accumulation in the liver; CDKN1A-induced cell-cycle arrest leads to hepatocellular stagnation, corresponding to liver dysfunction; TNFRSF12A mediates bile-acid-triggered pyroptosis and cholestatic injury, echoing “dampness obstructing bile flow;” whereas HGFAC activates hepatocyte growth factor to facilitate hepatic regeneration, functionally akin to reinforcing the spleen to restore transformation and transportation.Collectively, these gene–symptom correspondences forge a mechanistic bridge between modern molecular pathology and traditional concepts, highlighting why formulas that clear damp-heat, soothe the liver, and fortify the spleen are rational strategies for correcting NAFLD-related gene dysregulation. The two core genes and five components identified in this study align with the therapeutic principles for this disease. Notably, source tracing analysis and molecular docking visualization indicate that Polygonum cuspidatum exhibits the highest correlation with this disease. Polygonum cuspidatum topped the herb–gene network and its signature stilbenes (resveratrol, polydatin) showed the strongest docking to NAFLD hub genes CDKN1A and CXCL2, while its classical function of “clearing damp-heat and invigorating blood” mirrors the disease’s inflammatory-metabolic block. Resveratrol suppresses SREBP-1c mRNA and upregulates PPAR-α, thereby diminishing lipid accumulation and hepatocellular pathology in fatty liver disease (Mo et al., 2018). Furthermore, Resveratrol effectively protects against high-fat-induced liver injury by inhibiting endoplasmic reticulum stress and promoting autophagy (Chen et al., 2020), thus offering new avenues for the treatment of non-alcoholic fatty liver disease. The component–gene associations obtained from the HERB database in this study are merely statistical inferences based on similarities between transcriptomic signatures; such inference cannot establish a definitive regulatory relationship between the five active TCM compounds and the hub genes. Future confirmation will require ligand-binding assays, CRISPR-based perturbations, and validation at the transcriptomic or proteomic level.