Submit an Article
Online ISSN: 1945-4589
Impact Journals, LLC. is the publisher of Aging: www.impactjournals.com.
Impact Journals is a member of the Wellcome Trust List of Compliant Publishers.
Impact Journals is a member of the Society for Scholarly Publishing.
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 14, Issue 21 pp 8783—8804

Multilevel data integration and molecular docking approach to systematically elucidate the underlying pharmacological mechanisms of Er-Zhi-Wan against hepatocellular carcinoma

Shaoyan Zheng1,2, , Botao Pan1, ,

  • 1 Affiliated Foshan Maternity and Child Healthcare Hospital, Southern Medical University, Foshan 528000, P.R. China
  • 2 Traditional Chinese Medicine Department, Affiliated Foshan Maternity and Child Healthcare Hospital, Southern Medical University, Foshan 528000, P.R. China

Received: August 17, 2022       Accepted: October 27, 2022       Published: November 7, 2022      

https://doi.org/10.18632/aging.204369
How to Cite

Copyright: © 2022 Zheng and Pan. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

As a multicomponent, multitarget empirical therapy, traditional Chinese medicine (TCM) has been used clinically in Asia for thousands of years. Due to this unique feature, TCM therapy is considered a promising therapeutic strategy for the treatment of hepatocellular carcinoma (HCC). Er-Zhi-Wan (EZW), a well-known TCM formula containing two herbs, Fructus Ligustri Lucidi (FLL, Nü-Zhen-Zi) and Ecliptae Herba (EH, Mo-Han-Lian), is commonly used in clinical practice to prevent and treat liver diseases. Modern pharmacological studies have shown that both EH and FLL can inhibit HCC proliferation. However, the pharmacological mechanism, potential targets, and clinical value of EZW in inhibiting HCC have not been fully elucidated. We used multilevel databases (Gene Expression Omnibus (GEO), Traditional Chinese Medicine Systems Pharmacology (TCMSP), High-throughput Experiment- and Reference-guided database (HERB), and SwissTargetPrediction) to show that EZW suppresses HCC through 19 active components acting on 66 potential targets. Enrichment analysis revealed that EZW mainly regulates HCC progression through various metabolic pathways, the cell cycle, and cellular senescence. Furthermore, we used The Cancer Genome Atlas (TCGA)-LIHC database to analyze the expression patterns and clinical characteristics of cellular senescence-related genes and identified CDK1, CDK4, CHEK1, and G6PD as key therapeutic molecular targets in EZW-suppressed HCC. Molecular docking revealed that EZW could exert its anti-HCC effect by binding various active components to the above cellular senescence-related genes and regulating their activities. In conclusion, we systematically revealed the potential pharmacological mechanisms and molecular targets of EZW against HCC based on multilevel data integration and a molecular docking strategy.

Introduction

Worldwide, liver cancer is one of the most common cancers, and its morbidity and mortality are on the rise [1, 2]. It has a mortality-to-morbidity ratio of 0.91 and is 2.3 times more common in men than in women [3]. Notably, the outlook for patients is even grimmer in Asia, where 72% of new cases are reported to be diagnosed (over 50% in China), with five-year survival rates as low as 12% [4, 5]. Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer and often develops from chronic liver disease caused by hepatitis B virus or hepatitis C virus infection, alcoholism, or metabolic syndrome [6]. Currently, surgical resection, liver transplantation, and locoregional therapy (including radiofrequency ablation) are recommended as curative treatments for only one-third of HCC patients. The remaining 60%–70% of patients receive noncurative treatments, such as molecularly targeted agents, monoclonal antibodies, or immune checkpoint inhibitors, as initial therapy [7]. However, due to the heterogeneity and complexity of HCC, most patients are diagnosed at an advanced stage; therefore, systemic therapy is often recommended as the standard of care [8]. Although it is more effective than monotherapy, it is only suitable for use in a small number of patients and is associated with severe toxicity [9].

The pathological mechanism of HCC is very complex, and multiple targets and signaling pathways are involved in its development [10]. This complexity requires the development of a therapeutic strategy that modulates multiple targets in HCC to improve patient outcomes. As a multicomponent and multitarget empirical therapy, traditional Chinese medicine (TCM) has been recognized worldwide in recent years for its multitarget synergistic intervention effect on HCC. Clinically, as an adjuvant drug, TCM has shown good efficacy in HCC patients, significantly reducing the incidence, preventing the recurrence, and improving the overall survival of HCC patients [11–14]. In addition, modern pharmacological studies have shown that TCM or TCM-derived natural medicines can effectively suppress the development and progression of HCC in vitro and in vivo [15–18]. Due to this unique feature, TCM is considered a promising therapeutic strategy for the treatment of complex diseases, including liver cancer.

Er-Zhi-Wan (EZW), a well-known TCM formula, was first recorded in the ‘Fu Shou Jing Fang’ in the Ming Dynasty [19]. It consists of an equal weight mixture of two herbs, Fructus Ligustri Lucidi (FLL, Nü-Zhen-Zi) and Ecliptae Herba (EH, Mo-Han-Lian). According to Chinese medicinal theory, EZW is commonly used clinically in China for the treatment of liver-kidney yin deficiency syndrome (LKYDS), which is a pathological and diagnostic pattern caused by an imbalance of yin and yang and is more common in primary liver cancer, diabetes, and hypertension [20]. Therefore, EZW is commonly used clinically in China to prevent or treat various liver and kidney diseases. Yao et al. [21] uncovered the hepatoprotective effect of EZW based on a metabolomic strategy. Hu et al. [22] revealed that FLL extract induced apoptosis and cellular senescence in human hepatoma cells by upregulating p21, confirming that FLL is a potential anticancer herb for the treatment of HCC. Moreover, our previous study showed that EH extract could inhibit the proliferation of HCC cells by inhibiting PI3K-AKT signaling [18]. However, few studies have comprehensively investigated the molecular mechanisms involved in EZW in the treatment of HCC.

Since TCM prescriptions comprise many kinds of herbs and contain many kinds of ingredients, it is difficult to systematically and comprehensively study the pharmacological mechanism of TCM with the existing experimental methods. With the successful establishment of multiple biological databases and the rapid development of systems biology, the emergence of network pharmacology has brought great opportunities for breakthroughs in TCM research [23, 24]. To date, this method has been successfully used to elucidate the multitarget efficacy of TCM in the treatment of various diseases, effectively bridging the gap between Western medicine and Chinese medicine [25, 26]. In this study, we employed various biological databases and biocomputational approaches to investigate the pharmacological network involved in EZW in the treatment of HCC to predict potential molecular targets and pharmacological mechanisms. The overall research flowchart is shown in Figure 1.

Flowchart of the analytical procedures of the study.

Figure 1. Flowchart of the analytical procedures of the study.

Results

Identification of pathological genes in HCC

To identify which genes are involved in the progression of HCC, we analyzed the Gene Expression Omnibus (GEO) dataset GSE84402 to identify genes that were differentially expressed in 14 pairs of HCC tissues and corresponding noncancerous tissues. As shown in the volcano plot, a total of 1199 differentially expressed genes (DEGs) were identified in these 14 pairs of liver tissues, of which 632 genes were upregulated and 567 genes were downregulated in cancerous tissues compared with noncancerous liver tissues (Figure 2A and Supplementary File 1). It is speculated that the progression of HCC involves extensive and complex pathological gene regulation.

Identification and enrichment analysis of differentially expressed genes in 14 pairs of HCC tissues and corresponding noncancerous tissues. (A) The expression patterns of the DEGs are shown in volcano plots. Red and blue points represent upregulated genes (log2FC ≥ 1) and downregulated genes (log2FC ≤ -1), respectively, while gray represents genes with no significant difference in expression (P.adj B, C) Bubble plot showing the top 20 GO (B) and KEGG (C) enrichment analysis results. The larger the ordinate value in the bubble chart, the more significant the corresponding GO or KEGG result is. The abscissa represents the normalized upregulation and downregulation value (the ratio of the difference between the number of upregulated genes and the number of downregulated genes to the total number of differentially expressed genes). The higher the value is, the higher the number of upregulated genes enriched in the GO/KEGG pathway results; conversely, the lower the value is, the higher the number of downregulated genes enriched in the GO/KEGG pathway results. (D, E) Heatmaps showing the expression patterns of genes involved in the cell cycle (ko04110) or metabolic pathways (ko01100).

