3.1 Several oncomiRCs are dysregulated in HCC

We initially compared the miRNA expression between RNAs from control and HepG2 or Huh7 HCC cells to identify the differentially expressed miRNAs in HCC. The nCounter miRNA assay revealed 81 DEmiRNAs (upregulated-27 and downregulated-54) in HepG2 cells. Further, among the upregulated miRNAs, 12 clustered miRNAs (miR-17-5p, miR-18a-5p-miR-17/92-a cluster, miR-106b-5p, miR-93-5p, miR-25-3p-miR-106b/25 cluster, miR-660-5p-miR-188/660 cluster, miR-103a-3p-miR-103a/103b cluster, miR-421-miR-374b/421 cluster, miR-7-5p-miR-7-2/3529 cluster, miR-182-5p-miR-96/183 cluster, miR-106a-5p-miR-18b/363 cluster, and miR-221-3p-miR-221/222 cluster) were upregulated in HCC (Fig. 2a). Furthermore, miR-106b-5p, miR-25-3p belonging to miR-106b/25 cluster was also upregulated in Huh7 cells (Supplementary Table 1). Additionally, small RNA sequencing data revealed that these miRNAs were upregulated in HCC samples compared with adjacent normal liver samples in the TCGA-LIHC cohort (Fig. 2b). While we identified several downregulated clustered miRNAs, this study aimed only to study oncomiRCs in detail. The complete list of all the DEmiRNAs and oncomiRCs identified in this study are provided in Tables 1 and 2.

Fig. 2

Differential expression of oncogenic miRNA clusters in HCC. a Heatmap of differentially expressed miRNAs showing the panel of 12 miRNAs upregulated in HepG2 cells compared with control RNA (Note: The n-Counter probes for miR-17-5p and miR-106b-5p are of equal sensitivity due to shared sequence design towards their respective targets). b Differential expression of the 12 miRNAs mapped to their respective miRNA clusters, including the miR-106b/25 cluster, in TCGA-LIHC compared with adjacent normal tissues

As all the constituent miRNAs of miR-106b/25 oncomiRC were differentially upregulated, we decided to characterize the gene targets and evaluate their clinical significance in detail.

3.2 miR-106b/25 oncomiRC targets multiple genes that act as tumor suppressors

We cultured HCC spheroids using HepG2 and Huh7 cells by combining the hanging drop and forced suspension methods (Fig. 3a). We observed that HepG2 cells, which are a relatively less aggressive subtype of HCC, conveniently formed self-assembled spheroids. However, such tumoroids could not be derived using more aggressive Huh7 cells, which instead formed loose aggregates (data not shown). This prompted us to characterize the cluster targets using HepG2 spheroids. We observed that spheroids derived using HepG2 cells reached maximum growth on day 5 at an optimum seeding density of 5k cells (Fig. 3b–d).

Fig. 3

Generation and molecular characterization of HCC spheroids. a Diagrammatic representation of the methodology used for spheroid generation. b Estimation of seeding density for spheroid culture. c Estimation of the number of days at different seeding densities for spheroid generation. d Representative image of HepG2-derived HCC spheroids. e Volcano plot showing all the differentially upregulated and downregulated genes in HCC spheroids compared with normal liver samples, highlighting established markers of HCC, such as overexpressed AFP and underexpressed CRP levels in spheroids

