An Investigation of the Effect of the Traditional Naxi Herbal Formula Against Liver Cancer through Network Pharmacology, Molecular Docking and In Vitro Experiments

1. Introduction

Cancer is the leading cause of human deaths worldwide. Moreover, cancer ranks as the first or second leading cause of death before the age of 70 in 112 out of 183 countries and ranks third or fourth in another 23 countries according to estimates from the World Health Organization (WHO) for 2019 [1,2]. Hepatocellular carcinoma (HCC) was the most common type of primary liver cancer, and it is the third leading cause of cancer-related deaths worldwide [1,3,4]. Although the exact numbers change from year to year, global cancer statistics for 2020 still rank liver cancer as the third highest cause of mortality (8.3%), with lung cancer (18%) and colorectal cancer (9.4%) being the first and second most common types of cancer for that year, respectively [1,5]. At present, there are many ways to treat liver cancer, including surgery, chemotherapy, radiotherapy, targeted therapy, and immunotherapy [6,7]. However, the use of chemotherapy and radiotherapy is leading to an increase in the prevalence of drug resistance in cancer. Moreover, immunotherapy is a limited treatment option that is only applicable to some patients. Moreover, chemotherapy and radiotherapy often produce adverse effects such as mucositis, neurotoxicity, and extravasation [7,8]. Therefore, further research is needed to develop effective cancer drugs with low toxicity and high efficiency. Traditional Chinese medicine (TCM) is often used as an adjunct therapy to mitigate the adverse reactions to radiotherapy and chemotherapy [7,9,10,11,12], and the natural products contained in TCM can be used to improve the therapeutic effect of chemotherapy, thereby improving the quality of life and prolonging patient survival [13,14,15].

The balance-regulation theory is used to regulate the integrity of the human body in TCM [16]. Specifically, traditional Chinese herbal formulas (Fang-ji, Yao-fang, or Fu-fang in Chinese) composed of several herbs are considered a cultural treasure, and have been passed down from generation to generation and used in China for thousands of years [17]. Many formulas have been reported to have effective anti-liver and anti-lung cancer effects, such as “Fei Yan Ning” (肺炎宁), traditionally used as an adjuvant therapy for lung cancer [7]. Another preparation named “Yangzheng Xiaoji Jiaonang” (养正消积胶囊) is also used to treat liver cancer [18]. Naxi medicine is one kind of TCM, which was widely used by Naxi people in Northwestern Yunnan, China. However, in many cases, Naxi medicine is prescribed without a clear understanding of the underlying mechanism. The Chong-Lou-Yao-Fang (CLYF, 重楼药方) herbal medicinal formulation is recorded in traditional Naxi medicine and its composition has been passed down through generations; CLYF has been effectively used to treat diseases and symptoms related to liver cancer for hundreds of years in the Northwestern Yunnan, China. CLYF comprises 11 herbs (Table S1), Paris polyphylla var. yunnanensis (Franch.) Hand.-Mazz. (CL), Panax bipinnatifidus Seem. (ZZS), Panax notoginseng (Burkill) F. H. Chen ex C. H. Chow (SQ), Fritillaria cirrhosa D. Don (CBM), Pleione bulbocodioides (Franch.) Rolfe (DSL), Psammosilene tunicoides W. C. Wu & C. Y. Wu (JTS), Panax quinquefolius L. (XYS), Engleromyces sinensis M.A. Whalley, Khalil, T.Z. Wei, Y.J. Yao & Whalley (ZJ), Cynanchum otophyllum C. K. Schneid. (QYS), Glycyrrhiza yunnanensis S. H. Cheng & L. K. Dai ex P. C. Li (GC), Gastrodia elata Bl. (TM). CL is an herb with potential anticancer properties in our previous research [12,19,20].

Formulas from Chinese medicine consist of multiple herbs whose active ingredients have multiple potential targets, making it challenging to study and understand the pharmacological basis of anticancer properties. The concept of “network pharmacology” was first proposed by the British pharmacologist A. L. Hopkins [21,22], based on which a four-tiered interaction network “drug-component-target-disease” was constructed to predict the mechanism of TCM and its formulas [23,24,25,26,27,28,29,30]. Network pharmacology has been instrumental in elucidating the key compounds present in the formula Le-Cao-Shi (LCS), which is active against hepatitis B, as well as the target proteins involved [16]. The main active compounds in the leaves of Eucommia ulmoides Oliver were predicted by network pharmacology for renal protection [24]. Molecular docking is used to evaluate the likelihood of interaction between a compound and its putative target protein in the human body, which involves calculating the docking energy of the compound and predicting its binding form within its target [31,32]. The interaction between the candidate compounds in LCS and the target protein in the human body was verified by a molecular docking [16]. However, underlying mechanism of studying Naxi medicine through network pharmacology and molecular docking had not been reported.

In this study, a network pharmacology strategy was used to uncover the active compounds and potential targets of traditional Naxi CLYF formula in the treatment of liver cancer. The enrichment analyses of these putative targets for Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were conducted to explore the possible mechanisms by which the CLYF formula exerts its anticancer properties. The interaction between the main compounds of the CLYF formula and their key putative therapeutic targets was validated by molecular docking simulations combined with microarray data analysis and survival analysis and confirmed their effects in vitro using a liver cancer cell line, offering a validation of the network pharmacological approach. In a summary, the purpose of this study is to validate the anticancer activity for pro-apoptosis and cytotoxic effects of traditional Naxi formula CLYF, which provides scientific basis for the application and optimization of the formula. The specific flow chart (Figure 1) is used to achieve the above objectives.

2. Results

2.1. Screening for the Active Ingredients in CLYF and Their Therapeutic Targets

The compounds for CLYF formulation were obtained through the TCMSP and BATMAN-TCM databases, as well as literature searches. The active ingredients included five chemicals in P. polyphylla var. yunnanensis (CL), three in P. bipinnatifidus (ZZS), 41 in P. notoginseng (SQ), 19 in F. cirrhosa (CBM), seven in P. bulbocodioides (DSL), one in P. tunicoides (JTS), 40 in P. quinquefolius (XYS), 12 in E. sinensis (ZJ), nine in C. otophyllum (QYS), 39 in G. yunnanensis (GC), and 14 in G. elata (TM). Eight compounds were presented in two herbs, and another three compounds were presented in three herbs. A total of 176 chemical compounds, and 1,130 related targets were found (Table S2). The Cytoscape software was employed to visualize an herb–compound–putative target network diagram, which comprised 1,314 nodes and 4,284 edges (Figure S1). Multiple compounds could be connected to a shared target; conversely, multiple targets may be connected to a unique compound in one or multiple herbs. For instance, the herb JTS only contains one active compound, psammosilenin A, but this compound is connected to 41 putative targets. These results suggest that a pure compound may have anticancer effects through multiple targets. Or several compounds may have the same anti-cancer effect. The three herbs (ZZS, SQ, XYS) from the same Panax genus may also be pharmacologically similar.

2.2. Targets for Liver Cancer

The therapeutic targets were searched by searching the GeneCards, OMIM, TTD, and DrugBank databases, using “Liver cancer” as the keyword. The search results were combined from all four databases and reduced the protein list to a non-redundant set, resulting in 4,928 putative therapeutic targets. Then the union targets of CLYF compounds and disease were mapped through the Venn diagram, 583 intersection targets of liver cancer and CLYF prescription were obtained (Figure 2).

2.3. Protein–Protein Interaction Network of the Putative Therapeutic Targets of CLYF for Treating Liver Cancer