Figure 2. Identification and enrichment analysis of differentially expressed genes in 14 pairs of HCC tissues and corresponding noncancerous tissues. (A) The expression patterns of the DEGs are shown in volcano plots. Red and blue points represent upregulated genes (log2FC ≥ 1) and downregulated genes (log2FC ≤ -1), respectively, while gray represents genes with no significant difference in expression (P.adj < 0.05). (B, C) Bubble plot showing the top 20 GO (B) and KEGG (C) enrichment analysis results. The larger the ordinate value in the bubble chart, the more significant the corresponding GO or KEGG result is. The abscissa represents the normalized upregulation and downregulation value (the ratio of the difference between the number of upregulated genes and the number of downregulated genes to the total number of differentially expressed genes). The higher the value is, the higher the number of upregulated genes enriched in the GO/KEGG pathway results; conversely, the lower the value is, the higher the number of downregulated genes enriched in the GO/KEGG pathway results. (D, E) Heatmaps showing the expression patterns of genes involved in the cell cycle (ko04110) or metabolic pathways (ko01100).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs in HCC

To further investigate the underlying molecular mechanisms involved in HCC progression, we performed enrichment analysis on 1199 DEGs, including GO and KEGG analyses. We used bubble plots to display the top 20 GO (Figure 2B) and KEGG (Figure 2C) enrichment analysis results. GO functions can be divided into three categories: biological process (BP), molecular function (MF), and cellular component (CC), and the top 20 GO annotation results show that these DEGs are mainly enriched in the BP category. These BPs mainly involve metabolic and cellular processes, including the following: “carboxylic acid metabolic process (GO:0019752)”, “organic acid metabolic process (GO:0006082)”, “oxoacid metabolic process (GO:0043436)”, “mitotic cell cycle (GO:0000278)”, “mitotic cell cycle process (GO:1903047)”, “cell cycle (GO:0007049)”, and “cell division (GO:0051301)”.

In addition, the results of the top 20 KEGG pathways indicated that the progression of HCC is closely related to 6 categories, including metabolism, cellular processes, human diseases, genetic information processing, organismal systems, and environmental information processing. In particular, 12 of the 20 pathways were involved in the metabolism category, including “Metabolic pathways (ko01100)”, “Fatty acid degradation (ko00071)”, “Retinol metabolism (ko00830)”, “Tryptophan metabolism (ko00380)”, “Tyrosine metabolism (ko00350)”, “Alanine, aspartate and glutamate metabolism (ko00250)", "Valine, leucine and isoleucine degradation (ko00280)”, “Drug metabolism – cytochrome P450 (ko00982)”, “Caffeine metabolism (ko00232)”, “Metabolism of xenobiotics by cytochrome P450 (ko00980)”, “Carbon metabolism (ko01200)”, and “Biosynthesis of amino acids (ko01230)”. Notably, metabolic pathways (ko01100) were enriched with multiple downregulated DEGs, including multiple metabolic enzymes (PFKFB1, AKR1D1, CYP1A2, XDH, ALDH8A1, LDHD, GCDH, ADH6, ADH4, ABAT, and HOGA1). The heatmap results showed that abnormal changes in various metabolic enzymes occurred during the progression of HCC, which suggested that abnormal expression of metabolic enzymes affected the development of HCC (Figure 2D). Furthermore, HCC progression is regulated by the cell cycle (ko04110), and multiple upregulated DEGs were found to be significantly enriched in this pathway, including multiple cell cycle-related regulatory genes (CCNB3, PCNA, CDK1, CDK4, CCNB1, CCNB2, and CHEK1). The combined heatmap analysis indicated that the cell cycle is abnormally activated in HCC, leading to the indefinite growth of HCC (Figure 2E).

Screening of active compounds and potential therapeutic targets of EZW

EZW is produced by mixing Fructus Ligustri Lucidi and Ecliptae Herba in equal proportions. According to the two criteria of drug-likeness (DL) ≥ 0.18 and oral bioavailability (OB) ≥ 30% [27, 28], a total of 9 active ingredients in EH and 10 active ingredients in FLL were identified in TCMSP. Among them, studies have shown that both oleanolic acid and ursolic acid are active. components of FLL, but they were not included in the TCMSP screening analysis because they did not meet the DL and OB conditions [29]. Therefore, to increase the credibility of the study, these two active ingredients were included in the list of active ingredients of FLL By combining the active ingredients of EH and FLL, it was revealed that EZW comprised a total of 19 active ingredients, of which the two herbs each contained quercetin and luteolin (Supplementary File 2). The chemical structures of these active compounds are shown in Figure 3.

Chemical structures of 19 active ingredients of EZW.

Figure 3. Chemical structures of 19 active ingredients of EZW.

We identified potential targets for these 19 active compounds through two target prediction databases, High throughput Experiment- and Reference-guided (HERB) and SwissTargetPrediction. In the HERB database, flavonoids (luteolin, quercetin, and kaempferol) possessed more targets, but targets failed to be identified for 4 active compounds (linarin, lucidumoside D, lucidusculine, and olitoriside) (Figure 4A and Supplementary File 3). Similarly, in the SwissTargetPrediction database, flavonoids (luteolin, quercetin, and kaempferol) also possessed more targets, while targets for taxifolin, lucidusculine, and olitoriside failed to be identified (Figure 4B and Supplementary File 3). After removing duplicate values, 215 potential targets of 9 active compounds in EH were identified from the two databases. Likewise, 379 potential targets of 12 active compounds in FLL were identified. Finally, 446 potential targets of 19 active ingredients in EZW were identified.

Prediction and screening of potential targets of EZW for the treatment of HCC. (A, B) Prediction and collection of potential targets for EZW based on the HERB and SwissTargetPrediction databases. (C) Venn diagram identifying 66 potential therapeutic targets for EZW in the treatment of HCC (35 targets were downregulated and 31 were upregulated in HCC). (D) Heatmap analysis of the expression patterns of 66 potential therapeutic targets of EZW in the treatment of HCC in the GSE84402 dataset.

Figure 4. Prediction and screening of potential targets of EZW for the treatment of HCC. (A, B) Prediction and collection of potential targets for EZW based on the HERB and SwissTargetPrediction databases. (C) Venn diagram identifying 66 potential therapeutic targets for EZW in the treatment of HCC (35 targets were downregulated and 31 were upregulated in HCC). (D) Heatmap analysis of the expression patterns of 66 potential therapeutic targets of EZW in the treatment of HCC in the GSE84402 dataset.

Target screening, network, and topological analysis of EZW in the treatment of HCC

To explore which pathological targets EZW acts on to treat HCC, we performed Venn diagram analysis of 1199 DEGs in HCC with potential therapeutic targets of the active components in EZW. As shown in Figure 4C, a total of 66 targets overlapped; presumably, EZW may suppress HCC by regulating these 66 genes. Notably, both herbs in EZW act on HCC through 26 common targets, which is speculated to be the reason FLL and EH can exert synergistic anti-HCC pharmacological effects. Moreover, 35 of these 66 genes were abnormally low expressed in HCC and mainly involved metabolism-related genes, including CYP3A4, XDH, ARG1, ADRA1B, and ALDH2. In contrast, there were 31 genes expressed at abnormally high levels, mainly including genes involved in the cell cycle and related to proliferation such as CHEK1, CCNA2, CDK4, CCNB2, CCNB1, CDK1, PCNA, and MMP9 (Figure 4D). EZW may treat HCC by reversing the expression patterns of these genes, but further confirmation is needed.

To further understand the interconnections between the herbs, active compounds, and potential therapeutic targets for HCC, we generated an herb-compound-target (H-C-T) network (Figure 5A). The results revealed that the two herbs exerted anti-HCC effects on multiple targets mainly through active ingredients such as quercetin, luteolin, kaempferol, demethylwedelolactone, wedelolactone, oleanic acid, and ursolic acid. Moreover, the results indicated that ESR1, AR, CCNA2, PTGS2, and CA2 were most regulated by multiple active components of EZW. The H-C-T network revealed an intricate molecular network involved in the suppression of HCC by EZW. To further investigate the intrinsic connectivity of the therapeutic targets of EZW against HCC, we generated a protein-protein interaction (PPI) network and performed Minimal Common Oncology Data Elements (MCODE) analysis and annotation on this network. The PPI network contained 60 nodes, 136 connections, and 3 MCODE networks (Figure 5B). We performed pathway and process enrichment analysis for each MCODE component and retained the top 3 terms with the lowest P values as functional descriptions of the corresponding components (Supplementary File 4) MCODE1 was used to annotate cell cycle-related pathways, while MCODE3 was used to annotate metabolism-related pathways.