Whole-transcriptome sequencing of pools of 96 spheroids in biological duplicates revealed a total of 2902 upregulated and 3529 downregulated genes in spheroids compared with those in normal liver. Furthermore, known markers of HCC, such as elevated AFP and reduced CRP, have also been reported in HCC spheroids (Fig. 3e). Next, only the DEGs whose expression was negatively correlated with ocomiRC expression were considered for target gene identification. The target genes of oncomiRCs were identified by computing the overlaps between the 3529-underexpressed DEGs identified in spheroids, 3942 downregulated genes from the TCGA-LIHC cohort, 1148-predicted targets of the miR-106b/25 cluster, and 1217-potential coding tumor suppressor genes, which ultimately yielded a total of 14 differentially downregulated target genes of the miR-106b/25 oncomiRC which included-Caveolin 1 (CAV1), DNAJ Heat shock protein family B4 (DNAJB4), Protein Tyrosine Phosphatase Receptor type Delta (PTPRD), Major Facilitator Superfamily Domain-containing protein 2 A (MFSD2A), Transcription Factor 4 (TCF4), Kruppel-like Factor 6 (KLF6), Mutated in Colorectal Cancer or Colorectal Mutant Cancer protein (MCC), Cytochrome B5A (CYB5A), Estrogen Receptor 2 (ESR2), Nuclear Receptor subfamily 4 group A member 3 (NR4A3), Protein Kinase C Beta (PRKCB), Ras Association Domain Family Member 2 (RASSF2), Thioredoxin-interacting Protein (TXNIP), and Superoxide Dismutase 2 (SOD2) (Fig. 4a and b). Clustering analysis of these genes via unsupervised k-means clustering identified at least six different 6 clusters based on their expression patterns. We observed Cluster 1-CAV1, DNAJB4, and PTPRD; Cluster 2- MFSD2A, TCF4, and KLF6; Cluster 3-MCC, CYB5A; Cluster 4-ESR2, NR4A3, PRKCB, and RASSF2; Cluster 5-TXNIP; and Cluster 6-SOD2 as the different clusters on the basis of their differential expression (Fig. 4c).

Fig. 4

Identification of miR-106b/25 cluster targets in HCC. a Venn diagram showing 14 common gene targets of the miRNA cluster derived from overlapping gene sets from different databases and experimental whole-transcriptomic analysis. b Protein-protein interaction network for the 14 target genes of the miR-106b/25 cluster. c Heatmap of the different miR-106b/25 cluster gene targets clustered into six gene clusters on the basis of differential expression patterns

Furthermore, our transcriptomic analysis revealed that the lncRNA SOD2-overlapping transcript 1 (SOD2-OT1), like the miR-106b/25 cluster, was overexpressed in HepG2 tumuroids (Supplementary Fig. 1a–1c). Furthermore, we predicted a cis-interaction and a negative correlation between SOD2-OT1 lncRNA and SOD2 mRNA, suggesting that both the SOD2-OT1 and the miR-106b/25 cluster can target SOD2 mRNA. The mean binding energies obtained from the BiBiServe RNA Hybrid 2 tool for each of these RNA interactions are given in Supplementary Fig. 1d–1f.

Table 1 Differentially upregulated and downregulated miRNAs in HepG2-HCC cells3.3 miR-106b/25 oncomiRC target genes regulate critical gene ontologies and cell cycle pathways associated with HCC

We then queried the miR-106b/25 oncomiRC targets for gene signature (GSA) and enrichment (GEA). Compared with those in normal tissues, the expression of these genes was consistently downregulated in primary HCC and metastatic tissues (Fig. 5a). Furthermore, GEA revealed that the genes were involved in critical biological processes, such as the regulation of peptidyl serine phosphorylation, the regulation of fatty acid metabolism, the regulation of cellular ketone metabolism, carbohydrate transport, enzyme-linked receptor protein signaling, the hormonal response, the positive regulation of cellular metabolism and biosynthesis, and the negative regulation of cellular signal transduction (Fig. 5b). The key molecular functions included nuclear estrogen receptor activity, superoxide dismutase activity, oxidoreductase activity, TFIIB-class transcription factor binding, nuclear steroid receptor binding and activity, ion channel regulator activity, and Cis-regulatory region-specific DNA binding by RNA pol II (Fig. 5c).

Fig. 5

Gene signature and enrichment analysis. a Violin plots showing the signature expression of the 14 gene targets of the miR-106b/25 cluster across metastatic, tumor, and normal liver tissues in independent RNA-Seq and microarray data using TNMplot. Enriched gene ontologies in terms of the top 15 b Biological processes, c Molecular functions and d Key pathways associated with the fourteen miR-106b/25 cluster gene targets