The list of 583 proteins identified above was submitted to the String database and downloaded their putative interactors. The resulting dataset was visualized in Cytoscape: the protein–protein interaction (PPI) network comprised 880 nodes and 21,023 edges (Figure 3). Putative key therapeutic targets were identified based on their degree value ranking, representing the number of connections. The average node degree across the entire network was 47.780, with an average local clustering coefficient of 0.425, as calculated by Cytoscape. The top 20 targets were AKT1, TP53, TNF, IL6, CTNNB1, SRC, MYC, EGFR, VEGFA, JUN, MAPK3, IL1B, STAT3, HSP90AA1, CASP3, PTEN, ESR1, HIF1A, EGF, CCND1 (Table S3, including their full names), which are mainly associated with apoptosis in liver cancer cell lines.

2.4. Determining the Main Active Compounds of CLYF against Liver Cancer through Target–Pathway Network Analysis

The top 20 putative therapeutic targets were selected, then a new network was constructed in Cytoscape between them and their associated compounds extracted from each herb of the CLYF formulation (Figure 4). The active compounds were ranked according to their degree value in this new network, resulting in a list of 35 compounds. Specifically, the top active ingredients were SQ7 (quercetin), CL5 (20-hydroxyecdysone), CF3 (ginsenoside Rh2), DSL4 (blestriarene A), DSL5 (pleionesin C), JTS1 (psammosilenin A), XYS18 (phlegmariuine-N), ZJ3 (19,20-epoxycytochalasin D), ZJ6 (jiangxienone), ZJ8 (3,5,9-trihydroxyergosta-7,22-dien-6-one), ZJ11 (epoxycytochalasin D), TM12 (sucrose), QYS9 (arjunolic acid), QYS8 (glycosmisic acid), CL2 (diosgenin), ZZS1 (4'-hydroxywogonin), ZJ1 (cytochalasin D), ZJ10 (cytochalasin C), CL3 (pennogenin), DSL6 (shanciol H), ZJ4 (trichothecin), TM5 (M-hydroxybenzoic acid), CF2 (beta-sitosterol), CBM9 (verticinone), CF1 (liquiritigenin), SQ13 (sandaracopimarinol), CF10 (gamma-sitosterol), CBM16 (solanidine), XYS24 (pulegone), ZJ12 (cerevisterol), QYS3 (Baishouwubenzophenone), GC12 (2-Methyl-1,3,6-trihydroxyanthraquinone), GC16 (glycyrrhetol), GC37 (3-methyl-6,7,8-trihydropyrrolo[1,2-A] pyrimidin-2-one), ZJ2 (rosenonolactone) (Table S4). As a result, these active compounds as being the core active ingredients of CLYF with potential to treat liver cancer.

2.5. GO and KEGG Enrichment Analyses Identifies Key Signaling Pathways in Treatment of Liver Cancer by CLYF

GO and KEGG pathway enrichment analyses were performed using the top 20 putative therapeutic targets of CLYF through the Metascape database (Figure 5). As a result, 771 significantly enriched GO terms for the biological process (BP) category, 39 for the cellular component (CC) category, and 45 for the molecular function (MF) category were identified. In addition, 128 significantly enriched pathways were obtained by KEGG analysis. We extracted the top 20 terms for each of the BP, CC, and MF categories (Figure 5, Table 1).

The top BP terms mainly included ‘positive regulation of epithelial cell proliferation’, ‘positive regulation of transferase activity’, and ‘gland development’, which are closely related to the regulation of apoptosis, proliferation, phosphorylation, enzyme activity, DNA, microRNA metabolism, and transcription (Figure 5A). The top CC terms mainly involved ‘transcription regulator complex’, ‘transcription repressor complex’, ‘membrane raft’, ‘membrane microdomain’, and ‘euchromatin’ (Figure 5B). The top MF terms were mainly related to ‘DNA-binding transcription factor binding’, ‘kinase binding’, ‘transcription factor binding’, and ‘nitric oxide synthase regulator activity’ (Figure 5C).

According to the KEGG pathway enrichment analysis, the top 20 putative therapeutic targets of CLYF were associated with four main signaling pathways: ‘PI3K-Akt’, ‘mitogen-activated protein kinase’ (MAPK), ‘thyroid hormone’, and the ‘Advanced glycation end products (AGE)−Receptor of AGE (RAGE) signaling pathway in diabetic complications. Six KEGG pathways were also related to cancer: ‘cancer’, ‘breast cancer’, ‘colorectal cancer’, ‘prostate cancer’, ‘endometrial cancer’, and ‘bladder cancer’. Three other KEGG pathways were related to viral infections: ‘human cytomegalovirus infection’, ‘kaposi sarcoma-associated herpesvirus infection’, and ‘human papillomavirus infection’. In addition, the KEGG pathway ‘EGFR tyrosine kinase inhibitor resistance’ was also among the significantly enriched pathways. Notably, KEGG pathway enrichment analysis also returned an association between putative therapeutic targets related to CLYF and hepatitis B and C (Figure 5D, Table 1). Based on the above results, the anticancer properties of CLYF may involve action on multiple pathways, multiple targets, multiple biological processes, multiple molecular functions, and cell components.

2.6. Network Pharmacology Analysis through Generating a Compound–Therapeutic Target–pathway Interaction Network

In order to visualize the interactions between the compounds, their putative therapeutic target proteins, and associated pathways, a compound–target–pathway network was assembled in Cytoscape consisting of the top 20 putative therapeutic targets, the top 20 associated KEGG pathways, and the top 35 active compounds of the CLYF formulation (Figure 6). This new network comprised 75 nodes and 294 edges. The target–pathway network was also assembled in Figure S2. Complex interactions among different compounds, therapeutic targets, and pathways were observed. Among them, there were clear visual examples of the following relationships: one target–multiple pathways; one pathway–multiple targets; one target–multiple compounds; and one compound–multiple targets. (Figure 6).

Following a network topology analysis, degree value and betweenness centrality were used to screen hub targets. The core therapeutic targets in this network were MAPK3 (degree=23, betweenness centrality=0.0621), CCND1 (degree=23, betweenness centrality=0.1129), EGFR (degree=21, betweenness centrality=0.0869), AKT1 (degree=20, betweenness centrality=0.0340), TP53 (degree=17, betweenness centrality=0.0235). Of the 35 active compounds included in the network, the compound with the top degree value order was SQ7 (quercetin, degree=13, betweenness centrality=0.0378), followed by CL3 (pennogenin, degree=5, betweenness centrality=0.0033), ZJ8 (3,5,9-trihydroxyergosta-7,22-dien-6-one, degree=4, betweenness centrality=0.0151), DSL4 (blestriarene A, degree=4, betweenness centrality=0.0166), and ZZS1 (4'-hydroxywogonin, degree=3, betweenness centrality=0.0064). The core pathways related to liver cancer with the highest degree order were those related to ‘hepatitis B’ (degree=10, betweenness centrality=0.0181), ‘hepatitis C’ (degree=11, betweenness centrality=0.0191), ‘PI3K-Akt signaling pathway’ (degree=11, betweenness centrality=0.0116), ‘MAPK signaling pathway’ (degree=11, betweenness centrality=0.0224), ‘thyroid hormone signaling pathway’ (degree=9, betweenness centrality=0.0242), ‘EGFR tyrosine kinase inhibitor resistance’ (degree=9, betweenness centrality=0.0089), ‘proteoglycans in cancer’ (degree=14, betweenness centrality=0.0740), ‘pathways in cancer’ (degree=17, betweenness centrality=0.0865).

In the “compound–target–pathway” network, multiple examples of one pathway–multiple targets-multiple compounds were observed. Indeed, 10 potential therapeutic targets (AKT1, CASP3, IL6, JUN, MYC, MAPK3, SRC, STAT3, TNF, TP53) are directly related to the KEGG pathway ‘hepatitis B’ and are directly connected to 21 compounds (SQ7, ZJ11, ZJ3, ZJ4, CBM9, CF2, CF3, QYS8, QYS9, ZJ8, CL3, CL5, ZZS1, CL2, DSL4, ZJ10, ZJ6, DSL5, JTS1, TM12, XYS18, Table S2). Similarly, a set of 11 potential therapeutic targets (AKT1, CCND1, CASP3, CTNNB1, EGF, EGFR, MYC, MAPK3, STAT3, TNF, TP53) was directly related to the KEGG pathway ‘hepatitis C’ and these 11 targets were directly connected to 23 compounds (SQ7, ZJ11, ZJ3, ZJ4, CBM9, CF2, CF3, QYS8, ZJ8, CL3, CL5, ZJ6, DSL4, DSL5, JTS1, TM12, XYS18, CL2, ZJ10, ZJ12, ZJ2, DSL6, ZZS1, Table S2). In addition, 17 potential therapeutic targets (AKT1, CCND1, CASP3, CTNNB1, EGF, EGFR, ESR1, HIF1A, HSP90AA1, IL6, JUN, MYC, MAPK3, PTEN, STAT3, TP53, VEGFA) were directly related to the KEGG pathway ‘cancer’ and were directly connected to 30 compounds (SQ7, ZJ11, ZJ3, ZJ4, CBM9, CF2, CF3, QYS8, QYS9, ZJ8, CL3, CL5, ZJ6, CL2, ZJ10, ZJ12, ZJ2, DSL4, DSL5, DSL6, ZZS1, CF1, CF10, GC12, GC16, GC37, QYS3, SQ13, XYS24, CBM16, Table S2). The PI3K-Akt [33,34,35,36] and MAPK [37] signaling pathways play an anti-liver cancer role by modulating the activity of the cell cycle. In agreement, there were 11 potential therapeutic targets (AKT1, CCND1, EGF, EGFR, HSP90AA1, IL6, MYC, MAPK3, PTEN, TP53, VEGFA) related to the KEGG pathway ‘PI3K-Akt signaling’ that are connected to 15 compounds (SQ7, QYS9, ZJ8, CL3, CL5, CL2, ZJ10, ZJ12, ZJ2, DSL4, DSL5, DSL6, ZZS1, ZJ11, ZJ3, Table S2). Likewise, the MAPK signaling pathway was related to 11 potential therapeutic targets (AKT1, CASP3, EGF, EGFR, IL1B, JUN, MYC, MAPK3, TNF, TP53, VEGFA), themselves connected to 20 compounds (SQ7, ZJ11, ZJ3, ZJ4, CBM9, CF2, CF3, QYS8, ZJ8, CL3, CL5, DSL4, DSL5, JTS1, TM12, XYS18, ZJ6, DSL6, ZZS1, TM5, Table S2).

There were several examples of one target–multiple pathways-multiple compounds (Figure 6). Three therapeutic targets were selected from the top 20 core targets for characterization by molecular docking: EGFR [38,39], TP53 [40,41], and AKT1 [16,34]. These three putative therapeutic targets showed the highest association with liver cancer treatment, based on the literatures. EGFR was connected to 15 pathways: ‘cancer’, ‘proteoglycans in cancer’, ‘human cytomegalovirus infection’, ‘chemical carcinogenesis-receptor activation’, ‘colorectal cancer’, ‘breast cancer’, ‘hepatitis C’, ‘endometrial cancer’, ‘EGFR tyrosine kinase inhibitor resistance’, ‘bladder cancer’, ‘prostate cancer’, ‘MAPK signaling pathway’, ‘focal adhesion’, ‘human papillomavirus infection’, and ‘PI3K-Akt signaling pathway’. EGFR was also connected to six active compounds (CL3, DSL4, DSL5, DSL6, SQ7, ZZS1, Table S2). One compound (SQ7) and 16 pathways were connected with TP53: ‘cancer’, ‘proteoglycans in cancer,’ ‘human cytomegalovirus infection’, ‘Kaposi sarcoma-associated herpesvirus infection’, ‘colorectal cancer’, ‘breast cancer’, ‘hepatitis C’, ‘endometrial cancer’, ‘lipid and atherosclerosis’, ‘bladder cancer’, ‘hepatitis B’, ‘prostate cancer’, ‘MAPK signaling pathway’, ‘thyroid hormone signaling pathway’, ‘human papillomavirus infection’, and ‘PI3K-Akt signaling pathway’. Similarly, AKT1 was connected to one compound (SQ7) and to 19 pathways: ‘cancer’, ‘proteoglycans in cancer’, ‘human cytomegalovirus infection’, ‘Kaposi sarcoma-associated herpesvirus infection’, ‘chemical carcinogenesis-receptor activation’, ‘colorectal cancer’, ‘breast cancer’, ‘hepatitis C’, ‘AGE-RAGE signaling pathway in diabetic complications’, ‘endometrial cancer’, ‘lipid and atherosclerosis’, ‘EGFR tyrosine kinase inhibitor resistance’, ‘hepatitis B’, ‘prostate cancer’, ‘MAPK signaling pathway’, ‘focal adhesion’, ‘thyroid hormone signaling pathway’, ‘human papillomavirus infection’, and ‘PI3K-Akt signaling pathway’.

These results reflect the involvement of multiple components, multiple therapeutic targets, and multiple pathways, suggesting a complex mechanism in the action of CLYF for treating liver cancer.

2.7. Microarray data analysis and survival analysis

The expressions differences of genes between liver cancer and normal groups were compared by using top 20 targets in three GEO databases (GSE136247, GSE76427, GSE87630). All top 20 genes were found in these three GEO databases (Figure 7). There were ten genes (TP53, CTNNB1, SRC, MYC, VEGFA, MAPK3, HSP90AA1, CASP3, HIF1A and EGF) significantly upregulated and ten genes (AKT1, TNF, IL6, EGFR, JUN, IL1B, STAT3, PTEN, ESR1 and CCND1) significantly downregulated in GSE136247. In GSE76427, seven upregulated (CTNNB1, SRC, MAPK3, HSP90AA1, PTEN, HIF1A and EGF) and 13 downregulated (AKT1, TP53, TNF, IL6, MYC, EGFR, VEGFA, JUN, IL1B, STAT3, CASP3, ESR1 and CCND1) genes were observed. In GSE87630, nine upregulated (AKT1, TP53, CTNNB1, SRC, VEGFA, MAPK3, STAT3, CASP3 and EGF) and 11 downregulated (TNF, IL6, MYC, EGFR, JUN, IL1B, HSP90AA1, PTEN, ESR1, HIF1A and CCND1) genes were observed. Both EGFR and TNF were significantly downregulated in three GEO databases. Both TP53 and AKT1 were significantly downregulated in GSE76427 and upregulated in GSE87630. In GSE136247, TP53 was significantly upregulated and AKT1 was significantly downregulated. The microarray data analysis confirmed the reliability of screens results in network pharmacology and prepared for further survival analysis.

The Kaplan-Meier curve is the main expression of the survival function, showing the likelihood of an event occurring at a specific time [27]. The prognostic value of these core genes was analyzed through Kaplan-Meier survival curve (Figure 8). The results showed that EGFR, TP53 and TNF were more relevant to the prognosis of liver cancer, which may be regarded as CLYF anti-liver-cancer foci in the further research.

2.8. Molecular Docking Validation of Key Targets with Active Compounds

TCM usually studies the potential connections or effects between certain compounds in a formula and the target proteins in the human body. Active chemicals in a drug are commonly referred to as ligands, and their target proteins can be considered as their receptors. Molecular docking can be used to evaluate the probability of interaction between a ligand and a putative target protein, which can be quantified as the binding energy. If the binding energy is negative, the ligand will spontaneously bind to the target protein. If the binding energy is less than −5.0 kcal/mol, the ligand shows good binding affinity to its target protein. If the binding energy is less than −7.0 kcal/mol, reflecting stronger activity for binding configuration [22,42].

Generally, CLYF formula is orally taken through water decoction and powder in clinical application. Therefore, active ingredients in CLYF should conform the rules for ADME screening, which further predicted the drug similarity of compounds with certain biological activity. However, in the part of collecting active chemical compounds in CLYF, all the active compounds were performed by ADME screening. According to the results for survival analysis and literatures search [16,34,38,39,40,41], three core targets of EGFR, TP53, and AKT1 were selected for molecular docking tests. All three proteins were associated with the signaling pathways related to cancer (MAPK, PI3K-Akt, and hepatitis B); TP53 and AKT1 were connected to the signaling pathway related to Hepatitis C (Figure 6). SDF file for three out of the 35 most active compounds in CLYF could not been find (Table S4). Therefore, molecular docking tests for the 32 active ingredients in CLYF with available SDF files were performed with each of the three core target proteins; the binding energies are listed in Tables S5, S6, and S7. These 32 active ingredients showed strong binding to each of the three core target proteins.

A lower binding energy typically reflects a more stable binding conformation. For EGFR (Table S5), the five active compounds with the strongest binding ability were CBM16, followed in decreasing binding energy by CL2, CL3, ZJ2, and CBM9. For TP53 (Table S6), the active compounds with the strongest binding were CBM16, CL2, GC16, CBM9, and ZJ2; and for AKT1, the active compounds with the strongest binding were CBM16, ZJ2, CBM9, ZJ1, and GC16 (Table S7). Importantly, three compounds (CBM16, CBM9, and ZJ2) ranked among the top five in terms of binding ability to all three targets. The molecular docking energy of these three compounds were summarized to the three putative therapeutic targets in Table 2. All compounds were located in the target’s active pocket, most forming hydrogen bonds with the target. PyMol software was employed to visualize the molecular docking of these three active compounds and three core targets (Figure S3). Among them, EGFR had the strongest binding ability with CBM9, CBM16 and ZJ2 (Table 2), while the binding energy of these three compounds with EGFR was less than −5.0 kcal/mol. According to the ADME characteristics analysis, these three compounds conformed to the Lipinski rules, while the Gl absorption for these three compounds were high. The rotatable bonds for these three compounds (CBM9, CBM16 and ZJ2) were 0, 0 and 1, respectively. The H-bond acceptors for these three compounds (CBM9, CBM16 and ZJ2) were 4, 2 and 3, respectively. The H-bond donors for these three compounds (CBM9, CBM16 and ZJ2) were 2, 1 and 0, respectively. Therefore, the docking hydrogen bond and its length are visualized by PyMol software (Figure 9). Generally, only one hydrogen bond formed between EGFR and CBM9, named SER-912 (length of hydrogen bond=3.5). Likewise, EGFR formed one hydrogen bond with CBM16 named GLU-931 (length of hydrogen bond=2.7), while CYS-797 also formed hydrogen bond (length of hydrogen bond=1.9) around the EGFR pocket (Figure 9).

2.9. CLYF Inhibits the Proliferation of HepG2 cells

The activity of CLYF extracts on cell proliferation was assessed by using a cell counting kit in the presence of different concentrations of ethanolic, aqueous, and residual extracts. The results showed that the ethanol extract (CLYF-A) significantly inhibits the growth of tumor cells in a dose-dependent manner (Figure 10). 0.1 μM DOX is equivalent to 0.058 μg/mL. Therefore CLYF-A appeared not to be more effective than the positive control chemotherapeutic agent doxorubicin (DOX). The IC₅₀ value for the CLYF-A extract was 0.2509 mg/mL (Table 3). By contrast, the aqueous extract (CLYF-W) had no activity toward cell proliferation. These findings suggest that CLYF suppresses the proliferation of liver cancer cells by inhibiting the MAPK and PI3K-Akt signaling pathways.

CLYF can be administered in traditional Naxi medicine as a water decoction (equivalent to CLYF-W) or as a powder (equivalent to CLYF-A). In order to make consistency for the effective ingredients between the traditional edible method and modern experimental extraction method, CLYF-W was designed for the part of water decoction and CLYF-A was for powder part. The results in Figure 10 and Table 3 showed that there was no activity in the part of CLYF-W. However, CLYF-A part had inhibitory effects on HepG2 cells, which was significantly dose-dependent (Figure 10A). These indicated that the edible method of powder was better than water decoction. The top core active compounds based on the binding energy all displayed significantly inhibitory effects on the proliferation of HepG2 cell lines, with obvious dose-dependent (Figure S4). Table 3 displayed the IC₅₀ values of CLYF extracts and six core active compounds (Figure 11) for the HepG2 cell lines.

Six core active compounds in CLYF-A had been determined by HPLC (High Performance Liquid Chromatography) compared with six standards. The percentage content of the six core active compounds (CBM9, CBM16, SQ7, CL2, CL3, CF3) is 0.166 %, 0.126 %, 0.030 %, 0.006 %, 0.037 %, 0.005%, respectively. The results of these six hub active ingredients were consistent with the prediction of network pharmacology.

2.10. Effects of CLYF on the Cell Cycle and Apoptosis

The results of cell proliferation indicated that HepG2 cells may experience cell cycle arrest or other types of cell death. Flow cytometric was used to test whether inhibition of the growth of HepG2 cells also affected cell cycle progression (G0/G1, S, and G2/M). Cells were treated with different concentrations of CLYF-A (0.125, 0.25, 1.0 mg/mL) and positive control (DOX, 0.5 μM) for 24 h. The same content of DMSO was used as a negative control, which named control group in Figure 12. Compared with control group, CLYF-A treatment decreased the number of cells in G0/G1 phase while increasing cells in the G2/M phase, suggesting that CLYF-A causes G2/M cell cycle arrest, similar with the effects of DOX group (Figure 12A, Table 4). Double staining with annexin V-FITC and PI was used to determine the cell apoptosis. As shown in Figure 12B and Table 5, the cells in early-stage apoptosis (annexin V-FITC+ and PI−, Q3) were increased with the increasing concentration of CLYF-A treatment. However, the cells of medium-dose group (0.25 mg/mL) were higher than the groups of high-dose (1.0 mg/mL) and low-dose (0.125 mg/mL) group in late-stage apoptotic (annexin V-FITC+ and PI+, Q2). Comparing with the control group, the differences among three concentrations of CLYF-A were statistically significant. As shown in Figure 12, Table 4 and Table 5, CLYF-A rapidly induced both cell cycle arrest and apoptosis in HepG2 liver cancer cell lines with a dose-dependent manner.

2.11. Effects of CLYF on Mitochondrial Membrane Potential

The red fluorescence emitted by JC-1 probe indicates normal mitochondrial membrane potential, while turning to green fluorescence indicated a decrease in membrane potential, indicating that the cell is in the early stage of apoptosis. The control group showed red fluorescence (Figure 13), which used as negative control, indicating that the mitochondrial membrane potential was normal. Cells were treated with 0.5 μM DOX as a positive control, which resulted in strong green or yellow-green fluorescence, indicating that the mitochondrial membrane potential decreased, and the cells were in the early stage of apoptosis. Comparing with the negative and positive control groups, the treatment of HepG2 cells induced significant and dose-dependent decreases in mitochondrial membrane potential (Figure 13).

3. Discussion

Naxi herbal formulas have unique advantages and potential in the treatment of complex diseases, which are based on the characteristics of multi-components and multi-targets. However, the complexity of the Naxi herbal formula also limits its application and development. The Naxi formula CLYF is guided by the compatibility principles of properties, flavors, and efficacy. In the CLYF formula, the herbs XYS, GC, and QYS are used to invigorate qi and stimulate the spleen. Another five herbs are used to clear away heat and detoxify the body. The herbs CL, ZZS, SQ, CBM, DSL, JTS, and TM contribute to reducing swelling, dissipating blood stasis, relaxing the meridian, and clearing collaterals. Modern pharmacological studies have confirmed that CLYF consists of 11 herbs, ten of them having antitumor effects, and SQ showing no significant effect [12,43,44,45,46,47,48,49,50,51]. CL had different pharmacological anti-tumor, hemostatic, anti-inflammatory, analgesic, and anti-fungal activity [12,18,52]. The root and rhizome for ZZS could be used for anti-tumor, anti-inflammatory activity, analgesic sedation function, protection of cardiovascular and cerebrovascular system and protection of liver [18,46,53,54,55]. CBM was an important Chinese herbal medicine, which had the effects of relieving cough, relieving asthma, expelling phlegm, anti-inflammation, anti-tumor, lowering blood pressure, neuroprotection, analgesia and anti-oxidation [45,56]. The extracts and purified compounds from DSL had the effects of anti-cancer, hepatoprotective, anti-inflammatory, neuroprotective and antioxidant potentials [51,57,58,59]. The extracts and compounds of JTS had anti-tumor, antioxidant, analgesia and anti-rheumatoid arthritis activity, among them, carboline alkaloids had anti-tumor activity and antioxidant maltol glycosides [18,50,60,61]. XYS had a variety of pharmacological anti-tumor, antioxidative, antidiabetic, anti-stress, anti-aging, anti-fatigue, anxiolytic anti-inflammatory, neuroprotective, immunomodulatory activities [18,48,62,63,64,65]. It is used for the treatment of antibiotics, anti-inflammatory, and rheumatoid arthritis [32,40]. The crude extracts and compounds coming from in ZJ had pharmacological anti-cancer, antioxidant, and anti-inflammatory activity [49,67] while ZJ was used as antibiotics, anti-inflammatory, and rheumatoid arthritis in indigenous books [60,66]. QYS had pharmacological anti-cancer, anti-inflammatory, anti-hepatitis, anti-viral, immunomodulatory, neuroprotective effects, epilepsy, and rheumatism activity [43,60,68,69,70]. GC had pharmacological anti-cancer, anti-oxidation, anti-inflammatory, anti-virus, anti-bacterial activities and could treat the disease for respiratory, cardiovascular, memory enhancement and neuroprotective [44,60,71,72,73]. TM had different pharmacological activities including antivirus, antitumor, anti-vertigo, anti-anxiety, antidepressant, antiepileptic, anticonvulsive, anti-inflammatory, anti-aging, antipsychotic and could treat disease for sedative and hypnotic activities, neuroprotective activities, circulatory system, and cardio-cerebral-vascular [18,47,73,74,75]. According to these, CLYF has anti-cancer efficacy, but its underlying mechanism is complex with multi-components and multi-targets.

Network pharmacology analyzes the effects of drugs at a systemic level, revealing the synergistic mechanisms of drugs acting on multiple components, targets, and pathways in the human body, which is consistent with the systemic holistic view of TCM theory. At the same time, it also provides new methodological support for the transition of Naxi formula from traditional experience to theoretical science [21,22,23,24,25]. In this study, we used network pharmacology to construct herbs–compounds–targets network to systematically reveal the substance basis and underlying mechanism of CLYF for the treatment of liver cancer.

Network pharmacology revealed 176 active ingredients in the CLYF formulation and 583 potential therapeutic targets for liver cancer treatment. From these targets associated with herb-disease-target intersection, the top 20 targets were selected to assemble a PPI network, followed by GO and KEGG enrichment analyses. According to the interaction information for “compound-target” in the herb–compound–putative target network diagram (Figure S1) and degree values, one compound interacting with AKT1, TP53, CTNNB1, MYC, VEGFA, JUN, PTEN, HIF1A, and EGF, two compounds interacting with IL6 and STAT3, three compounds interacting with MAPK3, IL1B and HSP90AA1, five compounds interacting with SRC, six compounds interacting with EGFR, seven compounds interacting with CCND1, eight compounds interacting with CASP3, 12 compounds interacting with TNF, 14 compounds interacting with ESR1, were found in CLYF formula. After merging and deduplicating these compounds, there were a total of 35 compounds interacted with top 20 core targets in CLYF. Active chemical compounds in CLYF were collected and performed different screening rules based on ADME rules through TCMSP database, BATMAN-TCM database, literatures. Three core targets were EGFR, TP53, and AKT1, which were also focused on molecular docking analysis [16,34,38,39,40,41]. EGFR is an important regulator of angiogenesis [76] that is mainly involved in the maturation of neovascularization and the initiation of tumor angiogenesis, closely related to the transitional proliferation of liver cancer blood vessels [77,78]. Therefore, the predicted antitumor mechanisms of EGFR are cell cycle arrest, apoptosis promotion, anti-tumor invasion and metastasis, and anti-angiogenesis. Its activated signaling pathway can regulate tumor cell proliferation, differentiation, and survival; as well as cell cycle progression and angiogenesis [79]. TP53 is a tumor-suppressor protein with high tumor correlation; it regulates the cell cycle and induces cell apoptosis. TP53 is present at low levels in normal cells but accumulates to high levels in tumor cells [80,81]. In addition, TP53 is prone to genetic mutations in cancer cells, which can disrupt tumor-suppressor factors. Therefore, it has a regulatory effect on defective cells or the transformation from normal cells to cancer cells [39,82]. AKT is a core factor in the PI3K/AKT signaling pathway, with three members in the AKT family: AKT1/PKBα, AKT2/PKBβ, and AKT3/PKBγ. AKT1 regulates cell metabolism, apoptosis, proliferation, and glucose metabolism. AKT1 has been reported to be closely related to the differentiation, invasion, and metastasis of early tumor cells. Therefore, inhibiting AKT1 activity may inhibit tumor cell proliferation, promote tumor cell apoptosis, disrupt the energy metabolism of cancer cells and thus exert anticancer effects [83,84]. According to the results of microarray data analysis, differential expression of the top 20 core genes was found in three different datasets retrieved from three GEO databases (GSE1111, GSE1111, and GSE11111). Based on survival analysis, EGFR, TP53 and TNF were involved in the survival and prognosis of liver cancer. These targets may be specific receptors that CLYF acts on in the human body during cancer treatment. However, CLYF may affects the secretion of specific substances by acting on these specific receptors to achieve therapeutic effects.

GO and KEGG enrichment analyses revealed a significant enrichment in many cancer-related GO functions and KEGG signaling pathways. There were 128 pathways identified in the KEGG enrichment analysis, including PI3K-Akt signaling pathway, MAPK signaling pathway, hepatitis B signaling pathway, hepatitis C signaling pathway, and AGE–RAGE signaling pathway in diabetic complications. The PI3K-Akt signaling pathway modulates the cell cycle, which can be harnessed to treat cancer [36]. In addition, the PI3K-Akt pathway plays an important role in the occurrence, proliferation, induction, cell cycle progression, and apoptosis of epithelial mesenchymal transition of tumors, and contributes to drug resistance of tumors [82,84]. Importantly, the PI3K-Akt and MAPK signaling pathways have also been reported to be closely related to liver cancer treatment [33,34,35,37,85,86]. PI3K-Akt signaling pathway, MAPK signaling pathway, and AGE–RAGE signaling pathway in diabetic complications are mainly related to oxidative stress, immune regulation, and inflammatory response. Among them, the AGE-RAGE signaling pathway is closely related to inflammation, activating the MAPK pathway, interfering with immune and oxidative stress responses [37]. These signaling pathways were among the top 20 signaling pathways in the PPI network, providing a possible mechanism for treating liver cancer by CLYF. Therefore, the key to effective treatment of liver cancer with CLYF is based on these core targets and pathways. Specifically, the therapeutic effect of CLYF on liver cancer is mainly achieved by acting on specific receptors or enzymes to influence the secretion of core compounds.

This study constructed a compound-herb-key target interaction network (Figure 4) and a compound–target–pathway interaction network (Figure 6), which helped us identify the most promising active compounds in CLYF and their putative therapeutic targets. The binding between core compounds and receptors is evaluated by the binding energy of molecular docking. Specifically, this study narrowed our list of compounds to 32 core active components and three core therapeutic targets; then predicted the binding of each compound to each of these three therapeutic targets by molecular docking. These 32 active ingredients exhibited strong binding affinity toward EGFR (Table S5). Except for TM12 (sucrose) in TM, the remaining 31 core active compounds also had strong binding activity toward TP53 (Table S6) and AKT1 (Table S7). These molecular docking results are consistent with the results obtained by network pharmacological analysis.

To test the effects of the CLYF formulation and the top six active compounds, based on binding energy predictions, in vitro cell experiments were performed. CLYF-A inhibited the cell proliferation of HepG2 cells, as did each of the six individual compounds. The percentage content of these six compounds in CLYF-A had been determined by HPLC, which was consistent with the prediction for network pharmacology and molecular docking. In addition, CLYF-A also induced cell apoptosis and cell cycle arrest in HepG2 cells, which was associated with the loss of mitochondrial membrane potential. Based on this, CLYF regulates these intersected targets and pathways through these core compounds, thereby achieving therapeutic effects on liver cancer.

The main objectives in this study were the CLYF formula with clinical application in folk was reported the composition and ratio for the first time. The anticancer activity for pro-apoptotic and cytotoxic effects of CLYF was verified through scientific methods including network pharmacology, microarray data analysis, survival analysis, molecular docking and in vitro experiments, which were provided research base and methods for TCM formulas that had been clinically applied in folk. In this study, an approach of network pharmacology was used to preliminarily predict the effective ingredients and potential targets of CLYF for treating liver cancer. Further a multi-component and multi-target network was constructed through network pharmacology to better elucidate the complex interactions between active compounds and target proteins. Microarray analysis, survival analysis, and molecular docking were used as molecular verification to further validate the predictive results for targets and pathways of network pharmacology. In vitro experiments (cell proliferation assay, cell cycle and apoptosis assay, determination of mitochondrial membrane potential) combined with compounds determination were used to verify the anticancer activity for pro-apoptotic and cytotoxic effects of CLYF. However, the deeper molecular mechanism of CLYF was unclear, which need to further research through more pharmacological experiments, molecular biology experiments and in vivo experiments. Additionally, the current methods used in this study have some limitations, as the active compounds in the samples/formulas were only searched from databases or literatures, which means that the range of discovered compounds is limited to known compounds. However, there may be other undiscovered novel anti-liver cancer compounds in herbs or formulas. Using metabolomics to analyze the main active compounds in herbs or formulas, combined with pharmacology, pharmacodynamics, pharmacokinetics and other methods, to screen for more effective components and analyze of their molecular mechanisms, will be more convincing.

4. Materials and Methods

4.1. Medicinal Materials, Chemicals, Reagents, Cell Culture and Treatment

The herbs that comprise CLYF were provided by Lijiang Yunxin Green Biological Development Co., Ltd., in Northwestern Yunnan Province, China, in May 2022. The 11 herbs were stored in Key Laboratory of Economic Plants and Biotechnology, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, China and identified in May 2022 according to their morphological features by Professor Heng Li (CAS Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, China). The voucher number for each herb was showed in Table S1. Ethanol was analytical grade. A centrifuge (L530, Changsha High and New Technology Industrial Development Zone, Xiangyi Centrifuge Instrument Co., Ltd., China) was used to separate solids from liquids in suspension. Quercetin (SQ7) was purchased from ShangHai D&B Biological Science and Technology Co. Ltd. (purity ≥97%). Solanidine (CBM16) was purchased from Beijing Wokai Biotechnology Co., Ltd. (purity ≥99%). Diosgenin (CL2), pennogenin (CL3), ginsenoside Rh₂ (CF3), and verticinone (CBM9) were purchased from Sichuan Weikeqi Biological Technology Co., Ltd. (purity ≥98%).

Human liver cancer cell lines (HepG2) were purchased from BNCC (BeNa culture Collection, Beijing, China) and cultured in DMEM medium (Procell Life Science & Technology Co., Ltd.) supplemented with 10% (v/v) fetal bovine serum (Procell Life Science & Technology Co., Ltd.) and 1% (v/v) penicillin-streptomycin solution (stocks of 10,000 U/mL penicillin G sodium salt; 10 mg/mL streptomycin, Procell Life Science & Technology Co., Ltd.) in an incubator set to 37 ℃ and 5% (v/v) CO₂.

4.2. CLYF Preparation and Extraction

The formula CLYF was the clinical experience of Naxi people in the treatment of liver cancer with significant therapeutic effects, which was recorded in the Naxi language in local medical books, passed on at least three generations. In order to prepare the powder that is the basis of the CLYF formulation, the dried parts of 11 herbs were crushed into powder with a pulverizer. The relative proportion of each herb is listed in Table S1. The extraction was performed based on traditional methods of CLYF consumption in Naxi by using ethanol and deionized water as solvents [12]. The CLYF powder (50.0 g) was first extracted three times in 75% (v/v) ethanol (400 mL each extraction) and refluxed at 100 ºC for 4 h. The filtrate was filtered through filter paper, merged, and concentrated by rotary evaporator as ethanol extract (CLYF-A, 12.1 g). The residual ethanol was allowed to evaporate naturally from the filtered residue, which was then extracted with deionized water (H₂O, 3 × 400 mL) by refluxing at 100 ºC for 4 h. The liquid portion was collected by centrifugation (4,000 r/min, 10 min, 25 ºC), and combined, concentrated by rotary evaporator as the aqueous CLYF extract (CLYF-W, 8.9 g). The drugs after filtration were dried and taken as the residual portion of CLYF (CLYF-S, 26.0 g) for further research.

4.3. Collection of Active Chemical Compounds and Their Corresponding Targets in CLYF

CLYF is composed of 11 herbs. The active chemical compounds in each herb were obtained from the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) (https://tcmspw.com/tcmsp.php (accessed on 4 May 2023)), the Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine (BATMAN-TCM) (http://bionet.ncpsb.org.cn/batman-tcm/index.php/ (accessed on 4 May 2023)), and literatures searches from Web of Science (https://www.webofscience.com/wos/alldb/basic-search (accessed on 4 May 2023)). In TCMSP, the ADME property value was set to Oral Bioavailability (OB) ≥ 30% and Drug Likeness (DL) ≥ 0.18. In BATMAN-TCM, the input parameter score cutoff was set to 20, and the adjusted P-value was 0.05. The chemical structures and SMILES (Simplified Molecular Input Line Entry System) code of these active compounds were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/ (accessed on 6 May 2023)). The Swiss Institute of Bioinformatics (SIB) database was used to perform ADME screening (http://www.swissadme.ch/ (accessed on 7 May 2023)). The ADME characteristics were Lipinski (yes), Gl absorption (high), the rotatable bonds (less than or equal to 5), the H-bond acceptors (less than 10), the H-bond donors (less than 5). The compound targets were predicted with the Swiss Target Prediction tool (http://www.swisstargetprediction.ch/ (accessed on 7 May 2023)). The official gene names of all target genes were confirmed or corrected by querying the UniProt database (https://www.uniprot.org/ (accessed on 8 May 2023)).

4.4. Collection of Disease-Related Proteins and Therapeutic Targets

Therapeutic targets for treating liver cancer were collected from four databases: GeneCards (https://www.genecards.org/ (accessed on 12 May 2023)), Online Mendelian Inheritance in Man (OMIM, https://omim.org/ (accessed on 12 May 2023)), Therapeutic Target (TTD, http://db.idrblab.net/ttd/ (accessed on 12 May 2023)), and DrugBank (https://go.drugbank.com/ (accessed on 12 May 2023)). The repetitive region of gene targets was removed. The overlap of therapeutic targets across the four databases was visualized by drawing a Venn diagram. The official names of all therapeutic targets were confirmed or corrected by querying the UniProt database. Then, the overlap between the putative targets for each active ingredient and the therapeutic targets was determined using the online tool Venny (https://bioinfogp.cnb.csic.es/tools/venny/ (accessed on 19 May 2023)). The intersection between the two sets of proteins were considered as potential targets for CLYF formula in anticancer treatment.

4.5. Network Construction and Analysis

The putative CLYF targets for cancer treatment, as defined above, were used as input in the String database (https://string-db.org/ (accessed on 24 May 2023)) to obtain a list of interacting proteins and estimate the strength of the putative interaction. The resulting file was then loaded into Cytoscape 3.7.1 software to define the underlying PPI network. The Network Analyzer tool within Cytoscape was used to analyze the node degree, and the top 20 therapeutic targets and their associated core active ingredients of CLYF with anticancer properties were extracted. Cytoscape software was used to construct the medicine–component–target diagrams.

4.6. GO and KEGG Pathway Enrichment Analyses

GO and KEGG pathway enrichment analyses were carried out on core therapeutic targets by the database Metascape (https://metascape.org (accessed on 14 June 2023)), with the species set to “Homo sapiens” and the P-value cutoff set to 0.01 [87]. The GO categories Biological Process (BP), Cellular Component (CC), and Molecular Function (MF), and KEGG pathways were selected for enrichment analysis. The top 20 GO terms obtained from the BP, CC, and MF categories, and for KEGG pathways were selected based on their P-values and the results were visualized on the website https://www.bioinformatics.com.cn (accessed on 16 June 2023).

4.7. Microarray data analysis

Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) is a high-throughput gene expression database that includes microarray chips, second-generation sequencing data, and other forms of high-throughput genomic data [27]. Three microarray datasets GSE136247, GSE76427 and GSE87630 were downloaded from GEO database (accessed on 2 August 2024). Differentially expressed genes were screened through the “Limma” package of R language (R version 4.4.0). Differential expression information was obtained by adding annotations and deduplicating through the “dplyr” package of R language. The “ggplot2” package of R was used to construct volcano maps and visualize changes in gene expression. P<0.05 was considered as statistically significant. The microarray data analysis of top 20 targets was performed to screen out genes that were significantly upregulated or downregulated in liver cancer compared with the normal group for later survival analysis.

4.8. Survival analysis

The online tool for Kaplan-Meier Plotter（https://kmplot.com/analysis/ (accessed on 3 August 2024)）was used to survival analysis of hub genes and evaluation of the clinical importance of a particular gene based on gene expression levels. The type for disease was selected as liver cancer and the organism was set to “Homo sapiens”. Kaplan-Meier survival curves were constructed for survival analysis and examining the relationship between core genes and overall survival of liver cancer patients. P<0.05 was considered as statistically significant. Whether high or low expression of the core genes was beneficial for liver cancer patients was analyzed through the changes in gene expression levels during clinical survival.

4.9. Molecular Docking between CLYF Active Compounds and Their Target Proteins

Molecular docking was performed to evaluate the interaction between target proteins and the active compounds in CLYF [22]. The structure data (SDF) files of the 3D structure of key active components present in CLYF were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/ (accessed on 28 June 2023)). Then, all SDF files were transformed to MOL2 file format through the Open Babel 3.1.1 software (https://openbabel.org/wiki/Main_Page (accessed on 28 June 2023)). Lastly, the active compounds were pre-processed with full hydrogenation, set up as ligands, and automatically assigned a charge. The detection and settings of rotatable bonds were completed and saved as PDBQT format for storage and further use.

The protein structures of epidermal growth factor receptor (EGFR) [88], cellular tumor antigen p53 (TP53) [89], and RAC-alpha serine/threonine-protein kinase (AKT1) [90] were downloaded from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) (https://www.rcsb.org/ (accessed on 6 July 2023)). The organism was set to “Homo sapiens” and the resolution was set to less than 2.5 Å. The original water molecules and small molecular ligands were removed by PyMol 2.4.0 software and saved as PDB format. In order to completing the selection of receptors, hydrogens and total Gasteiger charges were added, and then saved as PDBQT format for molecular docking analysis. The Autogrid was used to set the grid box to cover the entire molecular to identify the docking active center. After defining the binding site of targets, the first docking calculation was set through the Autodock 4.2.6 software. The method of Local Search Parameters was used to calculate the second docking through Autodock. Last, all binding energies were quantified to verify the predicted results of network pharmacology. The smaller the binding energy was, the better the binding activity between the target protein and core active compounds. The final files were saved as PDBQT format for further visualization by PyMol software.

4.10. Determination of Six Core Active Compounds Content

Samples were dissolved in methanol and filtered through a 0.45 μm membrane before use for HPLC analysis. Six core active compounds (SQ7, CL2, CL3, CF3, CBM9, CBM16) were divided into two groups to analyze and determine the content in the extract for CLYF-A according to compound property by HPLC (Agilent 1260 series system, Zorbax SB-C18 column, 5 μm, ϕ 4.6 × 250 mm). Four compounds (SQ7, CL2, CL3, CF3) were analyzed by HPLC-DAD (High Performance Liquid Chromatography-Diode Array Detection) and prepared at the concentration of 0.2695, 0.2719, 0.2681, 0.2703 mg/mL. Two compounds (CBM9 and CBM16) were analyzed by HPLC-ELSD (High Performance Liquid Chromatography-Evaporative Light Scattering Detector) and prepared at the concentration of 0.5401 and 0.5389 mg/mL.

The mobile phase for SQ7 consisted of (B1) methanol and (A1) 0.3% formic acid in water (v/v), with a flow rate of 1 mL/min, at 30°C, with injection volume in 10 μL, with the following gradient: 0–10 min, 10%–20% B1; 10–25 min, 20%–25% B1; 25–30 min, 25% B1; 30–35 min, 25%–50% B1; 35–40 min, 50% B1; 40–45 min, 50%–80% B1; 45–48 min, 80%–0% B1. The mobile phase for three compounds (CL2, CL3, and CF3) consisted of (B2) acetonitrile and (A2) water, with a flow rate of 1 mL/min, at 30°C, with injection volume in 10 μL, with the following gradient: 0–20 min, 20% B2; 20–40 min, 20%–95% B2; 40–65 min, 95% B2. The mobile phase for CBM9 consisted of (B2) acetonitrile and (A3) 0.03% triethylamine in water (v/v), with a flow rate of 1 mL/min, at 30°C, with injection volume in 2 μL, with the following gradient: 0–40 min, 70% B2. Both the evaporator and nebulizer temperature of ELSD were 40°C, and the gas flow rate was 1.4 L/min. The mobile phase for CBM16 consisted of (B2) acetonitrile and (A4) 0.1% trifluoroacetic acid in water (v/v), with a flow rate of 1 mL/min, at 35°C, with injection volume in 10 μL, with the following gradient: 0–10 min, 5%–100% B2; 10–13 min, 100% B2. Both the evaporator and nebulizer temperature of ELSD were 75°C, and the gas flow rate was 2.5 L/min.

4.11. Cell Proliferation Assay

The effects of the CLYF extracts on cell proliferation were evaluated by the CCK-8 (Cell Counting Kit-8) assay on human liver cancer cell lines [91,92,93,94,95]. Each well of a 96-well plate (Corning Life Sciences [Wujiang] Co., Ltd.) was seeded with 90 μL cell culture medium (adherent cells 5×10⁴/mL and suspension cells 9×10⁴/mL) per well, and then cultured at 37℃ and 5% (v/v) CO₂ for 24 h. The cell suspension contained 10% (v/v) fetal bovine serum (Procell Life Science & Technology Co., Ltd.). CLYF extracts (10 μL) were added to each well at the following concentrations, each in three replicates: CLYF-A (0.1, 0.5, 1.0, 5.0, 10.0, 20 mg/mL) and doxorubicin (DOX, 0.001, 0.01, 0.1, 1.0, 5.0, 10.0 μM, Solarbio) were dissolved in dimethyl sulfoxide (DMSO, Sigma); CLYF-W (0.1, 0.2, 0.5, 1.0, 2.0, 5.0 mg/mL) was dissolved in sterile water. The following six compounds were dissolved in DMSO: CBM9 (1.0, 10.0, 50.0, 100.0, 200.0, 500.0 μM), CBM16 (0.01, 0.1, 10.0, 50.0, 200.0, 500.0 μM), SQ7 (0.1, 1.0, 50.0, 100.0, 200.0, 500.0 μM), CL2 and CL3 (0.1, 10.0, 50.0, 100.0, 200.0, 500.0 μM), CF3 (0.01, 1.0, 50.0, 100.0, 200.0, 500.0 μM). 1.0 μM compound is equivalent to 1.0 × M × 10^-6 mg/mL. M is the molecular weight of the compound in g/mol. After culturing for 48 h in an incubator at 37℃ and 5% (v/v) CO₂, old culture medium and drug solution was removed from adherent cells and 100 μL CCK-8 solution (diluted 10 times with alkaline medium) was added. For suspension cells, 10 μL CCK-8 stock solution was directly added to each well. After an incubation of 1–4 h in the dark and real-time observation, the absorbance at 450 nm was measured by a microplate reader (Multiskan MK3, Thermo). The IC₅₀ value was calculated with Graphpad Prism 8 (version 8.0.2, from GraphPad Software Inc). DOX was used as a positive control in the cell experiments.

4.12. Cell Cycle and Apoptosis Assays

Flow cytometry (BeckmanCytoFLEX, Beckman, United States) was used to measure the effects of CLYF on cell cycle and apoptosis by staining with PI (propidium iodide, Beyotime Biotech. Inc., Shanghai, China) and annexin V-FITC/PI (annexin V-fluorescein isothiocyanate/ propidium iodide, Beyotime Biotech. Inc., Shanghai, China), respectively. HepG2 cells were cultivated in six-well plates at a density of 10×10⁵ cells per well for 24 h before the treatment with different concentrations of CLYF-A (0.125, 0.25, 1.0 mg/mL) and positive control (0.5 μM) for 24 h. Doxorubicin (DOX, Solarbio) was used as the positive control, which was the chemotherapy drug in cancer treatment. Adherent cells were treated with trypsin digestion fluid (Cell Cycle and Apoptosis Analysis Kit, Beyotime Biotech. Inc., Shanghai, China). After then, all cells were collected, centrifuged, and washed with phosphate-buffered saline (PBS). For the determination of cell cycle, cells were fixed in 70% (v/v) ice-cold ethanol at 4 ℃ for 2 h, and then centrifuged for 5 min to precipitating cells, washed again in pre-ice-cold PBS. The cells were stained by 0.5 mL PI for each tube at 37 ℃ for 30 min. The cell cycle stage was analyzed by flow cytometry. Red fluorescence at 488nm excitation wavelength and light scattering were detected. For the determination of cell apoptosis, the cells were incubated with 195 μL binding buffer of annexin V-FITC, after which 5 μL annexin V-FITC and 10 μL PI were added at 37 ℃ for 20 min before determination by flow cytometry.

4.13. Determination of Mitochondrial Membrane Potential

24-well white transparent bottom cell plates were seeded with 1×10⁵ cells per well and cultured for 24 h. The cultured medium was removed. And then the cells were treated with different concentrations of CLYF-A (0.125, 0.25, 0.5 mg/mL) and positive control (DOX, 0.5 μM) for 24 h. After removing the cultured medium, the cells were added 0.5 mL JC-1 dyeing solution (Annexin V-FITC Apoptosis Detection Kit, Beyotime Biotech. Inc., Shanghai, China) at 37 ℃ for 20 min. The supernatant was removed after incubation, and the cells were washed twice with 1 mL JC-1 buffer (×100). Cell culture medium (0.5 mL) was added to each well, and the stained cells were observed with inverted fluorescence microscope (IX70, Olympus, Japan). The wavelengths for excitation were 490 nm (monomers for JC-1) and 525 nm (polymer for JC-1). And the wavelengths for emission were 530 nm (monomers for JC-1) and 590 nm (polymer for JC-1).

4.14. Statistical Analysis

All data are presented as means ± standard deviation (SD) from three parallel experiments. One-way ANOVA was used to evaluate significant differences followed by LSD, Tukey’ s-b, and Waller-Duncan tests by SPSS software (Version 26, IBM). P<0.05 was considered as statistically significant.

5. Conclusion

In summary, the anticancer activity for pro-apoptotic and cytotoxic effects for Naxi formula CLYF was verified through scientific methods including network pharmacology, microarray data analysis, survival analysis, molecular docking and in vitro experiments in this study. Firstly, the approach of network pharmacology was to predict the effective ingredients and potential targets (mainly were EGFR, TP53, and AKT1) of CLYF against liver cancer and construct a multi-component and multi-target network that better elucidated the complex interactions between active compounds and target proteins. Secondly, CLYF tend to treat liver cancer through PI3K-Akt signaling pathway, MAPK signaling pathway, hepatitis B signaling pathway, and hepatitis C signaling pathway to exert drug effects. The help of microarray data analysis, survival analysis, and molecular docking were further used as molecular verification for the predictive targets and pathways of network pharmacology. In addition, cell experiments indicated that CLYF-A suppressed cell proliferation and induced cell apoptosis, cell cycle arrest in HepG2 cells, which was associated with the loss of mitochondrial membrane potential. In conclusion, this study systematically verified the antitumor efficacy of the CLYF formulation, which is based on Naxi indigenous knowledge. This study validated the key components, target therapeutic proteins, and signaling pathways through network pharmacology, microarray data analysis, survival analysis, molecular docking, and in vitro experiments, providing new possibilities for clinical treatments and CLYF formulations. This study also provides methods and case support for the scientific verification of the efficacy of many ethnic folk formulas. Of course, there are some limitations to this study: firstly, the deeper molecular mechanisms underlying the action of CLYF need to more pharmacological experiments, molecular biology experiments. Secondly, this study did not test the efficacy or safety of CLYF in vivo. Finally, more antitumor cases, toxic effects, and adverse events should be analyzed in clinical studies.