Network and topological analyses of 66 potential therapeutic targets for EZW in HCC. (A) Herb compound-target network analysis. (B) Protein-protein interaction network and gene clustering analysis (Metascape web tool). (C) Identification of the top 10 hub genes in the PPI network by different topological calculation methods.

Figure 5. Network and topological analyses of 66 potential therapeutic targets for EZW in HCC. (A) Herb compound-target network analysis. (B) Protein-protein interaction network and gene clustering analysis (Metascape web tool). (C) Identification of the top 10 hub genes in the PPI network by different topological calculation methods.

To explore hub genes in the PPI network, we performed a topological analysis of this network. The PPI network was entered into Cytoscape software to calculate the topological parameters of the nodes in the network. We calculated the top 10 hub genes in the PPI network using a different approach via the CytoHubba plugin and visualized the network via Cytoscape. These calculation methods included the degree method, maximum neighborhood component (MNC), maximal clique centrality (MCC), and density of maximum neighborhood component (DMNC). The results showed that regardless of the calculation method, the genes in the MCODE1 and MCODE2 components of the top ten hub genes were highly enriched, and mainly included CDK1, CDK4, CCNB1, CCNB2, CCNA2, PCNA, AURKA, and AURKB (Figure 5C).

Functional enrichment analysis of the potential therapeutic targets of EZW for the treatment of HCC

To fully reveal the mechanism underlying the treatment of HCC by EZW, we performed GO and KEGG pathway analyses. The top 25 GO terms revealed the functions of the potential therapeutic targets of EZW for the treatment of HCC. Annotation analysis of CCs revealed that these therapeutic targets were mainly involved in “cyclin-dependent protein (GO:0000307)”, “serine/threonine protein kinase complex (GO:1902554)”, and “condensed chromosome (GO:0000793)” (Figure 6A). For MF analysis, these therapeutic targets were mainly involved in “histone kinase activity (GO:0035173)”, “enzyme binding (GO:0019899)”, and “catalytic activity (GO:0003824)” (Figure 6B). Moreover, the top 25 GO enrichment analysis results showed that these therapeutic targets were mainly enriched in the following BPs: “response to lipid (GO:0033993)”, “cell proliferation (GO:0008283)”, “G2/M transition of mitotic cell cycle (GO:0044839)”, “regulation of cell cycle (GO:0051726)”, and “fatty acid derivative metabolic process (GO:1901568)” (Figure 6C). Furthermore, we performed a secondary classification of all enriched GO terms, which showed that within the BP category, the GO terms were mainly involved in cellular processes, metabolic processes, and biological regulation (Figure 6D).

GO enrichment analysis of 66 potential therapeutic targets for EZW in HCC. (A) Cellular components. (B) Molecular functions. (C) Biological processes. (D) Secondary classification chart of GO enrichment terms.

Figure 6. GO enrichment analysis of 66 potential therapeutic targets for EZW in HCC. (A) Cellular components. (B) Molecular functions. (C) Biological processes. (D) Secondary classification chart of GO enrichment terms.

In the KEGG pathway enrichment analysis, the results indicated that EZW exerted its effect on HCC mainly through the following pathways: “Cell cycle (ko04110)”, “Progesterone-mediated oocyte maturation (ko04914)”, “Steroid hormone biosynthesis (ko00140)”, “p53 signaling pathway (ko04115)”, “Cellular senescence (ko04218)”, “TNF signaling pathway (ko04668)”, “Metabolism of xenobiotics by cytochrome P450 (ko00980)”, “IL-17 signaling pathway (ko04657)”, “Hepatitis B (ko05161)”, and “FoxO signaling pathway (ko04068)” (Figure 7A). These top 25 KEGG pathways were mainly divided into 6 categories, including metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases (Figure 7B). To gain a more comprehensive understanding of the mechanism by which EZW suppresses HCC, we next performed a secondary classification of all KEGG pathway annotations, and the results are shown in Figure 7C. In the metabolism category, the KEGG pathways were mainly enriched in “Global and overview maps”, “Lipid metabolism”, and “Amino acid metabolism”. For the cellular process category, the KEGG pathways were mainly enriched in “Cell growth and death”, “Cellular community – eukaryotes”, and “Cell motility”.

KEGG enrichment analysis of 66 potential therapeutic targets for EZW in HCC. (A) Top 25 KEGG pathways. (B) Secondary classification of the top 25 KEGG pathways. (C) Secondary classification of all KEGG pathways. (D) Network of enriched terms colored by cluster ID analyzed by the Metascape database, where nodes that share the same cluster ID are typically close to each other. (E) Network of enriched terms colored by P value, where terms containing more genes tend to have a more significant P value.

Figure 7. KEGG enrichment analysis of 66 potential therapeutic targets for EZW in HCC. (A) Top 25 KEGG pathways. (B) Secondary classification of the top 25 KEGG pathways. (C) Secondary classification of all KEGG pathways. (D) Network of enriched terms colored by cluster ID analyzed by the Metascape database, where nodes that share the same cluster ID are typically close to each other. (E) Network of enriched terms colored by P value, where terms containing more genes tend to have a more significant P value.

To further identify the molecular mechanisms involved in the suppression of HCC by EZW, we captured these therapeutic targets through the Metascape database for pathway and process enrichment analysis. We present the correlations between different enriched terms in a network graph, and the results revealed that “PID FOXM1 pathway”, “Cell cycle”, “steroid metabolic process”, “SUMOylation”, and “cell cycle G2/M phase transition” were significantly enriched (Figure 7D, 7E).

Identification of cellular senescence-induced/inhibited genes involved in the effects of EZW on HCC and prognostic analysis

Based on the above enrichment analysis, we speculate that cell cycle pathways play a key role in the effects of EZW treatment on HCC. However, recent studies have shown that cellular senescence permanently inhibits the proliferative capacity of cells and induces irreversible cell cycle arrest, which is considered a promising strategy for the treatment of cancer [30]. Our KEGG enrichment analysis results revealed that EZW inhibition of HCC involves the cellular senescence pathway. Thus, we believe that the cellular senescence pathway plays a key role in the suppression of HCC by EZW.

To fully characterize which senescence-related genes are involved as therapeutic targets in the treatment of HCC by EZW, we first identified 153 cellular senescence-induced genes and 121 cellular senescence-inhibited genes from the cellular senescence gene database CellAge (Supplementary File 5). Subsequently, 66 therapeutic targets were assessed by Venn diagram analysis, and 4 cellular senescence-induced genes (IGFBP3, CHEK1, AXL, and AR) and 6 cellular senescence-inhibited genes (CDK4, CDK1, FOS, G6PD, AURKA, and MMP9) were assessed (Figure 8A). Moreover, we comprehensively characterized the expression patterns of these 10 senescence-related genes in HCC (Figure 8B). The results showed that the cellular senescence-induced genes AR, AXL, and IGFBP3 were abnormally expressed at low levels in HCC, while CHEK1 showed the opposite trend. Regarding the cellular senescence-inhibited genes, except FOS, which was abnormally expressed at low levels in HCC samples, they were all abnormally highly expressed.

Prognostic analysis of cellular senescence-related genes and establishment of a prognostic model. (A) Venn diagram identified cellular senescence-related genes among 66 therapeutic targets. (B) Heatmap analysis of the expression patterns of 10 cellular senescence-related genes in 14 pairs of adjacent nontumor liver tissues and hepatocellular carcinoma tissues. (C) Forest plot of univariate Cox analysis of 10 cellular senescence-related genes. (D) Correlation network of 10 cellular senescence-related genes. (E) LASSO coefficient profiles of 10 cellular senescence-related genes. (F) Cross-validation for tuning parameter selection in LASSO regression. (G) The distribution of risk scores, gene expression levels, and survival status of LIHC patients in the training cohort. (H) Kaplan–Meier curves of the OS of all LIHC patients in the TCGA cohort based on the risk score. (I) Time-dependent ROC curve analysis of the prognostic model (1, 3, and 5 years).