Furthermore, the key pathways regulated by the cluster targets included several cancer-related pathways, such as Wnt signaling, TGF-beta signaling, VEGFA-VEGFR2 signaling, and EGF/R signaling pathways. Other pathways included fetal androgen synthesis, arsenic metabolism and ROS generation, serotonin and anxiety-related events, development of pulmonary dendritic cells and macrophage subpopulations, development of mammary glands during pregnancy and lactation, congenital lipodystrophy, the CRH signaling pathway, male infertility, and focal adhesion (Fig. 5d).

3.4 The miR-106b/25 cluster and its interactome are associated with the overall survival of HCC patients

We performed multivariate Cox regression analysis for the miRNAs of miR-106b/25 oncomiRC and its 14 gene targets. Notably, miR-93-5p was significantly associated with OS (p = 0.0246) (Fig. 6a). Moreover, HCC patients with higher expression of miR-93-5p had poorer OS than HCC patients with lower expression of this miRNA [HR = 0.72, 95% CI (0.54–0.96)]. Similarly, TCF4 [HR = 0.66, 95% CI (0.49–0.91), p = 0.0106], DNAJB4 [HR = 1.29, 95% CI (1.04–1.61), p = 0.0214], MCC [HR = 1.35, 95% CI (1.03–1.75) p = 0.0268], and CYB5A [HR = 0.77, 95% CI 0.59–0.99), p = 0.0423] were found to affect the OS of HCC patients (Fig. 6b).

Fig. 6

Analysis of the survival significance of the miR-106b/25 cluster and gene targets in HCC. a miRNAs of the miR-106b/25 cluster are significantly associated with overall survival in patients with HCC. b Specific gene targets of the miR-106b/25 cluster that are significantly associated with overall survival in HCC (highlighted in red)

Furthermore, we constructed a combined prognostic model for the miR-106b/25 cluster and its 14 gene targets to assess its actual clinical utility. This panel of 17 genes (3 oncogenic miRNAs and 14 tumor suppressor targets) could classify patients into high-risk and low-risk groups and predict their OS (p < 0.0001), with a high overall sensitivity of 90% and specificity of 94%, respectively (Fig. 7a and b). The top five genes affecting this model included PRKCB, miR-106b-5p, NR4A3, MFSD2A, and KLF6 (Fig. 7c).

Table 2 Oncogenic miRNA clusters, their component miRNAs, and their expression in HepG2-HCC cellsFig. 7

Construction of a combined prognostic model for the miR-106b/25 cluster and the 14 target genes in HCC. a Combined survival analysis and b distribution plots for the miR-106b/25 cluster and its interactome. c Bar plot showing the contribution of prominent miRNAs or target genes of the cluster that affect overall survival the most

3.5 Drug‒gene interaction analyses identified known drugs for repurposing in HCC

PanDrug analysis of 14 differentially expressed gene targets of miR-106b/25 cluster identified 69 drugs, 44 of which were approved for various cancers and 25 drugs were already at different stages of clinical trials for cancer treatment (Fig. 8a and b). The major drug families identified are estrogen receptor agonists, selective modulators of estrogen receptors, protein kinase inhibitors, estrogen receptor antagonists, CDK inhibitors, aromatase inhibitors, inhibitors of tubulin polymerization, PI3K inhibitors, KIT inhibitors, and DNA alkylating agents (Fig. 8c). Furthermore, drug-gene interaction network analysis revealed different approved drugs that can be repurposed for targeting HCC (Fig. 8d). The complete details of the PanCancer drug-gene interaction along with the interaction scores is provided in Supplementary Table 2.

Fig. 8

Pancancer drug-gene interaction analysis. a Representative plot showing the best available drugs to target 14 genes of the miR-106b/25 cluster. b Pie chart showing the experimental, approved, or under-trial status of different drugs in different cancers. c Distribution of families of drugs targeting the 14 gene targets of the miR-106b/25 cluster. d STITCH drug-gene interaction network with approved status