Introduction
Population-based cancer incidence data from the Global Cancer Observatory indicate a rising trend in liver cancer incidence.1 Globally, 4.7% of the 19.3 million new cancers and 8.3% of cancer-related mortalities are associated with liver cancer. 2 Nearly half of all liver cancer cases and liver cancer deaths worldwide occur in China (45.3% and 47.1%, respectively). 3 The epidemiology of liver cancer is dominated by hepatocellular carcinoma (HCC), which accounts for about 70% of primary liver cancers. 4 HCC is one of the most fatal malignancies, usually diagnosed at advanced stages with very limited therapeutic options, leading to a dismal prognosis. 5 Current survival rates for liver cancers remain low (21%). 1 The therapeutic outcomes for patients with hepatobiliary cancer are unsatisfactory due to the bleak prognosis, leaving significant room for improvement.
The etiologies and geographical trends of HCC are complex. Molecular heterogeneity contributes to the complexity of HCC at different levels. 6 Recently, next-generation sequencing-based technologies have been widely used for multiomic and single-cell analyses. These advances have enabled the creation of high-resolution atlases of tumor genomes and epigenomes, accelerating our understanding of tumor spatial and temporal heterogeneity. 7 HCC involves a multitude of genetic and epigenetic aberrations. 8 A deeper understanding of cancer epigenetics in tissue samples and biopsies is increasingly valuable in predicting patient prognosis and treatment response in clinical settings. Epigenetic mechanisms have revealed that alterations in ncRNAs, such as miRNAs and circRNAs, promote uncontrolled cell growth, metastasis, and HCC progression. 9,10 These ncRNAs have also emerged as a critical area of research in HCC-targeting therapies. 11
In the present research, bioinformatics analyses were used to identify differentially expressed circRNAs, miRNAs, and HCC-related genes in HCC. The target miRNAs of the circRNAs and mRNAs of the miRNAs were predicted and intersected with the differentially expressed miRNAs and HCC-related genes to identify common targets. A network was subsequently established, and the hsa_circ_0001726/miR-140-3p/KRAS axis was further investigated and validated.
Methods
RNA expression data retrieval and processing
RNA-seq RNA expression data for HCC and adjacent noncancerous liver samples were obtained from the Gene Expression Omnibus (GEO). The GSE235991 dataset contains five high-invasive and five low-invasive HCC tissues for circRNA expression profiling. The GSE121714 and GSE216115 datasets provide circRNA microarray analysis data for HCC and nontumor tissues. GSE82303, a dataset of blood exosomal miRNA profiles, includes two subsets, GSE82301 and GSE82302, with 72 patients with normal livers and 32 patients with HCC.
Construction of the circRNA-miRNA-gene network
The GSE235991, GSE121714, and GSE216115 datasets were analyzed for differentially expressed circRNAs in HCC using the GEO2R online tool (logFC > 1; p < 0.05). The common circRNAs among these datasets were identified using a Venn diagram. Potential downstream miRNAs were collected from ENCORI, the Circular RNA Interactome, and the Cancer-Specific CircRNA database. The GSE82303 dataset was analyzed for differentially expressed miRNAs in HCC using GEO2R (logFC > 1; p < 0.05). The downstream miRNAs were then intersected with the differentially expressed miRNAs in HCC using VennPainter to obtain the HCC-related, circRNA-targeting miRNAs. The miRNA targets were obtained from TargetScan and the Encyclopedia of RNA Interactomes (ENCORI/starBase). HCC-associated genes were retrieved from GeneCards using the keyword “hepatocellular carcinoma” and filtered by a twofold mean relevance score. The overlapping genes from miRNA-targeting genes and HCC-associated genes were used to construct the circRNA-miRNA-gene network, which was visualized using Cytoscape software.
Gene Ontology (GO) and pathway enrichment analysis
The genes involved were subjected to enrichment analysis of GO terms for biological process, cellular component, molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using OmicShare (https://www.omicshare.com/tools/ ).
Patient material
HCC patients included in this study (n = 133) were enrolled at the 900TH Hospital of the Joint Logistics Support Force. These patients had no history of cancer-related therapy or other serious diseases. All patients provided written informed consent, and the study was approved by the institutional ethical review board. Patient information, including characteristics and clinicopathological variables, is described in Table 1, following current clinical definitions in China.
Table 1Association of hsa_circ_0001726 expression with clinicopathological features in hepatocellular carcinoma (n = 133)
| Characteristics | Level | Low expression of hsa_circ_0001726 | High expression of hsa_circ_0001726 | p | Test |
|---|
| T stage | T1&T2 | 43 (32.33) | 31 (23.31) | 0.017 | Chi-Square |
| T3&T4 | 22 (16.54) | 37 (27.82) | | |
| N stage | N0 | 63 (47.37) | 63 (47.37) | 0.441 | Fisher’s Exact |
| N1 | 2 (1.50) | 5 (3.76) | | |
| M stage | M0 | 63 (47.37) | 65 (48.87) | 1.000 | Fisher’s Exact |
| M1 | 2 (1.50) | 3 (2.26) | | |
| BCLC stage | Stage I & Stage II | 50 (37.59) | 44 (33.08) | 0.122 | Chi-Square |
| Stage III & Stage IV | 15 (11.28) | 24 (18.05) | | |
| Gender | Female | 28 (21.05) | 25 (18.80) | 0.457 | Chi-Square |
| Male | 37 (27.82) | 43 (32.33) | | |
| Age | ≤61 | 34 (25.56) | 38 (28.57) | 0.679 | Chi-Square |
| >61 | 31 (23.31) | 30 (22.56) | | |
| BMI | ≤24 | 29 (21.80) | 35 (26.32) | 0.429 | |
| >24 | 36 (27.07) | 33 (24.81) | | |
| Histologic grade | G1 & G2 | 46 (34.59) | 41 (30.83) | 0.204 | Chi-Square |
| G3 & G4 | 19 (14.29) | 27 (20.30) | | |
| AFP | ≤356 ng/ml | 54 (40.60) | 53 (39.85) | 0.455 | Chi-Square |
| >356 ng/ml | 11 (8.27) | 15 (11.28) | | |
| Fibrosis Ishak score (vs.) | 0/1/2 | 49 (36.84) | 52 (39.10) | 0.884 | Chi-Square |
| 3/4/5/6 | 16 (12.03) | 16 (12.03) | | |
| OS event, n (%) | Alive | 55 (41.35) | 41 (30.83) | 0.002 | Chi-Square |
| Dead | 10 (7.52) | 27 (20.30) | | |
Cell culture
HCC cells were cultured in specific media (Gibco, Carlsbad, CA, USA) as recommended by the supplier (China Center for Type Culture Collection). For transfection experiments, the cells were preincubated with the indicated concentrations of inhibitors, mimics, or siRNAs for 24 h.
Establishment of lenvatinib-resistant cell lines
To determine the half-maximal inhibitory concentration (IC50) of lenvatinib, MHCC97-H, and Hep 3B2.1-7 cells were incubated with various concentrations of lenvatinib, and cell viability was measured after three days. For the establishment of lenvatinib-resistant cell lines, MHCC97-H and Hep 3B2.1-7 cells were incubated with lenvatinib at the approximate IC50 concentration, which was increased by 0.2 µmol/L per week for 6 months. These cells were termed MHCC97-H-LR and Hep 3B2.1-7-LR, respectively.
RNA isolation and real-time reverse transcription-quantitative polymerase chain reaction (RT-qPCR)
Fresh tumor samples were collected and preserved in RNAlater (Qiagen, Hilden, Germany). Total RNA was extracted from tissues and cells using the RNeasy kit (Qiagen, Valencia, USA), and RNA concentrations were measured using an ND-100 Nanodrop spectrophotometer. cDNAs were obtained through reverse transcription using the PrimeScript™ 1st Strand cDNA Synthesis Kit (TaKaRa Bio, Japan). RT-qPCR was performed using RT-Taq/SYBR Green QPCR reagents (Invitrogen, Carlsbad, USA) on a Light-Cycler96 thermocycler (Roche, Basel, Switzerland). Relative expression values were calculated using the 2–ΔΔCT method.
Dual-luciferase reporter assay
The wt-circ plasmid contained the synthesized 3′-untranslated region (UTR) wild-type sequence of hsa_circ_0001726, while the mut-circ plasmid contained the mutant sequence, along with the Fluc coding sequence. The WT KRAS luciferase reporter plasmid was generated using the full-length KRAS 3′ UTR vector, whereas the MUT KRAS plasmid was generated with the mutant KRAS 3′ UTR vector. The cells were subjected to lysis in Passive Lysis Buffer (Progma, USA), and luciferase activity was measured following the Dual-Luciferase® Reporter 1000 Assay System protocol (Promega).
Cell proliferation assay
The proliferation of transfected or treated cells was evaluated using the BeyoClick EdU Cell Proliferation Kit (APExBIO, Houston, USA). Cells were seeded into 96-well plates at 40% confluence. The kit reagent was added to each well on different days (zero, one, two, three, and four days), followed by measurement of absorption at 450 nm two hours later.
In vitro Transwell migration
HCC cells were subjected to an assessment of in vitro Transwell migratory ability using HTS Transwell-96 Permeable Support with an 8.0 µm pore polyester membrane (Corning Life Sciences - Axygen, Union City, USA). Cells were pretreated with serum-free medium and added to the top wells of Transwell inserts. The bottom wells were filled with 600 µL of media containing 10% fetal bovine serum (Gibco, Grand Island, USA). After 24 h of incubation, cells in the bottom wells were harvested and counted.
Statistical analysis
For group comparisons, continuous variables were compared using the Mann-Whitney U test, while categorical variables were compared using the chi-square test or Fisher’s exact test. Logistic regression analysis was used for Kaplan-Meier analysis with categorical variables. Univariate and multivariate Cox regression analyses were performed to evaluate the odds ratios for associations between survival risk factors and events. Statistical significance was set at a two-sided p-value < 0.05.
Results
Differentially expressed circRNAs in HCC based on microarray data
We retrieved the GSE235991, GSE121714, and GSE216115 datasets to identify differentially expressed circRNAs in HCC. A total of 328 differentially expressed circRNAs were identified from the GSE235991 dataset, 252 from the GSE121714 dataset, and 225 from the GSE216115 dataset. Among these differentially expressed circRNAs, hsa_circ_0001726, hsa_circ_0008558, hsa_circ_0005354, hsa_circ_0035436, hsa_circ_0035435, and hsa_circ_0006633 were common (Fig. 1A). The detailed expression levels of these six circRNAs are depicted in Figure 1B–D.
Differentially expressed miRNAs in HCC based on microarray data and circRNA target prediction
The GSE82303 dataset was retrieved and analyzed to identify 100 differentially expressed miRNAs in HCC. The prediction of downstream miRNAs was conducted via databases. Hsa-miR-1182, hsa-miR-1208, hsa-miR-1290, hsa-miR-1299, hsa-miR-140-3p, hsa-miR-494-3p, and hsa-miR-622 were identified as hsa_circ_0001726-targeting differentially expressed miRNAs in HCC, while hsa-miR-623 and hsa-miR-659-3p were identified as hsa_circ_0005354-targeting miRNAs in HCC (Fig. 2A). Hsa-miR-623 and hsa-miR-659-3p were also identified as hsa_circ_0006633-targeting miRNAs in HCC, and hsa-miR-1287-5p and hsa-miR-584-5p were identified as hsa_circ_0008558-targeting miRNAs in HCC (Fig. 2B). Three miRNAs—hsa-miR-557, hsa-miR-526b-5p, and hsa-miR-622—were observed to target differentially expressed miRNAs in HCC, while hsa-miR-623 was identified as the hsa_circ_0035436-targeting miRNA in HCC (Fig. 2C).
The circRNA-miRNA-mRNA regulatory network involved in HCC
HCC-related genes were obtained from GeneCards and intersected with the miRNA target genes, resulting in a total of 99 genes. After protein-protein interaction prediction, disconnected nodes in the network were hidden, and the interaction network was analyzed via Cytoscape software, which revealed 82 nodes and 426 edges (Fig. 3A). The first six biological process terms and the first six GO terms enriched by the genes included GO:0008283 cell population proliferation, GO:0042127 regulation of cell population proliferation, GO:0012501 programmed cell death, GO:0006915 apoptotic process, GO:0042981 regulation of the apoptotic process, and GO:0043067 regulation of programmed cell death (Fig. 3B). Pathways in cancer, apoptosis-multiple species, microRNAs in cancer, the ErbB signaling pathway, and the hepatocellular carcinoma pathway were among the top 20 KEGG pathways (Fig. 3C). Finally, the circRNA-miRNA-gene regulatory network involved in HCC was constructed (Fig. 3D), consisting of 34 circRNA-miRNA pairs and 194 miRNA-mRNA pairs.
Validation of hsa_circ_0001726, miR-140-3p, and KRAS expression in HCC
The microarray dataset analyses revealed the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC. Furthermore, using RT-qPCR, the expression of hsa_circ_0001726, miR-140-3p, and KRAS was evaluated in 133 HCC and matched adjacent normal tissues. We confirmed that hsa_circ_0001726 was significantly upregulated in HCC (Fig. 4A). miR-140-3p was downregulated in HCC (Fig. 4B), while KRAS mRNA was upregulated (Fig. 4C). Similar alterations were observed in HCC cells (Fig. 4D, F).
Validation of the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC
To further validate the relationships between hsa_circ_0001726 and miR-140-3p, as well as between miR-140-3p and KRAS, we performed dual-luciferase reporter assays and RNA pull-down analyses. Prediction revealed binding sites in hsa_circ_0001726 and miR-140-3p (Fig. 5A). Notably, hsa_circ_0001726 was negatively correlated with miR-140-3p (Fig. 5B). A pull-down experiment in MHCC-97-H cells revealed that hsa_circ_0001726 was more enriched in pull-down experiments involving miR-140-3p than in the negative control (Fig. 5C). A dual-luciferase reporter assay in Hep 3B2.1-7 cells confirmed these findings, revealing enhanced luciferase activity in miR-140-3p-inhibited cells and reduced activity in miR-140-3p-overexpressing cells (Fig. 5D). We also examined the relationship between miR-140-3p and KRAS (Fig. 5E–H).
Role of the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC progression
To investigate the role of the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC progression, we analyzed HCC cell proliferation and migration after modulating the axis. KRAS mRNA expression was significantly downregulated after hsa_circ_0001726 inhibition, but this effect was rescued by the miR-140-3p inhibitor and reversed by KRAS knockdown (Fig. 6A, B). Cell proliferation assays revealed that hsa_circ_0001726 knockdown inhibited cell growth, which was partly reversed by the miR-140-3p inhibitor and restored by KRAS knockdown (Fig. 6C, D). Similarly, hsa_circ_0001726 and KRAS knockdown inhibited HCC cell migration, while the miR-140-3p inhibitor promoted HCC cell migration (Fig. 6E, F). We next assessed the clinical role of hsa_circ_0001726 in the prognosis of HCC patients. hsa_circ_0001726 expression was significantly correlated with T stage and overall survival (Table 1). Cox regression analysis was used to calculate the risk of poor prognosis (hazard–risk ratio) for the high and low hsa_circ_0001726 populations. There was a 2.900-fold (95% CI: 1.402–5.998) increased risk of disease recurrence or death in high-hsa_circ_0001726-expressing patients compared with low-hsa_circ_0001726-expressing patients after univariate Cox proportional hazards analysis and a 2.402-fold (95% CI: 1.152–5.009) increased risk after multivariate analysis (Table 2). The cumulative survival of the enrolled HCC patients was estimated via Kaplan-Meier analyses, which revealed that HCC patients with high hsa_circ_0001726 levels had a lower mean survival rate than those with low hsa_circ_0001726 levels (Fig. 6G).
Table 2A multivariate and univariate Cox proportional hazards analysis of hsa_circ_0001726 expression in hepatocellular carcinoma overall survival
| Characteristics | Univariate analysis
| Multivariate analysis
|
|---|
| Hazard ratio (95% CI) | p-value | Hazard ratio (95% CI) | p-value |
|---|
| Hsa_circ_0001726 | 2.900 (1.402–5.998) | 0.004 | 2.402 (1.152–5.009) | 0.019 |
| T stage (T3 & T4 vs. T1&T2) | 2.507 (1.275–4.9300 | 0.008 | 1.302 (0.547–3.101) | 0.551 |
| N stage (N1 vs. N0) | 2.220 (0.679–7.256) | 0.187 | | |
| M stage (M1 vs. M0) | 1.536 (0.367–6.424) | 0.556 | | |
| BCLC stage (Stage III & Stage IV vs Stage I & Stage II) | 2.435 (1.270–4.668) | 0.007 | 1.793 (0.787–4.081) | 0.164 |
| Gender (Male vs Female) | 1.169 (0.609–2.241) | 0.639 | | |
| Age (>61 vs ≤61) | 1.173 (0.611–2.250) | 0.631 | | |
| BMI (≤24 vs >24) | 1.517 (0.786–2.926) | 0.214 | | |
| Histologic grade (G3 & G4 vs G1 & G2) | 2.018 (1.055–3.860) | 0.034 | 1.756 (0.895–3.446) | 0.102 |
| AFP (>356 ng/ml vs ≤356 ng/ml) | 2.580 (1.274–5.224) | 0.008 | 2.247 (1.099–4.595) | 0.027 |
| Fibrosis Ishak score (3/4/5/6 vs. 0/1/2) | 1.464 (0.734–2.920) | 0.279 | | |
Role of the hsa_circ_0001726/miR-140-3p/KRAS axis in lenvatinib resistance
To examine whether the hsa_circ_0001726/miR-140-3p/KRAS axis influences lenvatinib resistance, lenvatinib-resistant Hep 3B2.1–7 and MHCC97-H cells were established. Hsa_circ_0001726 was upregulated in lenvatinib-resistant HCC cells, particularly in lenvatinib-resistant MHCC97-H cells (Fig. 7A). After treatment with varying concentrations of lenvatinib, hsa_circ_0001726-knockdown MHCC97-H cells exhibited a lower IC50 value for lenvatinib compared with negative controls (Fig. 7B). As shown in Figure 7C, compared with negative control cells, hsa_circ_0001726-silenced MHCC97-H cells were more sensitive to lenvatinib and presented a decreased relative resistance factor. KRAS mRNA expression was subsequently examined in transfected MHCC97-H cells via RT-qPCR, which confirmed that lenvatinib-resistant MHCC97-H cells were successfully transfected (Fig. 7D). The viability of hsa_circ_0001726-knockdown cells resistant to lenvatinib was lower than that of negative control-transfected cells resistant to lenvatinib, but the miR-140-3p inhibitor and KRAS siRNA reversed this effect (Fig. 7E). The migration of lenvatinib-resistant cells was measured via Transwell assay, revealing less migration in hsa_circ_0001726-knockdown and lenvatinib-resistant cells compared with negative control cells (Fig. 7F).
Discussion
Advances in sequencing technologies have led to the discovery of ncRNAs, and increasing research has revealed their regulatory role in cellular processes and pathways in cancer. 12 It is challenging to uncover ncRNA functions in isolation. CircRNAs regulate the abundance of specific miRNAs through mechanisms known as sequestration; miRNAs, in turn, target the mRNAs of various genes through complementary combinations; and vice versa, each gene can be targeted by multiple miRNAs. Thus, circRNAs link associated miRNAs and genes into a complex regulatory network. 13 The highly complex characteristics of circRNA interactions support their regulatory roles in important cellular processes. In our study, we analyzed recent datasets to identify differentially expressed circRNAs in HCC. We combined dataset analysis and in silico target prediction approaches to construct a circRNA-miRNA-gene network in HCC patients. Additionally, cell line experiments were conducted to validate the role of the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC. We revealed that hsa_circ_0001726 promotes HCC progression and predicts HCC prognosis.
The decoy role of circRNAs is based on their ability to sequester miRNAs. 14 CircRNAs containing complementary sequences for miRNAs can function as sponges for several miRNAs. The hypothesis of a ceRNA network suggests that circRNAs regulate the expression of other genes and influence tumorigenesis and cancer progression. 15 Based on comprehensive analysis, Zhang et al. established a prognostic ceRNA network of 21 circRNAs, five miRNAs, and seven hub mRNAs for hepatocellular carcinoma. 16 Our experimental results revealed an HCC-related network consisting of six circRNAs, 13 miRNAs, and 88 mRNAs. Hsa_circ_0001726 has been identified as a back-spliced junction and a prognostic marker in esophageal squamous cell carcinoma. 17 Hsa_circ_0008558 modulates the maturation and exosomal dissemination of miR-17, enhancing colorectal carcinoma metastasis. 18 Hsa_circ_0008558 was shown to be downregulated in HCC. 19 Hsa_circ_0005354 has been linked to the risk of acute myeloid leukemia in pediatric patients. 20 Hsa_circ_0035436 exhibited significant differences between clear cell renal cell carcinoma cell lines and 293T cell lines. 21 Hsa_circ_0006633 is expressed at low levels in HCC patients, is associated with poor overall survival, and can sponge miR-545-3p to promote Smad7 expression, thereby inhibiting HCC progression. 22
For a long time after their discovery, ncRNAs were considered nonfunctional molecules or “transcriptional noise.” However, in recent years, ncRNAs have been recognized for their crucial biological roles as transcription and translation regulators in various diseases, including cancer. 23 This recognition has prompted scientists to explore their potential in clinical practice, considering them new targets or strategies for cancer treatment. We found that hsa_circ_0001726 promotes HCC by accelerating cell proliferation and migration, and its knockdown can inhibit HCC progression. A clear correlation between hsa_circ_0001726 expression and HCC T stage was observed. The predictive power of hsa_circ_0001726 for overall survival in HCC patients was also demonstrated. These findings suggest that hsa_circ_0001726, along with the axis it forms (hsa_circ_0001726/miR-140-3p/KRAS), has potential as a target for HCC therapy.
RAS is a member of the small GTPase family, and its gain-of-function mutations are among the leading causes of human oncogenesis. 24 However, attempts to develop RAS-targeted therapies have often failed. Therefore, the development of RAS inhibitors targeting upstream regulatory molecules offers a novel approach to overcoming anti-RAS challenges. 24 The KRAS family is encoded by KRAS genes in humans. The RAF/MEK/ERK pathway is the primary pathway involved in the anticancer effects of some promising anti-HCC tyrosine kinase inhibitors, such as lenvatinib. 25 In our study, we found that KRAS silencing inhibits HCC cell proliferation and migration, consistent with previous results. 26 Furthermore, we validated the role of the hsa_circ_0001726/miR-140-3p/KRAS axis in lenvatinib resistance in HCC. Our results showed that downregulation of hsa_circ_0001726 increases miR-140-3p expression and decreases KRAS mRNA expression in lenvatinib-resistant HCC cells, thereby inhibiting HCC progression. These results support our hypothesis that hsa_circ_0001726 is involved in lenvatinib resistance in HCC. However, this study lacks more in-depth mechanistic investigations of the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC and lenvatinib resistance.
Conclusions
This study established an HCC-related circRNA-miRNA-mRNA regulatory network and identified hsa_circ_0001726 as a novel potential biomarker for HCC prognosis. Additionally, we experimentally revealed the involvement of the hsa_circ_0001726/miR-140-3p/KRAS axis in HCC progression and lenvatinib resistance.
Declarations
Ethical statement
All analyses were approved by the 900TH Hospital of Joint Logistics Support Force Institutional Ethical Review Board (No. 2023016) and were conducted in accordance with the Declaration of Helsinki. All patients provided written informed consent.
Data sharing statement
The datasets generated and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Funding
This work was supported by the Fujian Natural Science Foundation Project (General) [Grant number: 2021J011265] and the Fund for Scientific Research Projects of the Organ Transplantation Section of the Joint Logistics Support Force Key Discipline of Joint Logistics Medicine [Grant number: LQZD-QG].
Conflict of interest
The authors have no conflict of interests related to this publication.
Authors’ contributions
Literature research, study conception (ZC), protocol development, ethical approval, patient recruitment, data analysis (JY, YX, ZL), the first draft writing of the manuscript (XC), manuscript review and editing, and approval of the final version of the manuscript (YC, LL). All authors have approved the final version and publication of the manuscript.