Figure 8. Prognostic analysis of cellular senescence-related genes and establishment of a prognostic model. (A) Venn diagram identified cellular senescence-related genes among 66 therapeutic targets. (B) Heatmap analysis of the expression patterns of 10 cellular senescence-related genes in 14 pairs of adjacent nontumor liver tissues and hepatocellular carcinoma tissues. (C) Forest plot of univariate Cox analysis of 10 cellular senescence-related genes. (D) Correlation network of 10 cellular senescence-related genes. (E) LASSO coefficient profiles of 10 cellular senescence-related genes. (F) Cross-validation for tuning parameter selection in LASSO regression. (G) The distribution of risk scores, gene expression levels, and survival status of LIHC patients in the training cohort. (H) Kaplan–Meier curves of the OS of all LIHC patients in the TCGA cohort based on the risk score. (I) Time-dependent ROC curve analysis of the prognostic model (1, 3, and 5 years).

Univariate Cox proportional hazards regression analysis revealed that seven cellular senescence-related genes were associated with the prognosis of HCC, and CDK1, CDK4, G6PD, AURKA, MMP9, IGFBP3, and CHEK1 were considered risk factors (P<0.01, HR>1) (Figure 8C). Furthermore, the expression levels of these 7 prognostic genes were strongly positively correlated with each other (Figure 8D).

Establishment of a prognostic risk scores with cellular senescence-related genes in the TCGA dataset

The above seven senescence-related genes were analyzed by least absolute shrinkage and selection operator (LASSO) Cox regression analysis to establish a cellular senescence-related signature for predicting survival. A 4-gene signature was constructed according to the optimum λ value (Figure 8E, 8F). We then established a risk score formula based on the expression of the four genes for patients with LIHC: risk score = (0.0199 + expression value of CDK1) + (0.0210 + expression value of CDK4) + (0.0665 + expression value of CHEK1) + (0.2228 + expression value of G6PD). The risk score of each patient was then calculated using this formula, and patients in the TCGA cohort were stratified into low- and high-risk groups according to the median value of the risk score. The distribution of the cellular senescence-related signature score, the survival status, and a heatmap exhibiting the expression profiles of the 4 genes in the high- and low-risk groups are presented in Figure 8G. Kaplan-Meier survival analysis demonstrated that patients in the high-risk group had a significantly shorter overall survival (OS) times than those in the low-risk group (Figure 8H, HR = 1.77 (1.25-2.50), P = 0.001). Subsequently, time-dependent receiver operating characteristic (ROC) analysis was performed, which showed that the risk score performed well in predicting 1-, 3-, and 5-year OS, with areas under the curve (AUCs) of 0.763, 0.662, and 0.614, respectively (Figure 8I).

Analysis of the clinical relevance of the 4 cellular senescence-related genes in LIHC

We analyzed the clinical characteristics of the 4 cellular senescence-related genes (CDK1, CDK4, CHEK1, and G6PD) involved in the treatment of HCC by EZW to further evaluate the clinical application value of EZW in HCC treatment. To investigate the role of these genes in LIHC, we assessed RNA-seq data obtained from 374 HCC patient tissues and 50 normal tissues using transcriptional data from the TCGA-LIHC database. The results showed that the levels of all four genes were significantly higher in tumor tissues than in normal liver tissues (P < 0.05), which was consistent with the previous GEO dataset results (Figure 9A). The transcription levels of the four cellular senescence-related genes in LIHC patients were significantly higher than those in the normal group and were closely related to clinical features such as T stage, pathological stage, and vascular invasion (Figure 9B–9D). Regarding tumor status, the transcript levels of these four genes were downregulated in the ‘tumor-free’ group relative to the ‘with tumor’, while CDK1 and CHEK1 showed significance (Figure 9E). Next, the ROC analysis results showed that the AUCs achieved using the expression levels of CDK1, CDK4, G6PD and CHEK1 were 0.976, 0.885, 0.949 and 0.951, respectively, indicating that these 4 cellular senescence-related genes exhibited adequate predictive performance (Figure 9F). OS analysis showed that patients with high expression levels of CDK4 (P < 0.001, HR = 1.92 (1.35-2.73)), CDK1 (P < 0.001, HR = 1.93 (1.36-2.74)), G6PD (P < 0.001, HR = 1.89 (1.33-2.68)), and CHEK1 (P = 0.001, HR = 1.81 (1.28-2.57)) had poorer OS than patients with low expression levels of these genes (Figure 9G).

Correlation analysis of the expression of four key cellular senescence-related genes with the clinical characteristics of LIHC patients. (A) The differential expression of CDK1, CDK4, CHEK1, and G6PD between normal and tumor tissues. (B, C) CDK1, CDK4, CHEK1, and G6PD mRNA expression in normal individuals or in patients with different T stages (T1&T2 and T3&T4) and pathologic stages (stage I&II and stage III&IV). (D, E) Differences in the expression of CDK1, CDK4, CHEK1, and G6PD mRNA according to vascular invasion and tumor status. (F) Analysis of the AUCs of the 4 cellular senescence-related genes in LIHC. (G) Kaplan-Meier curves of OS for different cell cycle-related genes. *, **, and *** represent P

Figure 9. Correlation analysis of the expression of four key cellular senescence-related genes with the clinical characteristics of LIHC patients. (A) The differential expression of CDK1, CDK4, CHEK1, and G6PD between normal and tumor tissues. (B, C) CDK1, CDK4, CHEK1, and G6PD mRNA expression in normal individuals or in patients with different T stages (T1&T2 and T3&T4) and pathologic stages (stage I&II and stage III&IV). (D, E) Differences in the expression of CDK1, CDK4, CHEK1, and G6PD mRNA according to vascular invasion and tumor status. (F) Analysis of the AUCs of the 4 cellular senescence-related genes in LIHC. (G) Kaplan-Meier curves of OS for different cell cycle-related genes. *, **, and *** represent P < 0.05, P < 0.01, and P < 0.001, respectively.

Molecular docking strategy to verify the multicomponent multitarget network of EZW in the treatment of HCC

The previous results suggest that CDK1, CDK4, CHEK1, and G6PD are key therapeutic targets through which EZW suppresses HCC. Based on the H-C-T network analysis, EZW inhibits HCC progression through multiple active components acting on these targets. To further reveal how these active ingredients act on these targets, we used a molecular docking strategy to simulate their binding modes and calculated the binding energies to infer the affinity of these active ingredients to the targets. According to the prediction results obtained using molecular docking, a variety of active ingredients (luteolin, quercetin, kaempferol, acacetin, wedelolactone, demethylwedelolactone, ursolic acid, oleanolic acid, pratensein, β-sitosterol, and 3’-O-methylorobol) could bind to these key targets (CDK1, CDK4, G6PD, and CHEK1) with low binding energies (Figure 10). As shown in the binding diagram, these active compounds can bind well in the pockets of these targets and form stable noncovalent interactions with the amino acid residues around the pockets. Presumably, because these active ingredients have multiple hydroxyl groups, they are good hydrogen bond donors or acceptors. Analysis of protein-ligand interactions based on a molecular docking strategy provides evidence for the hypothesis that EZW inhibits HCC by binding to these proteins and inhibiting their activity or expression levels.

Molecular docking analysis of different active components in EZW and four cellular senescence-related targets. Visualization of 3D binding diagrams for protein-ligand predictions based on PyMOL software. Cyan represents the surrounding amino acid residues in the binding pocket, and green represents the active ingredients.

Figure 10. Molecular docking analysis of different active components in EZW and four cellular senescence-related targets. Visualization of 3D binding diagrams for protein-ligand predictions based on PyMOL software. Cyan represents the surrounding amino acid residues in the binding pocket, and green represents the active ingredients.

Discussion

HCC is a common malignancy that usually arises in the context of chronic liver disease. Despite modern strategies for patient management, patients with advanced HCC have few treatment options and a poor prognosis. Moreover, the progression of HCC involves multiple mechanisms and aberrant changes in multiple pathological genes [31]. This complexity requires the identification of a therapeutic strategy that modulates multiple targets in HCC to improve patient outcomes. In recent years, the emergence and development of TCM has provided new opportunities for the treatment of HCC. Given that the pathological mechanism of HCC and the multicomponent and multitarget characteristics of EZW are very complex, a single experiment cannot be used to systematically and efficiently reveal the pharmacological mechanism through which EZW suppresses HCC. Through the joint analysis of multiple biological databases, the introduction of molecular docking technology and network pharmacology provides a feasible method for systematically studying the pharmacological mechanism of TCM. Therefore, this study attempted to utilize this comprehensive strategy to explore how EZW exerts pharmacological anti-HCC effects through a multiactive ingredient-multitarget network.

