Abstract
Background. Gecko has been widely documented in Chinese scientific literature as an anti-tumor agent for various illnesses for thousands of years, and more recently, it has been examined for its anti-tumor effects on several cancers. The effect of Gecko microRNAs (miRNAs) on hepatocellular carcinoma (HCC) has not yet been reported.
Objectives. This study was designed to identify miRNAs in Gecko through small RNA sequencing and utilize bioinformatics techniques to construct a potential regulatory network and explore the possible mechanisms of exogenous miRNAs involved in HCC.
Materials and methods. RNA was extracted from Gecko tablets, and we screened the Gecko miRNA expression dataset after high-throughput sequencing. Bioinformatics analysis was used to identify novel Gecko and HCC survival-related miRNA-mRNA cross-species regulation networks.
Results. miR-100-5p, miR-99a-5p and miR-101-3p were identified as critical for the role of Geckos in HCC. Nine downstream mRNAs (EZH2, KPNA2, LMNB1, LRRC1, MRGBP, SMARCD1, STMN1, SUB1, and UBE2A) were identified as target genes for critical miRNAs. A miRNA-mRNA regulatory network was constructed, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed these key mRNAs might be associated with both the suppression and progression of HCC. The novel network significantly correlated with the abundance of multiple immune cells, as determined with immune infiltration analysis.
Conclusions. These findings suggest that Gecko may inhibit progression and exert a therapeutic effect on HCC by targeting critical miRNA-mRNA networks for cross-species regulation. It also provides a reference for future research and development of traditional Chinese medicine (TCM).
Key words: hepatocellular carcinoma, miRNA, immunity, traditional Chinese medicine, Gecko
Background
MicroRNAs (miRNAs) are small, endogenous noncoding RNAs approx. 22 bases in length that primarily exert negative regulation on target protein-coding genes. They exhibit specific complementary binding to their targets, enabling a single miRNA to target numerous genes. The mechanism of this multi-target regulation coincides with the therapeutic theory of traditional Chinese medicine (TCM).1 As tumor suppressors, miRNAs are associated with various biological processes, including cell proliferation, differentiation and apoptosis. Meanwhile, research on using miRNAs for tumor therapeutics is proliferating and has received extensive attention and investigation.2 Recent studies have revealed that miRNAs derived from herbal medicines can enter mammalian cells upon ingestion, where they stabilize and modulate gene expression, thereby influencing the functions of tissues and cells.3, 4, 5
In recent years, most of the anti-tumor studies of exogenous miRNAs have originated from medicinal plant miRNAs or compounds, while relatively few studies have been conducted on the effects of animal-based miRNAs on tumors.6 Gecko is commonly used in TCM, and numerous clinical studies have shown that it exhibits good efficacy against cancers of the digestive tract, especially esophageal, liver and gastric cancers, inhibiting tumor development and improving patient immunity.7, 8 The aqueous Gecko extract can inhibit the growth of liver cancer cells, and the Gecko protein can induce apoptosis of cancer cells without adversely affecting normal cells.9, 10 Current studies have focused only on Gecko protein crude extracts and clinical efficacy, but there has been no extraction and identification of Gecko miRNAs and exploration of their molecular anti-tumor mechanisms. Herein, we acquired the miRNA sequences of Gecko through small RNA sequencing, and we used bioinformatics data integration analysis to construct a potential regulatory network of Gecko miRNAs affecting the prognosis of hepatocellular carcinoma (HCC), which was validated using multi-pathway enrichment analysis and immune infiltration correlation.
Objectives
The aim of the study was to screen the prognostic core miRNAs of Gecko and to establish a model of Gecko miRNA regulatory network by combining the survival analysis, differential expression and immune infiltration analysis, in order to better explore the mechanism of Gecko miRNAs in HCC.
Materials and methods
Gecko miRNA sequencing
Chinese medicinal geckos were purchased from the First Hospital of the Hunan University of Chinese Medicine (Changsha, China), dried after removing the viscera and processed by powdering. The total RNA of Gecko powder was extracted with the TRIzol method to estimate RNA concentration and integrity, and sRNA libraries were constructed using the Truseq Small RNA Library Prep Kit tool. The library was tested for quality and then sequenced using the Illumina HiSeq 2500 high-throughput sequencing platform (Beijing Novogene Technology Co., Beijing, China). The miRNAs with a read count value higher than 105 and a Transcripts per million (TPM) value higher than 7,000 were included in the miRNA sequencing data of Gecko in descending order of expression. The selected miRNA base sequences were compared with human miRNAs to obtain co-expressed miRNAs for subsequent studies.
Data collection
To screen out differential miRNAs in HCC, we downloaded human miRNA expression profiles GSE147889 (including 2565 miRNAs from 97 cancer samples and 97 adjacent standard samples) from the Gene Expression Omnibus database (GEO; http://www.ncbi.nlm.nih.gov/geo),11 which is based on the GPL21263 platform (3D-Gene Human miRNA V21_1.0.0).12 Furthermore, miRNA expression profiles of HCC patients (including 374 cancer and 50 normal samples) and clinical information of tumor samples (including 374 samples) were also obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov) database to improve the reliability of the results.13
Identification of key miRNAs
Gecko and human identically expressed miRNAs were intersected with prognostically relevant differential miRNAs in HCC patients to obtain Gecko and human co-expressed key miRNAs for follow-up studies. In addition, we further validated the expression patterns of key miRNAs in multiple cancers using the CancerMIRNone (CancerMIRNone; http://bioinfo.jialab-ucr.org/CancerMIRNome/) database.14 This database is a web-based tool that encapsulates a range of cutting-edge bioinformatics tools and machine learning algorithms that allow analysis of miRNAs of interest in multiple cancer types to identify divergent miRNAs and develop diagnostics, prognostic signatures, comparisons, and visualizations.
Enrichment analysis
The obtained miRNAs were enriched using miRTarBase (https://mirtarbase.cuhk.edu.cn),15 miRDB (http://www.mirdb.org) and TargetScan (http://www.targetscan.org)16 databases for downstream target gene prediction, while the predicted downstream target genes from the 2 databases were compared with the dataset and the intersection of target genes was obtained. The regulatory interactions between miRNAs and mRNAs helped to form a basic regulatory network of miRNA-targeting mRNAs. The miRNA-mRNA targeting network was viewed using Cytoscape (v. 3.8.2) software (https://cytoscape.org).17 To better understand the roles of key Gecko miRNAs, we analyzed these miRNAs for Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment using the R package (R Foundation for Statistical Computing, Vienna, Austria).
Validation of gene expression and survival analyses
Utilizing gene expression and clinical data from the TCGA
database, we performed univariate Cox regression analysis using the “survival” package to screen for target genes significantly associated with prognosis. Subsequently, we assessed the downstream mRNA expression levels targeted by miRNAs in HCC to identify prognostically relevant upregulated target genes for further analysis. The expression differences of these genes were illustrated with boxplot, and the clinical prognostic significance between high and low expression groups was depicted with Kaplan–Meier (KM) survival curves, with the hazard ratio outlined in a forest plot. All figures were generated using RStudio software (v.2023.12.1+402; https://posit.co/products/open-source/rstudio).
Survival-related miRNA network identification and annotation
Based on the aforementioned predictions and assessments, we constructed and visualized a novel miRNA-mRNA regulatory network associated with Gecko and HCC patient prognosis using Cytoscape. By further analyzing the role of Gecko miRNAs in this regulatory network, KEGG and GO enrichment analysis of downstream mRNAs were explored through the SangerBox database (http://sangerbox.com/Tool). The DisGeNET database (http://www.disgenet.org/) was utilized to investigate the molecular underpinnings of human diseases and their complications, analyze the features of disease-related genes, and evaluate the associations between genes and diseases.18
Immune infiltration analyses
In this study, we utilized CIBERSORT (https://cibersortx.stanford.edu/) to investigate the association between Gecko miRNA regulatory networks and immune cell infiltration in HCC. CIBERSORT is a deconvolution algorithm to predict the proportion of 22 tumor-infiltrating immune cells in each sample based on gene expression data.19 The immune cell infiltration matrix was divided into 2 groups according to the median expression of key genes, and thus, the association between key genes and tumor-infiltrating immune cells was analyzed using Spearman’s correlation analysis.
Statistical analyses
R software (v. 4.3.0; http://www.r-project.org/) was used to convert miRNA probe IDs into mature miRNA names, based on the platform annotation file. Adjusted p-values for differential expression of miRNAs and genes were computed as false discovery rate (FDR) using the “limma” package.20 Wilcoxon test was employed to compare differences between 2 groups containing non-parametric data. The KM method and a log-rank test were used to measure overall survival (OS) employing the “survival” package of R. Univariate Cox regression analysis was applied to test the risk factors. Martingale residuals were used for the linearity assumption testing. The proportionality of the hazard function was checked based on the Schoenfeld residuals using the R function “cox.zph”.21 The results of all hypotheses are presented in Supplementary Table 1 and Supplementary Fig. 1. The correlation between target gene expression and immune infiltration was analyzed using Spearman’s method. The ranked data from the Spearman’s correlation analysis are detailed in Supplementary Table 4, and the results of the Spearman correlation analysis are provided in Supplementary Table 5. Only enrichment terms with p < 0.05 were deemed statistically significant and reported in our analysis.
Results
Sequencing and screening of Gecko miRNAs
By sequencing Gecko miRNAs and comparing them with the same miRNA bases in human, a total of 19 co-expressed miRNAs were obtained after taking into account the miRNAs with read count values greater than or equal to 105 and TPM values greater than 7000 in the study, including miR-1a-3p, miR-143-3p, miR-148a-3p, miR-21-5p, miR-26-5p, miR-203-3p, let-7f-5p, miR-100-5p, miR-206-3p, miR-184-3p, miR-10b-5p, miR-99a-5p, miR-140-3p, let-7e-5p, miR-122-5p, miR-9-5p, let-7i-5p, miR-101-3p, and miR-200a-3p. The whole read counts, TPM values and sequence numbers are shown in Table 1.
Certification of DE-miRNAs in HCC
Differentially expressed miRNAs in the GEO and TCGA datasets were identified using empirical Bayes moderation in the “limma” package with adjusted p < 0.05 and |log2 FC| > 0.5, as depicted in the heat maps (Figure 1, Figure 2) and volcano plots (Figure 3A). In the TCGA dataset, 394 differential miRNAs were selected, while the GEO database identified 189 differential miRNAs. Using the Venn diagram intersection (Figure 3B), 92 miRNAs were differentially expressed in both datasets and were selected for subsequent analysis.
Key miRNA identification
To assess the prognostic value of differential miRNAs on the OS of HCC, clinical information of corresponding patients, including the survival status and survival time of 373 HCC patients, was gathered from the TCGA database. Survival-related univariate Cox analysis and KM survival analysis were performed, and 24 miRNAs out of 92 differentially expressed miRNAs were significantly associated with the prognosis of HCC patients (Figure 4A). Taking the intersection of the above 24 miRNAs with the 19 co-expressed miRNAs of Gecko showed that 4 differential miRNAs were most likely to be significantly associated with the prognosis of HCC patients (miR-9-5p: p = 0.008, miR-99a-5p: p = 0.001, miR-100-5p: p < 0.001, and miR-101-3p: p = 0.032). Meanwhile, we further compared the expression levels of these survival-related miRNAs between cancer and normal samples in the TCGA dataset (Figure 4B). Among them, miR-100-5p, miR-99a-5p and miR-101-3p exhibited low expression and were associated with poor prognosis in human HCC (p < 0.001) (Figure 4C). To enhance the reliability of the screening results, we performed a pan-cancer analysis in the CancerMIRNome database to validate these 3 differential miRNAs (Figure 4D). miR-100-5p, miR-99a-5p and miR-101-3p differed in multiple cancer types, including liver cancer. Therefore, these 3 miRNAs were selected for subsequent analysis.
Construction of miRNA regulatory network
To identify the downstream mRNAs that these 3 key miRNAs might target, we applied miRTarBase, targetScan and miRDB databases to improve the confidence of the predictive results. The 3 databases intersected a total of 198 mRNAs, of which 35 downstream mRNAs were targeted by miR-99a-5p, while there were 195 downstream mRNAs likely to be targeted by miR-101-3p and 68 downstream mRNAs likely to be targeted by miR-100-5p. Based on the above predictions, the miRNA regulatory network consisting of 3 miRNAs and 198 target mRNAs were established and visualized using Cytoscape (Figure 5A). The intersection targets of the databases are presented using a Venn diagram (Figure 5B).
Further enrichment analysis
To better understand these candidate target genes, GO function, and KEGG pathway enrichment analysis were performed using R software. A total of 33 KEGG pathways and 234 GO terms were enriched, showing each group’s top 10 GO-enriched terms and the top 10 KEGG-enriched pathways in terms of p-value. In the biological process (BP), the miRNA regulatory network was mainly enriched for regulating transmembrane receptor protein serine/threonine kinase signaling pathways, cellular response to transforming growth factor β stimulation, glycerolipid biosynthesis process, and regulation of the Wnt signaling pathway. In the cellular component (CC), the network was mainly enriched in ribonucleoprotein granules, fiber centers and cytoplasmic stress granules. It demonstrated significant enrichment in molecular functions (MF), including phosphatase activity, transcriptional co-regulatory activity, activin binding, and ubiquitin-protein transferase regulatory activity. These findings suggested that the miRNA network play a crucial role in cell apoptosic, migration and proliferation (Figure 5C). Furthermore, KEGG enrichment analysis revealed significant enrichment of miRNA regulatory networks in various cancer-related pathways, including the MAPK and PI3K-Akt signaling pathways, miRNA in cancer, HCC, Th17 cell differentiation, and transcriptional dysregulation in cancer pathways (Figure 5D).
Identification of key mRNAs
Based on miRNA regulation theory, miRNAs may be adversely correlated with downstream mRNAs. Consequently, we collected the expression profiles of 198 target mRNAs in the regulatory network from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) database and assessed the expression levels of these mRNAs using the “limma” package. Additionally, the prognostic value of these mRNAs on OS was evaluated using the “survival” package. Among them, 15 target genes were found to be significantly upregulated in LIHC (Supplementary Table 2). Additionally, 38 target genes were found to be significantly associated with HCC prognosis using univariate Cox survival model analysis (Supplementary Table 3). We finally took the intersection to obtain 9 upregulated genes (EZH2, KPNA2, LMNB1, LRRC1, MRGBP, SMARCD1, STMN1, SUB1, and UBE2A). The KM survival curves and forest plots indicated that 9 genes were significantly correlated with the prognosis of patients with HCC (Figure 6A,B). The boxplot showed differential expressions (Figure 6C).
Survival-related miRNA regulatory network and enrichment analysis
Based on the above validation, a survival-related miRNA regulatory network of Gecko was established, including 3 miRNAs and 9 target genes (Figure 7A). This regulatory network showed significant prognostic value, consistent with the miRNA hypothesis. Additionally, we investigated biological functions of the target genes using the Sangerbox database. Kyoto Encyclopedia of Genes and Genomes enrichment analysis revealed significant enhancement in cancer-related pathways, including miRNA in cancer signaling, MAPK pathway signaling and apoptosis signaling pathways. Gene Ontology enrichment in BP was supplemented in the regulation of the thrombin-activated receptor signaling pathway, multi-organism nuclear import and hepatocyte growth factor receptor signaling pathway. The CC pathway was enriched in a nuclear lumen, chromatin, nuclear part, and nucleoplasm, and the MF pathway was enriched in transcription coregulator activity, histone methyltransferase activity and primary miRNA binding (Figure 7B). In the DisGeNET database, the miRNA regulatory network was considerably enriched for several cancer diseases, including follicular lymphoma, solid tumors, small cell carcinoma, and gallbladder cancer. The Gecko miRNA regulatory network significantly affects both tumor pathways and diseases (Figure 7C).
Correlation of Gecko miRNA regulatory network with immune infiltration level in HCC
CIBERSORT analysis quantified immune cell proportions in tumor samples. Tumor specimens were categorized based on median gene expression, and the impact of differential gene expression on immune cell abundance was analyzed. The results indicated a significant influence on immune cell infiltration, including memory B cells, CD4+ memory-activated T cells, T follicular helper cells, M0 macrophages, and others on high and low expression of 4 core target genes (KPNA2, STMN1, EZH2, and LMNB1) when compared to other genes (Figure 8, Figure 9). Furthermore, we evaluated the correlation between the Gecko miRNA regulatory network and immune cell infiltration using Spearman’s correlation analysis in the HCC microenvironment. As expected, we observed a negative monotonic component of correlation with naïve B cells, monocytes, M2 macrophages, and resting mast cells, and a positive monotonic component with memory B cells, CD4 memory activated T cells, T follicular helper cells, T regulatory cells, and M0 macrophages (Figure 10). This suggested a close link between the Gecko miRNA regulatory network and immune infiltration in HCC, potentially influencing immune cell function to inhibit HCC progression and development.
Discussion
Recent studies have demonstrated the potential of miRNA to therapeutically target cancer-related genes. Better targeting is provided by miRNA-based medicines, which can also be combined with other treatment methods, including chemotherapy, radiation therapy and immunotherapy, to increase their efficacy.22, 23 However, addressing several downstream targets and identifying harmful side effects are still difficult tasks. The construction and delivery of endogenous miRNA mimics are required to increase clinical usability,24 while exogenous miRNAs derived from dietary sources exhibit potential as gene expression regulators with fewer negative consequences.25 The role of exogenous miRNAs in tumor regulation has, however, received relatively little research, and there is still potential for further investigation. Moreover, various etiologies may also affect the role of systemic therapies and response markers such as miRNAs, thereby limiting their effect in a clinical setting.26
As a traditional animal-based Chinese medicine, Gecko has clear effects in inhibiting tumor proliferation and metastases and promoting tumor apoptosis, especially in digestive system tumors such as esophageal, liver and gastric cancers.27 While current research on Gecko is mainly focused on herbal prescriptions and the utility of proteins, there is a lack of studies focusing on the extraction of small molecules such as miRNA. Furthermore, there are few studies on the mechanism of Gecko as an ingestible exogenous miRNA heterologous regulator of anti-tumor effects in humans. Therefore, exploring the regulatory network of Gecko-related miRNAs will aid our comprehensive understanding of the molecular mechanisms of gene interactions and regulation in human HCC, response to complex HCC and treatment through the multi-targeting effect of exogenous miRNAs in cells. It will also screen effective novel anti-tumor miRNA components, find alternative ways to synthesize the active ingredients of animal drugs, and provide new ideas for the inhibition mechanism of HCC. Finally, discovering circulating prognostic markers and predictive miRNAs in conjunction with efficacious drugs such as sorafenib is a future direction in the exploration of TCM.28, 29
In this investigation, we discovered that Gecko miRNA components such as miR-100-5p, miR-99a-5p and miR-101-3p are associated with HCC prognosis. These miRNAs may serve as biomarkers for the anti-HCC actions of Gecko. By targeting Polo-like kinase 1 (PLK1), miR-100 overexpression slows the growth of HCC and causes apoptosis.30 Moreover, by preventing lamellipodia and matrix metalloproteinase 2 (MMP2) activation, it also prevents HCC metastasis.31 Additionally, miR-100 regulates mTOR and insulin-like growth factor-1 receptor (IGF-1R) to promote autophagy and prevent the growth of HCC.32 MiR-99a, the 6th most abundant microRNA in normal human liver, is markedly reduced in HCC. Cholesterol-bound miR-99a mimics efficiently suppressed tumor growth and reduced α-fetoprotein levels in HCC-bearing mice upon intra-tumor injection.33 Moreover, miR-99a overexpression inhibited Argonaute 2 (AGO2) upregulation in HCC, suggesting that targeting AGO2 could significantly inhibit HCC growth in vivo.34 Furthermore, miR-101 functions as a liver tumor suppressor by regulating abnormal Nemo-like kinase (NLK) activity.35 Additionally, miR-101 overexpression inhibits proliferation, invasion, colony formation, and HCC progression in humans by directly targeting the EZH2 oncogene.36 The identified Gecko miRNA regulatory network is reliable, and potentially facilitates the regulation of HCC proliferation, metastasis and invasion, with therapeutic implications. Further exploration is needed as the focus on these 3 networks remains limited.
Although 9 critical prognostic downstream mRNAs regulated by miR-100-5p, miR-99a-5p and miR-101-3p were identified in HCC, only KPNA2, STMN1, EZH2, and LMNB1 were further related to progression as well as immunity in HCC. For instance, KPNA2 knockdown has been shown to reduce HCC cell migration and proliferation, while its overexpression increased malignancy of hepatoma carcinoma cell.37 The expression of STMN1 is significantly correlated with E2F1/TFPD1 and KPNA2 expression, as well as poor prognosis in HCC patients.38 EZH2, an epigenetic modifier, reduces the expression of the immunological checkpoint inhibitor PD-L1.39 Furthermore, LMNB1 promotes HCC progression by regulating the PI3K and MAPK signaling pathways.40 We constructed a Gecko miRNA network that is enriched in important tumorigenic pathways, such as PI3K-Akt41 and MAPK signaling pathways.42, 43 Additionally, it is associated with HCC, gallbladder cancer and other closely related oncological diseases. In addition, we observed a significant association between the above core targets of the Gecko miRNA regulatory network and immune cell infiltration, including B cells, T cells, macrophages, NK cells, and monocytes. The heterogeneous nature of immune cells in the tumor microenvironment compared to normal tissues highlights their importance in improving patient survival and the efficacy of immunotherapy for HCC.44 These reports improve the credibility of our study and further indicate the value of studying the effects of novel Gecko miRNA regulatory networks on HCC and immunotherapy.
Limitations
There is currently a lack of in vivo experimental evidence supporting the role of Gecko miRNAs. Furthermore, it is necessary to conduct further in-depth studies to determine whether Gecko have novel non-human miRNAs that also have anti-tumor effects.
Conclusions
We have constructed a miRNA-mRNA regulatory network of Gecko in HCC. This network exhibits significant correlations with patient prognosis and immune regulation in the tumor microenvironment. Our findings have opened up new insights into the use of herbal medicine to improve targeted therapy and immunotherapy for human HCC. In addition, this study provides a new direction for understanding the subsequent mechanism of tumor suppression by small molecules of TCM.
Supplementary data
The Supplemenary materials are available at https://doi.org/10.5281/zenodo.10624641. The package includes the following files:
Read_me_data.txt. A brief description of the database’s content sources.
Supplementary Fig. 1. Linearity assumption testing for univariate Cox regression models.
Supplementary Table 1. Proportional hazards hypothesis testing for univariate Cox analysis.
Supplementary Table 2. 15 target genes were found to be significantly upregulated in the regulatory network.
Supplementary Table 3. Downstream 38 target genes identified by univariate Cox regression analysis model were associated with prognosis.
Supplementary Table 4. The ranked data after transformation of sample expression.
Supplementary Table 5. The results of Spearman’s rank correlation coefficient.
Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.
Consent for publication
Not applicable.