Exploring the pathogenesis of HCC will make an important contribution to providing suitable treatment strategies in the future. In our study, HCC progression was found to be regulated by a complex gene network involving the aberrant expression of thousands of DEGs. Notably, functional enrichment analysis of these DEGs revealed that the pathogenesis of HCC is closely related to various metabolic pathways and processes of cell cycle regulation, which have been shown to be critical for the development of liver disease [32, 33]. Among them, the metabolic pathway (ko01100) in HCC involves abnormal changes in multiple metabolic genes, including high DNMT3A, G6PD, PGP, and EZH2 expression and low AKR1D1, SLC27A5, XDH, LDHD, and ADH3 expression. Chen et al. [34] showed that EZH2 promotes HCC progression by regulating the miR-22/galectin-9 axis. Zhou et al. [35] demonstrated that the long noncoding RNA HCP5 acts as a sponge for miR-29b-3p and promotes liver cancer cell growth and metastasis by upregulating DNMT3A. The study by Cao et al. [36] showed that knockdown of G6PD in HCC reduced tumor volume and tumor weight in vivo. Data from Nikolaou et al. [37] suggest that AKR1D1 may play an important role in regulating endogenous glucocorticoid action, which may be particularly relevant to physiological and pathophysiological processes affecting the liver. Chen et al. [38] showed that XDH downregulation promotes TGFβ signaling and the expression of cancer stem cell-related genes in HCC.

Furthermore, cyclins and cyclin-dependent kinases (CDKs) are typically involved in most metabolic processes, such as glucose metabolism, lipogenesis, amino acid metabolism and mitochondrial activity [39, 40]. KEGG pathway analysis revealed that HCC progression was regulated by cell cycle pathways, and multiple cyclins (CCNB1, CCNB2, and CCNB3) and cyclin-dependent kinases (CDK1 and CDK4) were significantly activated in HCC. Previous studies have shown that overexpression of CDKs leads to abnormal cell proliferation and requires CDK activity to respond to DNA damage during DNA replication [41, 42]. Based on these data, designing small-molecule compounds targeting CDKs is considered an effective strategy to treat cancer. Palbociclib, a selective inhibitor of CDK4/6, has been FDA-approved for the treatment of breast cancer and has been shown in multiple studies to be effective in the treatment of HCC [43]. These findings suggest that the key strategy in the treatment of HCC may lie in cell cycle regulation, and the multicomponent-multitarget feature of TCM may be a promising therapy.

The extract of EZW contains a variety of compounds, which undoubtedly increases the difficulty of systematically mapping the pharmacological mechanism of EZW in inhibiting tumors. Therefore, identifying the active components of EZW through databases can be used to more accurately and efficiently elucidate its mechanism. After screening, the key active components of EZW were mainly flavonoids, coumarins, sterols, and natural triterpenoid carboxylic acids. Among them, the flavonoids quercetin and luteolin were identified as the common active components of the two herbs in EZW. We speculate that these two active components may be an important basis for the synergistic anti-HCC effect of the two herbs in EZW. Natural flavonoids have been identified as one of the major classes of natural anticancer agents, exerting antitumor activity through cell cycle arrest as a major mechanism in various cancer cells [44]. Furthermore, multiple active components in EZW have been reported to exhibit substantial anticancer activity in various tumors. Pan et al. [18] showed that coumarin wedelolactone, a characteristic component of EH, inhibited the proliferation of HCC cell lines (HepG2 and Huh-7) by inhibiting the PI3K/AKT signaling pathway. Additionally, demethylwedelolactone, a coumarin component of EH, inhibited the lung metastasis of MDA-MB-231 breast cancer cells in a nude mouse model [45]. Two representative triterpenoids in FLL, oleanolic acid and ursolic acid, were observed to induce apoptosis in various human liver cancer cell lines, indicating that they are potent anticancer agents [46]. These results suggest that the anticancer activities of these active ingredients can serve as an important theoretical basis for EZW in the treatment of liver tumorigenesis.

According to the H-C-T network diagram, EZW may exert its anti-HCC effect through the action of multiple active components on multiple targets. Subsequent PPI network, MCODE, and topological analyses indicated that EZW suppresses HCC through the regulation of multiple cyclins (CCNA1, CCNB1, and CCNB2) and cyclin-dependent kinases (CDK1 and CDK4), which were defined as hub genes. Further KEGG enrichment analysis depicting the top 25 pathways revealed that the cell cycle pathway ranked first, and our aforementioned results indicated that this pathway plays a key role in HCC progression, suggesting that this pathway is an important mechanism by which EZW inhibits HCC. However, the pharmacological mechanism of EZW also involves the p53 signaling pathway, which has been shown to play a central role in regulating the cell division cycle [47]. These findings suggest that EZW may regulate the cell cycle through the p53 pathway, thereby inhibiting HCC. Although these speculations need further verification, they still provide directions for future research on the molecular mechanism of EZW.

Interestingly, the pharmacological mechanism by which EZW inhibits HCC also involves the process of cellular senescence. Cellular senescence constitutes a permanent state of cell cycle arrest in proliferating cells induced by different stresses and has been recognized in recent years as an important mechanism for preventing tumor cell proliferation [30]. TCM with multicomponent and multitarget characteristics is believed to induce cell senescence by activating or inactivating oncogenes, inducing SASP, and triggering DNA damage, thereby inhibiting the occurrence and development of tumors [48]. Leveraging these properties has become a new direction in antitumor research; however, the role of cellular senescence in the treatment of HCC by EZW has been largely underexplored, so a broader understanding of the links among HCC, EZW, and senescence is important. Importantly, we assessed 10 cellular senescence-related genes, including 6 cellular senescence-inhibited genes (CDK4, CDK1, FOS, G6PD, AURKA, and MMP9) and 4 cellular senescence-induced genes (IGFBP3, CHEK1, AXL, and AR). Among them, CDK1, CDK4, G6PD, AURKA, MMP9, IGFBP3, and CHEK1 were considered risk factors in LIHC, and 6 of the 7 genes (CDK1, CHEK1, MMP9, G6PD, CDK4, and AURKA) were abnormally highly expressed in HCC patients. Subsequently, the LASSO Cox regression model and Kaplan-Meier survival curve analysis showed that LIHC patients with high expression of senescence-related genes (CDK1, CDK4, CHEK1, and G6PD) had a poorer prognosis than patients with low expression of these genes. In addition, the expression levels of these genes were closely related to the progression of clinical features such as tumor status, T stage, pathological stage, and vascular invasion of LIHC.

These data strongly suggest that CDK1, CDK4, CHEK1, and G6PD are potential prognostic biomarkers in HCC and key therapeutic targets for EZW to suppress HCC based on the cellular senescence process. Wu et al. [49] reported that blocking CDK1/PDK1/β-Catenin signaling with the CDK1 inhibitor RO3306 could improve the efficacy of sorafenib in the treatment of HCC. A selective CDK4/6 inhibitor, palbociclib, has been shown in recent years to inhibit cell proliferation in human hepatoma cell lines by promoting reversible cell cycle arrest. Moreover, palbociclib alone or in combination with sorafenib, the standard treatment for HCC, impairs tumor growth in vivo and significantly improves survival [43]. Previous reports have demonstrated that it is overexpressed and associated with poor prognosis in HCC, suggesting that it is an oncogene and that CHEK1 is negatively regulated by miR-497 and providing a potential molecular target for HCC therapy [50]. Altered metabolism is one of the hallmarks of cancer cells. G6PD levels have been shown to be elevated in many cancers, and it has been shown that G6PD induces epithelial-mesenchymal transition by activating the signal transducer and activator of transcription 3 pathway, thereby promoting HCC migration and hepatoma cell invasion [51]. Although the mechanisms of these genes during cellular senescence remain to be revealed, they still provide potential molecular targets for HCC therapy.

Our study suggests that EZW can be used to treat HCC by acting on these molecular targets, but its intrinsic mechanism needs to be further elucidated. Therefore, in silico simulations were used in this study to elucidate the interaction mode between the active components of EZW and these molecular targets. Since the extract of EZW contains various compounds, systematically depicting a clear pharmacological mechanism of how EZW inhibits tumors remains a challenge. The results of molecular docking revealed that multiple active compounds of EZW could be well combined in the pockets of CDK1, CDK4, CHEK1, and G6PD. We speculate that these active compounds may inhibit their activities by targeting these molecular targets, thereby disturbing the cell cycle and metabolism to suppress HCC.

Conclusions

Overall, we elucidated the potential pharmacological mechanisms and molecular targets of EZWs in HCC therapy by integrating multiple databases and performing molecular docking analysis. EZW can regulate tumor progression through multiple metabolic pathways, the cell cycle, and cellular senescence. In particular, TCM-induced cellular senescence may become a promising cancer treatment strategy. TCGA data analysis showed that CXP may improve the prognosis and clinical outcomes of HCC patients by regulating cellular senescence-related genes, revealing its clinical application value. In addition, the CDK1, CDK4, CHEK1, and G6PD genes were identified as key therapeutic targets for EZW in the treatment of HCC. Our results suggest that EZW can suppress tumors by disrupting the cell cycle, senescence, and metabolism through these therapeutic targets. To our knowledge, this is the first systematic pharmacological study of EZW in HCC therapy. Therefore, although there are still limitations in this study, this study provides innovative research methods and breakthroughs for TCM research.

Materials and Methods

Screening the chemical components of EZW

We determined the chemical components of EZW through the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP, https://old.tcmsp-e.com/tcmsp.php). First, the Chinese terms “Mo-Han-Lian” and “ Nü-Zhen-Zi” of Ecliptae Herba and Fructus Ligustri Lucidi were separately entered into the database to determine their components and retrieve their pharmacokinetic data. Here, we selected two pharmacokinetic parameters as screening criteria to identify the active components in these herbs based on previous studies [26]. Chemical components meeting the criteria of oral bioavailability (OB) ≥ 30% and drug similarity (DL) ≥ 0.18 were considered active ingredients of EZW for subsequent analysis. Moreover, invalid components were removed, including duplicate structures and compounds that could not be retrieved by PubChem (https://pubchem.ncbi.nlm.nih.gov/).

Screening therapeutic targets of EZW in the treatment of HCC

First, we screened the potential therapeutic targets of EZW’s active ingredients based on the HERB and SwissTargetPrediction databases. HERB (http://herb.ac.cn/) is a high-throughput experiment- and reference-guided database of traditional Chinese medicine. SwissTargetPrediction (http://www.swisstargetprediction.ch/) is an online database that predicts targets of biologically active small molecules. To identify the pathological genes associated with HCC progression, we downloaded GSE84402 from the GEO (http://www.ncbi.nlm.nih.gov/geo/) database. The expression data obtained from 14 pairs of human HCC tissues and corresponding noncancerous tissues were analyzed using the “limma” R package of Bioconductor (https://bioconductor.org/packages/release/bioc/html/limma.html), and the differentially expressed genes of HCC were screened with false discovery rate (FDR) < 0.05 and |log2FC| > 1. Subsequently, the “ggplot2” and “ComplexHeatmap” packages in R language were used for volcano plot and heatmap visualization, respectively. Finally, the potential targets of EZW and the DEGs of HCC were subjected to Venn diagram analysis to obtain overlapping targets. These overlapping targets are the therapeutic targets of EZW in the treatment of HCC.

Publicly attainable expression datasets

The RNA sequencing (RNA-Seq) expression profile dataset of 374 HCC patients, which included data on clinicopathological characteristics and survival, was downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov). Statistical analysis and visualization of gene expression between adjacent nontumor tissues (n = 50) and tumor tissues (n = 374) were performed using the “ggplot2” R package. The list of cellular senescence-related genes was obtained from the cellular senescence gene database CellAge (https://genomics.senescence.info/cells/), including 153 senescence-induced genes and 121 senescence-inhibited genes.

GO and KEGG pathway enrichment analyses

OmicShare (https://www.omicshare.com/tools), an online data analysis and visualization platform for KEGG pathway and GO enrichment analyses, was used. OmicShare can be also used for Venn diagrams, heatmaps, network building, volcano map analysis, and more. In addition, we used the Metascape (https://metascape.org/gp/index.html#/main/step1) database for pathway and process enrichment analysis and generated network maps.

Network construction, key module selection, hub gene identification, and topology analysis

A list of 66 therapeutic targets for EZW in the treatment of HCC was entered into the Metascape database (species limited to “Homo sapiens”) to generate a protein-protein interaction (PPI) network. Moreover, the Molecular Complex Detection (MCODE) algorithm was used to identify densely connected network components if the network contained 3 to 500 proteins, as these MCODEs are likely to represent densely connected regions in large PPI networks of molecular complexes. Subsequently, we performed visualization, hub gene analysis, and topological analysis of the raw data (cys format file) of the PPI network obtained by Metascape analysis using Cytoscape software. Through the CytoHubba plugin of Cytoscape software, hub gene analysis was performed on the network, including analysis methods such as degree, MNC, MCC, and DMNC. In addition, the parameters of topological features can be calculated by the Cytoscape plugin Network Analyzer, including “degree”, “intercentrality”, “closeness centrality”, “clustering coefficient”, and “topological coefficient”. The H-C-T network was analyzed and displayed through the OmicShare tool.

Construction and validation of a prognostic model involving cellular senescence-related genes and LASSO Cox regression

Univariate Cox proportional hazards regression analysis was performed to identify cellular senescence-related prognostic genes (P < 0.01), and gene interactions were visualized by the R package “circlize”. Next, least absolute shrinkage and selection operator (LASSO) Cox regression was conducted with a random seed using the R package “glmnet” to construct the risk score model for predicting survival in the training cohort. Normalized expression matrices of candidate prognostic cellular senescence-related genes were set as independent variables in the regression, and the response variables were OS and patient status in the TCGA cohort. A risk score was determined for each patient based on the normalized expression level of each gene and its corresponding regression coefficient. The formula is as follows:

Risk score=∑i=1nCoefficient(i)∗expression level(i)

The patients were then divided into low-risk and high-risk groups based on the median risk score.

Molecular docking analysis

The X-ray crystal structures of the selected proteins were obtained from the Protein Data Bank (PDB, https://www.rcsb.org/), and water molecules and heteroatoms were removed by PyMOL 1.8 software. Moreover, the 3D chemical structures of the active ingredients of EZW were downloaded in SDF format from PubChem and converted to ‘pdb’ format by PyMOL 1.8. Next, the proteins and active ingredients were converted to ‘pdbqt’ format files by AutoDockTools (version 1.5.6), and the grid box feature of AutoDockTools was used to define specific docking pockets in the selected proteins to which the active ingredients could bind. Once all data were prepared, the command prompt was used to perform molecular docking analysis and visualize the docking results with PyMOL.

Statistical analysis

Data analysis and graph generation were all performed in R version 3.6.3, GraphPad Prism 7.0, and the OmicShare webtool. The statistical significance of normally distributed variables was analyzed by unpaired Student’s t test, and the Wilcoxon rank sum test was used to assess nonnormally distributed variables. ROC curves for 1-, 3-, and 5-year survival to evaluate the predictive efficacy of the risk score were generated using the “timeROC” R package. Moreover, the ROC curves of different genes were analyzed by the “pROC” package and visualized with the “ggplot2” package in R. Kaplan-Meier survival curves for overall survival (OS) analysis were drawn using the R package “survminer”. Two-tailed P values < 0.05 were considered statistically significant.

Data availability

Publicly available datasets were analyzed in this study. This data can be found here: TCMSP, TCGA, GEO, etc.

Abbreviations

HCC: breast cancer; EZW: Er-Zhi-Wan; FLL: Fructus Ligustri Lucidi; EH: Ecliptae Herba; TCM: traditional Chinese medicine; LKYDS: liver-kidney Yin deficiency syndrome; DEGs: differentially expressed genes; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; OS: overall survival; PPI: protein-protein interaction; ROC: receiver operating characteristic; AUC: areas under the curve; FDR: false discovery rate; LASSO: Least absolute shrinkage and selection operator; BP: biological process; CC: cytological component; MF: molecular function; H-AC: herb-active components; TCMSP: Traditional Chinese Medicine System Pharmacology Database; OB: oral bioavailability; DL: drug-likeness; LIHC: Liver hepatocellular carcinoma; TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus; MCODE: Molecular Complex Detection; H-C-T: Herb-Compounds-Targets.

Author Contributions

Botao Pan conceived and designed the study; Shaoyan Zheng, Botao Pan collected data, performed the data analysis; Shaoyan Zheng wrote the manuscript; Botao Pan reviewed the manuscript and provided comments. And all of the authors reviewed the manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

No funding was provided for this study.

References

  • 1. Ferlay J, Colombet M, Soerjomataram I, Mathers C, Parkin DM, Piñeros M, Znaor A, Bray F. Estimating the global cancer incidence and mortality in 2018: GLOBOCAN sources and methods. Int J Cancer. 2019; 144:1941–53. https://doi.org/10.1002/ijc.31937 [PubMed]
  • 2. Ferlay J, Colombet M, Soerjomataram I, Parkin DM, Piñeros M, Znaor A, Bray F. Cancer statistics for the year 2020: An overview. Int J Cancer. 2021. [Epub ahead of print]. https://doi.org/10.1002/ijc.33588 [PubMed]
  • 3. Sangro B, Sarobe P, Hervás-Stubbs S, Melero I. Advances in immunotherapy for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2021; 18:525–43. https://doi.org/10.1038/s41575-021-00438-0 [PubMed]
  • 4. Craig AJ, von Felden J, Garcia-Lezana T, Sarcognato S, Villanueva A. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2020; 17:139–52. https://doi.org/10.1038/s41575-019-0229-4 [PubMed]
  • 5. Torimura T, Iwamoto H. Treatment and the prognosis of hepatocellular carcinoma in Asia. Liver Int. 2022; 42:2042–54. https://doi.org/10.1111/liv.15130 [PubMed]
  • 6. Llovet JM, Zucman-Rossi J, Pikarsky E, Sangro B, Schwartz M, Sherman M, Gores G. Hepatocellular carcinoma. Nat Rev Dis Primers. 2016; 2:16018. https://doi.org/10.1038/nrdp.2016.18 [PubMed]
  • 7. Moon H, Choi JE, Lee IJ, Kim TH, Kim SH, Ko YH, Kim HB, Nam BH, Park JW. All-treatment array of hepatocellular carcinoma from initial diagnosis to death: observation of cumulative treatments. J Cancer Res Clin Oncol. 2017; 143:2327–39. https://doi.org/10.1007/s00432-017-2480-9 [PubMed]
  • 8. Villanueva A. Hepatocellular Carcinoma. N Engl J Med. 2019; 380:1450–62. https://doi.org/10.1056/NEJMra1713263 [PubMed]
  • 9. Sonbol MB, Riaz IB, Naqvi SAA, Almquist DR, Mina S, Almasri J, Shah S, Almader-Douglas D, Uson Junior PLS, Mahipal A, Ma WW, Jin Z, Mody K, et al. Systemic Therapy and Sequencing Options in Advanced Hepatocellular Carcinoma: A Systematic Review and Network Meta-analysis. JAMA Oncol. 2020; 6:e204930. https://doi.org/10.1001/jamaoncol.2020.4930 [PubMed]
  • 10. Alqahtani A, Khan Z, Alloghbi A, Said Ahmed TS, Ashraf M, Hammouda DM. Hepatocellular Carcinoma: Molecular Mechanisms and Targeted Therapies. Medicina (Kaunas). 2019; 55:526. https://doi.org/10.3390/medicina55090526 [PubMed]
  • 11. Liu X, Li M, Wang X, Dang Z, Yu L, Wang X, Jiang Y, Yang Z. Effects of adjuvant traditional Chinese medicine therapy on long-term survival in patients with hepatocellular carcinoma. Phytomedicine. 2019; 62:152930. https://doi.org/10.1016/j.phymed.2019.152930 [PubMed]
  • 12. Lu HL, Su YC, Lin MC, Sun MF, Huang ST. Integrating Chinese and Western medicines reduced the incidence of hepatocellular carcinoma in patients with diabetes mellitus: A Taiwanese population-based cohort study. Complement Ther Med. 2020; 49:102332. https://doi.org/10.1016/j.ctim.2020.102332 [PubMed]
  • 13. Sun C, Dong F, Xiao T, Gao W. Efficacy and safety of Chinese patent medicine (Kang-ai injection) as an adjuvant in the treatment of patients with hepatocellular carcinoma: a meta-analysis. Pharm Biol. 2021; 59:472–83. https://doi.org/10.1080/13880209.2021.1915340 [PubMed]
  • 14. Zhang W, Zhang B, Chen XP. Adjuvant treatment strategy after curative resection for hepatocellular carcinoma. Front Med. 2021; 15:155–69. https://doi.org/10.1007/s11684-021-0848-3 [PubMed]
  • 15. Yuan S, Gong Y, Chen R, Du J, Zhang H, Chen T. Chinese herbal formula QHF inhibits hepatocellular carcinoma metastasis via HGF/c-Met signaling pathway. Biomed Pharmacother. 2020; 132:110867. https://doi.org/10.1016/j.biopha.2020.110867 [PubMed]
  • 16. Hu Y, Wang S, Wu X, Zhang J, Chen R, Chen M, Wang Y. Chinese herbal medicine-derived compounds for cancer therapy: a focus on hepatocellular carcinoma. J Ethnopharmacol. 2013; 149:601–12. https://doi.org/10.1016/j.jep.2013.07.030 [PubMed]
  • 17. Li JJ, Liang Q, Sun GC. Traditional Chinese medicine for prevention and treatment of hepatocellular carcinoma: A focus on epithelial-mesenchymal transition. J Integr Med. 2021; 19:469–77. https://doi.org/10.1016/j.joim.2021.08.004 [PubMed]
  • 18. Pan B, Pan W, Lu Z, Xia C. Pharmacological Mechanisms Underlying the Hepatoprotective Effects of Ecliptae herba on Hepatocellular Carcinoma. Evid Based Complement Alternat Med. 2021; 2021:5591402. https://doi.org/10.1155/2021/5591402 [PubMed]
  • 19. Zhai Y, Xu J, Feng L, Liu Q, Yao W, Li H, Cao Y, Cheng F, Bao B, Zhang L. Broad range metabolomics coupled with network analysis for explaining possible mechanisms of Er-Zhi-Wan in treating liver-kidney Yin deficiency syndrome of Traditional Chinese medicine. J Ethnopharmacol. 2019; 234:57–66. https://doi.org/10.1016/j.jep.2019.01.019 [PubMed]
  • 20. Lee S, Park J, Lee H, Kim K. Development and validation of Yin-Deficiency Questionnaire. Am J Chin Med. 2007; 35:11–20. https://doi.org/10.1142/S0192415X07004576 [PubMed]
  • 21. Yao W, Gu H, Zhu J, Barding G, Cheng H, Bao B, Zhang L, Ding A, Li W. Integrated plasma and urine metabolomics coupled with HPLC/QTOF-MS and chemometric analysis on potential biomarkers in liver injury and hepatoprotective effects of Er-Zhi-Wan. Anal Bioanal Chem. 2014; 406:7367–78. https://doi.org/10.1007/s00216-014-8169-x [PubMed]
  • 22. Hu B, Du Q, Deng S, An HM, Pan CF, Shen KP, Xu L, Wei MM, Wang SS. Ligustrum lucidum Ait. fruit extract induces apoptosis and cell senescence in human hepatocellular carcinoma cells through upregulation of p21. Oncol Rep. 2014; 32:1037–42. https://doi.org/10.3892/or.2014.3312 [PubMed]
  • 23. Nogales C, Mamdouh ZM, List M, Kiel C, Casas AI, Schmidt HHH. Network pharmacology: curing causal mechanisms instead of treating symptoms. Trends Pharmacol Sci. 2022; 43:136–50. https://doi.org/10.1016/j.tips.2021.11.004 [PubMed]
  • 24. Zhang R, Zhu X, Bai H, Ning K. Network Pharmacology Databases for Traditional Chinese Medicine: Review and Assessment. Front Pharmacol. 2019; 10:123. https://doi.org/10.3389/fphar.2019.00123 [PubMed]
  • 25. He S, Wang T, Shi C, Wang Z, Fu X. Network pharmacology-based approach to understand the effect and mechanism of Danshen against anemia. J Ethnopharmacol. 2022; 282:114615. https://doi.org/10.1016/j.jep.2021.114615 [PubMed]
  • 26. Zhao M, Pan B, He Y, Niu B, Gao X. Elucidating the pharmacological mechanism by which Si-Wu-Tang induces cellular senescence in breast cancer via multilevel data integration. Aging (Albany NY). 2022; 14:5812–37. https://doi.org/10.18632/aging.204185 [PubMed]
  • 27. Tian S, Wang J, Li Y, Xu X, Hou T. Drug-likeness analysis of traditional Chinese medicines: prediction of drug-likeness using machine learning approaches. Mol Pharm. 2012; 9:2875–86. https://doi.org/10.1021/mp300198d [PubMed]
  • 28. Jia CY, Li JY, Hao GF, Yang GF. A drug-likeness toolbox facilitates ADMET study in drug discovery. Drug Discov Today. 2020; 25:248–58. https://doi.org/10.1016/j.drudis.2019.10.014 [PubMed]
  • 29. Chen B, Wang L, Li L, Zhu R, Liu H, Liu C, Ma R, Jia Q, Zhao D, Niu J, Fu M, Gao S, Zhang D. Fructus Ligustri Lucidi in Osteoporosis: A Review of its Pharmacology, Phytochemistry, Pharmacokinetics and Safety. Molecules. 2017; 22:1469. https://doi.org/10.3390/molecules22091469 [PubMed]
  • 30. Wang L, Lankhorst L, Bernards R. Exploiting senescence for the treatment of cancer. Nat Rev Cancer. 2022; 22:340–55. https://doi.org/10.1038/s41568-022-00450-9 [PubMed]
  • 31. Huang R, Liu J, Li H, Zheng L, Jin H, Zhang Y, Ma W, Su J, Wang M, Yang K. Identification of Hub Genes and Their Correlation With Immune Infiltration Cells in Hepatocellular Carcinoma Based on GEO and TCGA Databases. Front Genet. 2021; 12:647353. https://doi.org/10.3389/fgene.2021.647353 [PubMed]
  • 32. Piccinin E, Villani G, Moschetta A. Metabolic aspects in NAFLD, NASH and hepatocellular carcinoma: the role of PGC1 coactivators. Nat Rev Gastroenterol Hepatol. 2019; 16:160–74. https://doi.org/10.1038/s41575-018-0089-3 [PubMed]
  • 33. Matsuda Y, Wakai T, Kubota M, Takamura M, Yamagiwa S, Aoyagi Y, Osawa M, Fujimaki S, Sanpei A, Genda T, Ichida T. Clinical significance of cell cycle inhibitors in hepatocellular carcinoma. Med Mol Morphol. 2013; 46:185–92. https://doi.org/10.1007/s00795-013-0047-7 [PubMed]
  • 34. Chen S, Pu J, Bai J, Yin Y, Wu K, Wang J, Shuai X, Gao J, Tao K, Wang G, Li H. EZH2 promotes hepatocellular carcinoma progression through modulating miR-22/galectin-9 axis. J Exp Clin Cancer Res. 2018; 37:3. https://doi.org/10.1186/s13046-017-0670-6 [PubMed]
  • 35. Zhou Y, Li K, Dai T, Wang H, Hua Z, Bian W, Wang H, Chen F, Ai X. Long non-coding RNA HCP5 functions as a sponge of miR-29b-3p and promotes cell growth and metastasis in hepatocellular carcinoma through upregulating DNMT3A. Aging (Albany NY). 2021; 13:16267–86. https://doi.org/10.18632/aging.203155 [PubMed]
  • 36. Cao F, Luo A, Yang C. G6PD inhibits ferroptosis in hepatocellular carcinoma by targeting cytochrome P450 oxidoreductase. Cell Signal. 2021; 87:110098. https://doi.org/10.1016/j.cellsig.2021.110098 [PubMed]
  • 37. Nikolaou N, Gathercole LL, Kirkwood L, Dunford JE, Hughes BA, Gilligan LC, Oppermann U, Penning TM, Arlt W, Hodson L, Tomlinson JW. AKR1D1 regulates glucocorticoid availability and glucocorticoid receptor activation in human hepatoma cells. J Steroid Biochem Mol Biol. 2019; 189:218–27. https://doi.org/10.1016/j.jsbmb.2019.02.002 [PubMed]
  • 38. Chen GL, Ye T, Chen HL, Zhao ZY, Tang WQ, Wang LS, Xia JL. Xanthine dehydrogenase downregulation promotes TGFβ signaling and cancer stem cell-related gene expression in hepatocellular carcinoma. Oncogenesis. 2017; 6:e382. https://doi.org/10.1038/oncsis.2017.81 [PubMed]
  • 39. Leal-Esteban LC, Fajas L. Cell cycle regulators in cancer cell metabolism. Biochim Biophys Acta Mol Basis Dis. 2020; 1866:165715. https://doi.org/10.1016/j.bbadis.2020.165715 [PubMed]
  • 40. Aguilar V, Fajas L. Cycling through metabolism. EMBO Mol Med. 2010; 2:338–48. https://doi.org/10.1002/emmm.201000089 [PubMed]
  • 41. Kciuk M, Gielecińska A, Mujwar S, Mojzych M, Kontek R. Cyclin-dependent kinases in DNA damage response. Biochim Biophys Acta Rev Cancer. 2022; 1877:188716. https://doi.org/10.1016/j.bbcan.2022.188716 [PubMed]
  • 42. Roskoski R Jr. Cyclin-dependent protein serine/threonine kinase inhibitors as anticancer drugs. Pharmacol Res. 2019; 139:471–88. https://doi.org/10.1016/j.phrs.2018.11.035 [PubMed]
  • 43. Bollard J, Miguela V, Ruiz de Galarreta M, Venkatesh A, Bian CB, Roberto MP, Tovar V, Sia D, Molina-Sánchez P, Nguyen CB, Nakagawa S, Llovet JM, Hoshida Y, Lujambio A. Palbociclib (PD-0332991), a selective CDK4/6 inhibitor, restricts tumour growth in preclinical models of hepatocellular carcinoma. Gut. 2017; 66:1286–96. https://doi.org/10.1136/gutjnl-2016-312268 [PubMed]
  • 44. Ren W, Qiao Z, Wang H, Zhu L, Zhang L. Flavonoids: promising anticancer agents. Med Res Rev. 2003; 23:519–34. https://doi.org/10.1002/med.10033 [PubMed]
  • 45. Lee YJ, Lin WL, Chen NF, Chuang SK, Tseng TH. Demethylwedelolactone derivatives inhibit invasive growth in vitro and lung metastasis of MDA-MB-231 breast cancer cells in nude mice. Eur J Med Chem. 2012; 56:361–7. https://doi.org/10.1016/j.ejmech.2012.07.041 [PubMed]
  • 46. Yan SL, Huang CY, Wu ST, Yin MC. Oleanolic acid and ursolic acid induce apoptosis in four human liver cancer cell lines. Toxicol In Vitro. 2010; 24:842–8. https://doi.org/10.1016/j.tiv.2009.12.008 [PubMed]
  • 47. Engeland K. Cell cycle regulation: p53-p21-RB signaling. Cell Death Differ. 2022; 29:946–60. https://doi.org/10.1038/s41418-022-00988-z [PubMed]
  • 48. Liu Y, Yang S, Wang K, Lu J, Bao X, Wang R, Qiu Y, Wang T, Yu H. Cellular senescence and cancer: Focusing on traditional Chinese medicine and natural products. Cell Prolif. 2020; 53:e12894. https://doi.org/10.1111/cpr.12894 [PubMed]
  • 49. Wu CX, Wang XQ, Chok SH, Man K, Tsang SHY, Chan ACY, Ma KW, Xia W, Cheung TT. Blocking CDK1/PDK1/β-Catenin signaling by CDK1 inhibitor RO3306 increased the efficacy of sorafenib treatment by targeting cancer stem cells in a preclinical model of hepatocellular carcinoma. Theranostics. 2018; 8:3737–50. https://doi.org/10.7150/thno.25487 [PubMed]
  • 50. Xie Y, Wei RR, Huang GL, Zhang MY, Yuan YF, Wang HY. Checkpoint kinase 1 is negatively regulated by miR-497 in hepatocellular carcinoma. Med Oncol. 2014; 31:844. https://doi.org/10.1007/s12032-014-0844-4 [PubMed]
  • 51. Lu M, Lu L, Dong Q, Yu G, Chen J, Qin L, Wang L, Zhu W, Jia H. Elevated G6PD expression contributes to migration and invasion of hepatocellular carcinoma cells by inducing epithelial-mesenchymal transition. Acta Biochim Biophys Sin (Shanghai). 2018; 50:370–80. https://doi.org/10.1093/abbs/gmy009 [PubMed]
Download PDF
Corresponding Author

Botao Pan
pbtpku@163.com
https://orcid.org/0000-0002-9678-1512

Keywords
Er-Zhi-Wan (EZW) hepatocellular carcinoma (HCC) cell cycle cellular senescence molecular docking
Paper Sections