Background

Hepatocellular carcinoma (HCC) is a common reason for cancer-associated mortality, characterized by its highly metastatic and invasive abilities.1 Most HCC cases are diagnosed when the disease has progressed to an advanced stage. Despite significant advancements in the treatment of HCC, the 5-year survival rate is approximately 18%.2 Increasing evidence suggests that the development of HCC is associated with hepatic injury and cirrhosis due to chronic inflammation, which enhances immune tolerance and leads to immune escape.3 Moreover, the liver plays a crucial role in immune regulation and is a natural organ for immune tolerance.4 Immune checkpoint inhibitors (ICIs) have revolutionized HCC treatment, with anti-PD-1/PD-L1 therapies showing promise, particularly in combination with other ICIs or anti-angiogenic drugs.5,6 These combinations may synergize to enhance efficacy, overcoming the limitations of single-agent therapies. Despite this potential, the response rate to ICIs in HCC remains suboptimal, necessitating the exploration of combination therapies with agents such as anti-angiogenics to improve outcomes.7 Predictive biomarkers for ICI efficacy in HCC are still under investigation, with factors such as etiology, tumor markers, and gene signatures showing promise in patient selection for therapy.

MicroRNA (miRNAs), small non-coding RNAs consisting of approximately 22 nucleotides, play a vital role in regulating mRNA expression. Previous research has demonstrated that miRNAs are involved in the development, progression, and drug resistance of cancers.8 For example, the loss of microRNA-622 results in combined de-suppression of the MAPK14-ATF2-axis, leading to chemoresistance in HCC.9,10 Additionally, as a pivotal contributor to the tumor microenvironment (TME), miRNA plays a key participant part in regulating regulate natural killer (NK) cells, cytotoxic T lymphocytes, myeloid-derived suppressor cells, dendritic cells (DCs), as well as regulatory T cells through intercellular communication.11 In colon cancer patients, miR-448 enhances the response of CD8+ T cells by suppressing the expression of indoleamine 2.3-dioxygenase 1.12 Additionally, miR-1468-5p in exosome targets homeobox containing 1 (HMBOX1) and activates the JAK2/STAT3 pathway suppressor of cytokine signaling 1 (SOCS1), ultimately leading to immunosuppression and allowing immune escape of cervical cancer cells.13 Meanwhile, the 3’-untranslated region (3’-UTR) of Beta-actin (ACTB) mRNA significantly contributes to HCC progression by enhancing the migratory and invasive capabilities of HCC cells. This is achieved through the upregulation of the miR-1 target gene MET and the miR-29a target gene MCL1, as well as through the modulation of miR-1 and miR-29a degradation mediated by the Argonaute 2 (AGO2) protein.14 However, there are only a few studies investigating the prognostic value of immune-related miRNAs and mRNA among HCC patients and the associated regulatory mechanisms. Therefore, this study aims to explore immune-associated miRNAs and mRNAs, their prognostic effects and potential molecular regulatory mechanisms.

This study identified immune-associated miRNAs that correlated with the prognosis of HCC. Meanwhile, corresponding miRNA and mRNA risk score models were established, which may enhance our understanding of oncogenic and immune-related molecular pathways in HCC. These findings may also aid in identifying prognostic biomarkers as well as therapeutic targets for selecting immunotherapies.

Materials and MethodsData Source

The sequencing information and clinical data for the TCGA-HCC cohort were acquired from the TCGA database (https://xenabrowser.net) in the April of 2022. MiRNA expression matrix contained 50 normal and 370 hCC samples, of which 364 hCC samples had clinical and survival data. In the mRNA and lncRNA expression matrices, a total of 50 normal samples and 374 hCC samples were collected, among which 368 hCC samples had clinical and survival data, while 6 samples lacking survival information were excluded. Two HCC datasets, GSE31384 and GSE76427, were acquired from the GEO database (https://www.ncbi.nlm.nih.gov/). GSE31384 is a miRNA dataset comprising 166 hCC samples with survival data, which was utilized for external validation of the miRNA risk score. The external validation for the mRNA risk score model was conducted by using GSE76427, which includes 115 hCC samples with survival information.15 The flow chart of this study was shown in Supplementary Figure 1.

Variance Expression Analysis

The differentially expressed miRNAs and mRNAs were analyzed within the TCGA-miRNA and TCGA-mRNA expression matrices, respectively, using the R package ‘DESeq2’ (version 1.38.0).16 The Benjamini & Hochberg method was used to correct for multiple testing and obtain the corrected p-value, ie, the adjusted p-value. Differentially expressed lncRNAs, mRNAs and miRNAs (DE-lncRNAs, DE-mRNAs, DE-miRNAs) were mined according to the adjusted p-value < 0.05 and |log2FoldChange(FC)| > 0.5. The immunological score of each sample was calculated in the TCGA dataset using the R package ‘estimate’ (version 1.1.5)17. The correlation between DE-miRNAs differential miRNAs and immune scores was calculated by Spearman, and the corresponding p values and correlation coefficients cor were obtained, and the immune-associated differential miRNAs were screened to obtain the immune-associated differential miRNAs according to the threshold value of pvalue<0.05and|cor|>0.3.

Functional Enrichment Analysis

Functional enrichment analysis of miRNAs was performed using Funrich software with a significance threshold set at p < 0.05. The DAVID tool (version 6.8, https://david-d.ncifcrf.gov/) was used to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. A p-value < 0.05 and a count ≥ 2 were predefined as the thresholds for enrichment significance. The analysis included cellular component (CC), molecular function (MF), as well as biological process (BP) in GO analysis. The visualization of the enrichment analysis outcomes were accomplished with the application of the “ggplot2” library within the R programming environment.

Establishment of Immune-Associated miRNA and mRNA Risk Score Model in HCC

For the miRNA model, the 364 hCC cases containing survival data from the TCGA-HCC dataset, served as the training set. Univariate Cox and stepwise multivariate analyses were deployed to select prognostic miRNAs in the training set. For the mRNA model, the 368 hCC cases with survival data from the TCGA-HCC dataset were used as the training set. Univariate Cox and LASSO regression analyses were applied to identify prognosis-associated genes in the training set. Set the P-value threshold for univariate Cox proportional hazards regression analysis at 0.05. The mRNA selected through univariate Cox regression analysis underwent feature mRNA screening using the Lasso algorithm in the training set TCGA-mRNA. The glmnet package in R (Version 4.0–2) was utilized to perform regression analysis on the mRNA from the univariate results, with the parameters set as follows: family = “cox”, nfold = 20, indicating 20-fold cross-validation. The feature mRNA was screened accordingly. For model accuracy, lambda.min (0.06051506) was selected as the λ for model construction, leading to the final selection of prognosis-associated genes. Using the formula (Riskscore = , where coef represents the coefficient obtained from stepwise multivariate analysis or LASSO. Based on the median of risk score, the samples were categorized into two distinct risk categories (low risk [LR] and high risk [HR]). The GSE31384 dataset (166 hCC samples) and GSE76427 (115 hCC samples) datasets were used for external validation of the miRNA and mRNA risk score models, respectively. Receiver Operating Characteristic (ROC) analysis was performed using the “survivalROC” package (version 1.0.3), and risk curves and Kaplan-Meier (K-M) curves were plotted using the “ggplot2” (version 3.3.5) and “survminer” (version 0.4.9) packages, respectively, to demonstrate the predictive capacity of the models.

Independent Prognostic Analysis

Pathological information and risk scores were included in univariate and multivariate Cox analysis to authenticate independent prognostic predictors.

Immune Cell Infiltration Analysis

The differences in the abundance of 28 kinds of infiltrating immune cells between two risk subgroups were assessed using the ssGSEA algorithm from the “GSVA” package (version 1.34.0),18 with statistical significance determined by the Wilcoxon test (p < 0.05). The ssGSEA algorithm analyzes the gene expression matrix and immune cell gene sets to calculate the relative abundance score (ssGSEA score) for specific immune cell types in each sample. The height of the score reflects the enrichment level of the immune cell type. By comparing the scores of different samples on the same gene set, the differences in immune cell types can be assessed, and by combining the scores from multiple gene sets, a comprehensive evaluation of the heterogeneity of immune cells in the sample can be achieved. Additionally, to explore immune infiltration in HCC, the abundance and differences in 22 specific infiltrating immune cells between high and low risk subgroups were evaluated using the Cibersort algorithm, with significant differences identified at p < 0.05. Using the R package “psych” (version 2.1.6), we analyzed the Spearman correlation between the abundance of differential immune-infiltrating cells and prognostic genes. This analysis calculated the correlation coefficients (r-values) and significance levels (p-values) between them. Specifically, a larger absolute value of r (|r|) indicates a higher correlation, with positive values indicating a positive correlation and negative values indicating a negative correlation. A p-value less than 0.05 suggests a significant correlation between the two variables. The results for both analyses were visually presented using box plots generated by the R package “ggplot2” (version 3.3.5).19

Therapy Analysis Based on Risk Score Model

The sensitivity of the two risk groups to immune checkpoint blockade therapy was predicted and evaluated on the Tumor Immune Dysfunction and Exclusion (TIDE) website.20 Via the “pRRophetic” package (version 0.2),21 we determined the IC50 values for 138 chemotherapeutic agents across two patient risk subgroups of patients, with the aim of evaluating their responsiveness to chemotherapy treatments. The R package “ggplot2” (version 3.3.5) was used to visualize the results and the Wilcoxon test was used to screen for drugs with significant differences between groups. Ranked these drugs by their p-values and selected the top 5 drugs to display in a boxplot.

The mRNA, lncRNA Prediction and Construction of Competing Endogenous RNA (ceRNA) Network

Predictions of miRNA-regulated mRNAs were conducted using the miRwalk 3.0 database with default parameters (binding site: 3’ UTR, binding probability≥0.95) and confirmed by miRDB database. The predicted mRNAs were then intersected with DE-mRNAs, and final selection of downstream mRNAs was made by identifying those that exhibited inverse expression patterns relative to their corresponding miRNAs. The lncRNAs regulating miRNAs were predicted using the LncBase V2.0 with a threshold score > 0.6. The predicted lncRNAs were then intersected with DE-lncRNAs, and the final upstream lncRNAs were obtained by considering the opposite expression trends of miRNAs and lncRNAs. The final lncRNA-miRNA-mRNA interaction network was depicted with the aid of Cytoscape platform (version 3.6.1).22

Protein Level Analysis of Prognosis-Related Genes

To understand the expression levels of prognosis-related genes in different tissues, the “Multiple Genes Query” function offered by the Human Protein Atlas (HPA) (https://www.proteinatlas.org) database was utilized to conduct tissue-specific enrichment analysis on these prognosis-related genes. By examining the Transcripts Per Million (TPM) values of the genes across various tissues, the expression levels of each prognosis-related genes in different tissue locations were assessed.23

The Verification of Prognosis-Related mRNA Expressions Within Cell Lines

Human normal hepatocytes, QSG-7701, and three human HCC cell lines (HepG2, LM3, and HuH-7) were obtained from iCell Bioscience Inc (Shanghai, China). Quantitative polymerase chain reaction (qPCR) was performed according to the manufacturer’s guidelines. The primer sequences utilized in our experiment are presented in Table 1. The expression levels were standardized relative to the internal reference (GAPDH) and calculated utilizing the 2−ΔΔCq method.24

Table 1 The Primer Sequences for qPCR

Statistical Analysis

R language was used for all the bioinformatics analyses in this study. We compared data in different groups by using Wilcoxon test. The discrepancy in RT-qPCR results were compared by employing the Student’s t-test. P value < 0.05 represents statistically significant difference.

ResultsImmune-Related miRNAs in HCC

A total of 246 DE-miRNAs were identified between HCC and non-tumor samples, containing 90 up-regulated miRNAs and 156 down-regulated miRNAs (Figure 1A-B). The immune score for each sample was determined with the aid of the R package “estimate”. Spearman correlation analysis revealed correlations between each DE-miRNA and immune score. Based on p-values < 0.05 and |cor| > 0.3, a total of 32 immune-related miRNAs in HCC were confirmed (Supplementary Table 1). Functional enrichment analysis indicated that 17 GO entries (BP entries: 3, CC entries: 6, MF entries: 8) and 87 pathways were enriched (Supplementary Table 2). The enriched GO entries and top 10 pathways are shown displayed in Figure 1C-F with involved genes primarily related to “Signal transduction”, “Cell communication”, immune-related, and cancer-related pathways.

Figure 1 Identification of immune-related miRNAs in HCC. (A) The volcano map of DE-miRNAs in HCC. (B) The heatmap of DE-miRNAs in HCC. (C-F) Functional enrichment results of immune-related miRNAs.

Risk Score Models Based on Prognostic miRNAs Was Reliable

To identify immune-related miRNAs influencing the overall survival (OS) in HCC, the 32 immune-related miRNAs were analyzed utilizing univariate Cox analysis in the training set. Ten miRNAs significantly linked to OS (p<0.2) were chosen25 (Figure 2A). These miRNAs underwent stepwise multivariate analysis, resulting in five miRNAs, hsa-miR-145-3p, hsa-miR-150-3p, hsa-miR-153-3p, hsa-miR-223-3p, as well as hsa-miR-424-3p, being considered as prognostic indicators to predict the risk score in the model (Figure 2B). Among them, hsa-miR-153-3p, hsa-miR-223-3p, and hsa-miR-424-3p served as risk factors (Hazard ratio (HR)>1), while hsa-miR-145-3p and hsa-miR-150-3p were protective (HR < 1). This result was further confirmed by Kaplan-Meier (K-M) curves (Figure 2C-G). A miRNA risk signature was constructed on the basis of the risk score formula: RiskScore = (−0.296719854) ×expression (hsa-miR-145-3p) + (−0.187576368) ×expression (hsa-miR-150-3p) + 0.179018998×expression (hsa-miR-153-3p) + 0.179809095×expression (hsa-miR-223-3p) + 0.33640426×expression (hsa-miR-424-3p). The risk scores were computed for every HCC sample in the training set. These samples were then categorized into two risk subgroups (low- and high -risk) relying on median value (Figure 3A). As presented by the K-M curve, higher risk indicated significantly poorer OS compared to lower risk (Figure 3B). ROC curves were created to assess the predictive capacity of this model. In the training set, the AUC values for OS were 0.706 (1 year), 0.653 (3 years), and 0.644 (5 years), suggesting an ideal accuracy (Figure 3C). As demonstrated by the expression heatmap, the expressions of hsa-miR-153-3p, hsa-miR-223-3p, and hs-miR-424-3p were increased in the high risk score group, while the expressions of hsa-miR-145-3p and hsa-miR-150-3p were increased in group with low risk score (Figure 3A). To validate the reliability of the risk score model, the analyses mentioned above were performed by using the external validation set (GSE31384). As demonstrated by the K-M curve, a higher risk score was significantly associated with poorer overall survival (OS) compared to a lower risk score (Figure 3D-E). The predictive capacity of the model was assessed using ROC curves, which showed AUC values of 0.612 (1 year), 0.670 (3 years), and 0.636 (5 years) in the validation set, indicating good accuracy. The external validation results showed a comparable trend set (Figure 3F).

Figure 2 Screening for miRNAs related to the HCC prognosis. (A)Ten miRNAs were selected by the univariate analysis based on OS of HCC patients. (B) Five miRNAs were determined by stepwise multivariate analysis as prognostic miRNAs for risk score model establishment. (C-G) Five miRNAs related K-M survival analyses in HCC patient. * p < 0.05, ** p < 0.01.

Figure 3 Establishment of the immune-related miRNA risk score model for HCC. (A) Risk score curve for each HCC patient was computed based on miRNA risk signature and heatmap for the expression of prognostic miRNAs between two risk subgroups in the training set. (B) K-M curve between two risk subgroups in the training set. (C) Time-dependent ROC curve analysis in the training set. Distribution of risk score (D), heatmap of expression profiles of five miRNAs (D), K-M survival analyses (E) and time-dependent ROC curve analysis (F) was carried in validation set (GSE31384).

Immune Cell Infiltration, Immunotherapy Response, and Chemotherapy Sensitivity Predicted by a Prognostic miRNA Signature in HCC

Using the ssGSEA and Cibersort algorithm, the samples among the two subgroups were analyzed to obtain the score for each infiltrating immune cell type. The results of the ssGSEA algorithm analysis showed that the score of immature DCs, Type 17 T helper (Th) cells, together with Type 2 Th cells were higher in high-risk score group, while the score of Activated B/ CD8 T cells, CD56dim NK cells, immature B cells, monocytes, and Type 1 Th cells were superior in the patients with lower risk scores (Figure 4A). The results of the Cibersort algorithm analysis showed that the enrichment scores of 11 immune infiltrating cells such as B cells memory, Macrophages M0 and others were significantly different based on the high and low risk groups, and the enrichment scores were significantly higher in the high risk group than in the low risk group (Supplementary Figure 2A-B). The results of the correlation between prognosis-related genes and differential immune infiltrating cells showed that the differential immune cells that had a strong correlation (|r| > 0.3) with KCTD17 were Macrophages M0 (r = 0.39, p < 0.001) and Mast cells resting (r = −0.32, p < 0.001). Differential immune cells with a strong correlation (|r| > 0.3) with SFPQ were Macrophages M0 (r = 0.32, p < 0.001). Among the remaining prognostic genes, no differential immune cells with a strong correlation were found (Supplementary Figure 2C). The results of TIDE indicated that low risk score group had inferior TIDE score and T cell exclusion score, implying that patients with low risk may respond more actively to immunotherapy than patients with higher risk scores (Figure 4B-C). Dysfunction scores were higher in the low-risk scoring group than in the high-risk scoring group, whereas CD274 scores were not significantly different between the high- and low-risk groups (Figure 4D-E). Next, we investigated whether the immune-related miRNA signature could predict sensitivity to common chemotherapy drugs for HCC treatment. We noticed that patients with high risk were more sensitive to bleomycin, doxorubicin, and gemcitabine, and those who with low risk was more sensitive to gefitinib, imatinib, and sorafenib (Figure 4F-K).

Figure 4 The relationship between the immune-related miRNA signature and immune cell infiltration and therapy. (A) The box plots displayed the score of immune infiltration cells. (B-E) TIDE score, T cell exclusion score, dysfunction score and CD274 score were carried to identify the immunotherapy response. (F-K) The sensitivity of bleomycin, doxorubicin, gemcitabine (H), gefitinib, imatinib, and sorafenib between two risk subgroups were analysis. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Prognostic miRNA-Regulated mRNAs in HCC

To identify prognostic miRNA-regulated mRNAs, we screened 6265 DE-mRNAs between HCC and non-HCC samples from the TCGA-HCC cohort, including 5411 mRNAs with higher expression as well as 854 mRNAs with lower expression in HCC (Figure 5A). We predicted the mRNAs regulated by the five prognostic miRNAs using the miRwalk database, and took intersection with DE-mRNAs. As shown in Supplementary Figure 3, all five prognostic miRNAs were down-regulated in HCC group. Considering the opposite expression trends of miRNAs and mRNAs, we finally obtained 57 up-regulated prognostic miRNA-regulated mRNAs. The regulatory network was exhibited in Figure 5B. Functional enrichment analysis of these 57 mRNAs using DAVID revealed significant involvement in the AMPK signaling pathway, apoptosis, spindle organization, histone H3 deacetylation, peptidyl-tyrosine phosphorylation, chromatin organization, protein autophosphorylation, protein phosphorylation, and negative regulation of transcription, DNA-templated (Figure 5C).

Figure 5 The prognostic miRNA-regulated mRNAs were identified in HCC. (A) DE-mRNAs were identified by limma package in TGGA-HCC dataset. (B) The regulatory network showed the association between prognostic miRNAs and downstream mRNAs. (C) Functional enrichment analysis of mRNAs in the regulatory network. (D) 32 miRNAs were selected by the univariate Cox analysis based on OS of HCC patients. (E) Six mRNAs were determined as prognostic markers based on LASSO analysis for risk score model establishment.

Risk Score Models Based on Prognostic mRNAs Was Reliable

To identify the mRNAs related with the OS of HCC, we enrolled the 57 prognostic miRNA-regulated mRNAs into univariate Cox analysis of the training set. A total of 32 mRNAs were significantly related with OS of HCC patients (p value<0.05) (Figure 5D). These 32 mRNAs were further integrated into LASSO analysis. As revealed in Figure 5E, six mRNAs (KCTD17, MAFG, RAB10, SFPQ, TRMT6, and UBE2D2) were determined as prognostic mRNAs by 20-fold cross-validation. Next, we generated a mRNA risk signature according to the risk score formula: RiskScore = 0.026473031×expression (KCTD17) + (0.08181457) ×expression (MAFG) + 0.121674715×expression (RAB10) + 0.25088589×expression (SFPQ) + 0.082841154×expression (TRMT6) + 0.05988102×expression (UBE2D2). Afterwards, in the training set, the risk score of each HCC sample was calculated and divided into two risk subgroups (low- and high-risk) according to median value (Figure 6A). The K-M curve demonstrated that high risk group showed worse OS in comparison with low risk group (Figure 6B). In the training set, the AUC values of the ROC curves for OS were 0.782 (1 year), 0.673 (3 years), and 0.629 (5 years), suggesting an ideal accuracy (Figure 6C). The expression heatmap manifested that all six mRNAs were upregulated in the group with high risk (Figure 6A). To verify the reliability of the mRNA risk score model, the aforementioned analysis was repeated using the external validation set (GSE76427). The results from this external validation set were consistent with the trends observed in the initial analysis (Figure 6D-F).

Figure 6 Establishment of the immune-related mRNA risk score model for HCC. (A) Risk score curve for each HCC patient was computed based on mRNA risk signature and heatmap for the expression of prognostic mRNAs between two risk subgroups in the training set. (B) K-M curve between two risk subgroups in the training set. (C) Time-dependent ROC curve analysis in the training set. Distribution of risk score (D), heatmap of expression profiles of mRNAs (D), K-M survival analysis (E) and time-dependent ROC curve analysis (F) was carried out in GSE31384 dataset.

Immune Cell Infiltration, Immunotherapy Response, and Chemotherapy Sensitivity Predicted by a Prognostic mRNA Signature with Differential Expression in HCC Risk Groups

As revealed in Figure 7A, the fractions of activated DCs, CD4 T cells, central and effector memory CD4 T cells, NK T cells, plasmacytoid DCs, as well as Type 2 Th cells were significantly higher in the high risk score group, while the fractions of activated CD8 T cells, CD56dim NK cells, effector memory CD8 T cells, eosinophils, neutrophils, as well as Type 1 Th cell were notably higher in the low risk score group. The TIDE results revealed that low risk score group had inferior TIDE score and T cell exclusion score, suggesting that low risk was related with more positive response to immunotherapy compare with high risk (Figure 7B-C). In contrast, Dysfunction and CD274 scores were higher in the low-risk scoring group than in the high-risk scoring group. (Figure 7D-E). We then investigated whether the mRNA signature could predict sensitivity to common chemotherapy drugs for HCC treatment. We found that the group of high risk was more sensitive to bleomycin, cisplatin, doxorubicin, etoposide, gemcitabine, mitomycin C, rapamycin, and vinblastine, and the group of low risk was more sensitive to gefitinib (Figure 7F-N).

Figure 7 The relationship between the immune-related mRNA signature and immune cell infiltration and therapy. (A) The box plot displaying the score of immune infiltration cells. (B-E) TIDE score, T cell exclusion score, dysfunction score and CD274 score were deployed to identify the immunotherapy response between two risk subgroups. (F-N) The sensitivity between two risk subgroups to bleomycin, cisplatin, doxorubicin, etoposide, gemcitabine, mitomycin C, rapamycin, vinblastine and gefitinib between two risk subgroups were analysis. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Analysis on the Relationship Between Clinical Indicators and miRNA Risk Score

To discover the association between clinical indicators and the immune-associated miRNA signature, we proceeded to compare the risk scores of the subgroups of clinical indicators. In Supplementary Figure 4, the miRNA risk score was related with gender, pathological T stage, stage, and vascular tumor cell type (p value < 0.05), but not associated with age, histologic grade of neoplasm, N stage, and M stage.

The Immune-Related miRNA and mRNA Risk Score as Independent Prognostic Predictors

To further investigate the influence of clinicopathological characteristics on the prognosis of the risk model, we included clinical and pathological factors (age, gender, patho_M, patho_N, patho_T, tumor histological grade, tumor stage, vascular tumor cell type and RiskScore) of the 368 hCC samples in the TCGA-mRNA training dataset in the Cox regression analysis. Based on the results of univariate Cox analysis, factors with p-values less than 0.05 were selected for multivariate Cox analysis. The results showed that four clinical factors were significantly predictive: tumor stage, Patho_T, Patho_M and RiskScore. And by univariate and multivariate Cox analysis, we confirmed that both immune-associated miRNA and mRNA risk score were independent predictors of HCC prognosis (Supplementary Figure 5, p value<0.05).

The ceRNA Network of Prognostic miRNAs and mRNAs

Initially, 445 DE-lncRNAs between HCC and non-HCC samples were selected, containing 376 lncRNAs with high expression and 69 lncRNAs with low expression in the HCC samples (Supplementary Figure 6). Subsequently, we predicted the lncRNAs targeting the five prognostic miRNAs through the LncBase V2.0 database and intersected them with DE-lncRNAs. Considering the opposite expression trends of miRNAs and lncRNAs, we ultimately identified 63 up-regulated lncRNAs targeting the prognostic miRNAs. Thus, the regulatory network of lncRNA-miRNA-mRNA comprising 63 lncRNAs, 5 miRNAs, 57 mRNAs, and 822 lncRNA-miRNA-mRNA interaction pairs, was constructed and presented in Supplementary Figure 7 and Supplementary Table 3.

Validation of the Level of Prognosis-Related mRNAs in Cell Lines

The protein expression levels of the six prognosis-related genes were shown in Supplementary Figure 8, with KCTD17 expressing the highest in the Choroid plexus, MAFG, SFPQ, and TRMT6 expressing the highest in Bone marrow, and RAB10 and UBE2D2 expressing the highest in Skeletal muscle. The mRNA expression levels of prognosis-related genes were checked in human normal hepatocytes (QSG-7701) and three human HCC cell lines (HepG2, LM3, and HuH-7). Similarly, the expressions of KCTD17, MAFG, RAB10, SFPQ, TRMT6, and UBE2D2 were all increased in these HCC cell lines compared to non-HCC cells (Figure 8).

Figure 8 The expression of prognostic mRNAs in cell lines was measured by qPCR. Expression of KCTD17 (A), MAFG (B), RAB10 (C), SFPQ (D), TRMT6 (E), UBE2D2 (F) in QSG-7701 and three human hepatocellular carcinoma cell lines (HepG2, LM3 and HuH-7). * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Discussion

Recent research has shown that increasing antitumor immune response reveal a novel and promising approach for advanced HCC.26 Revolutions in high-throughput DNA sequencing have enabled the discovery of tumor-specific antigens or neoantigens that are exclusively expressed in cancer cells, absent in healthy cells, and capable of triggering an immune response. Consequently, these findings hold significant potential for cancer vaccine-based immunotherapy strategies.27 Growing evidence proved that ICIs could be the most beneficial and effective treatment method for malignancies.28 In recent years, there has been substantial progress in the development of prognostic models for HCC immunotherapy. The development of a novel prognostic gene signature and nomogram for HCC, facilitated by bioinformatics analysis of the drug target tanshinone IIA, which was validated for prognostic value and could inform immunotherapeutic decisions for HCC patients.29 Meanwhile, advancements in radiological imaging have demonstrated significant predictive value in the context of HCC immunotherapy.30,31 In particular, the use of multi-sequence magnetic resonance imaging (MRI) has demonstrated potential as a non-invasive imaging biomarker for predicting PD-1/PD-L1 expression levels, which could help assess responsiveness to immunotherapeutic treatments.32 However, the lack of effective immunological prognostic indicators remains a significant challenge.

Furthermore, evidence has shown that miRNAs play a necessary role in the immune system development, response, and programs activation, particularly in regulatory functions.33,34 Thus, identification of immune-related miRNAs is helpful in finding progressive biomarkers that could confirm the sensitivity of HCC patients to ICI therapy. In our research, 32 immune-associated miRNAs in HCC were identified. Subsequently, through Cox regression analysis, five prognostic miRNAs (hsa-miR-145-3p, hsa-miR-150-3p, hsa-miR-153-3p, hsa-miR-223-3p, and hsa-miR-424-3p) were identified to establish an immune-related miRNA signature. Among them, miR-145 is a common cancer suppressor. It has been demonstrated that miR-145 directly targets the 3′-UTR of PD‐L1 to suppress PD‐L1 expression, ultimately inhibiting invasion, epithelial-mesenchymal transition (EMT), and migration of HCC cells.35 In exosomes derived from cancer-associated fibroblasts (CAFs), miR-150-3p remarkably promotes the progression of HCC.36 Moreover, miR-153-3p is decreased in HCC cells and tissues and plays an ontogenetic role by mediating Rabl3.37 Previous studies have shown that miR-223-3p exerts an anti-HCC effect via targeting CXCR4.38 Interestingly, miR-223 has been reported to prevent nonalcoholic steatohepatitis (NASH) from progressing into liver malignancy through regulating certain inflammatory and oncogenic genes, including CXCL-10 and Taz.39 The expression level of serum miR-424 is decreased in patients with HCC and negatively related to tumor progression.40 The lncRNA MYLK-AS1 promotes angiogenesis of HCC tumors via directly targeting miR-424-5p.41 However, the function of these five prognostic miRNAs in mediating the immune microenvironment of HCC still needs further investigation.

To further understand the related regulatory pathways, 57 prognostic miRNA-regulated mRNAs in HCC were determined. Then, a total of six mRNAs (KCTD17, MAFG, RAB10, SFPQ, TRMT6, and UBE2D2) were authenticated to establish an mRNA signature via LASSO regression as well as univariate Cox analysis. Meanwhile, our research showed that the expressions of all of these mRNAs were increased in these HCC cell lines in contrast with to non-cancer cells. KCTD17, as an E3 ligase component protein for trichoplein, promotes hepatic steatosis via binding to phosphorylated PHLPP2 to regulate lipogenesis.42 Furthermore, mRNA levels of KCTD17 were found to decrease in relapsed HCC, which is related to cell proliferation, differentiation, and survival.43 The expression of MAFG has been found to be increased in HCC,44 but its specific features remain unknown. Previous research reported that Rab10 played a part in the progression of various cancers.45 In HCC, downregulation of Rab10 decreases the proliferation of cells and induces their autophagy and apoptosis.46 Rab10 could cooperate with LRRK2 to regulate the action of macrophages, mediating the immune response of phagocytes.47 Splicing factor proline and glutamine rich (SFPQ) was found to modulate platinum response via regulating SRSF2 activity in ovarian cancer.48 However, the functions of SFPQ, TRMT6, and UBE2D2 in HCC deserve further in-depth studies.

In HCC, different immune cells play crucial roles in tumor progression, including T cells, macrophages and NK cells, formation of the TME, and modulation of immune responses.49 Particularly, CD8+ T cell is important in the tumor microenvironment.50 Tumors with the infiltration of various kinds of tumor-specific CD8+ T cells can exhibit clinical responses to ICIs in multiple cancer types.51 In non–small cell lung cancer, tumor-infiltrating CD8+ T cells regulate the response to ICIs via mediating the expression level of programmed cell death ligand 1.52 NK cells, such as CD8+ T cells, provide anti-tumor immunity by releasing perforin and granulein, causing apoptosis.53 Besides, they also stimulate anti-cancer responses with cytokines and chemokines. Studies suggest that the depletion of NK cells to be a major reason for the failure of PD-1-based ICIs among treating patients with HCC.54 The results of TIDE demonstrated that group with low risk showed more positive response to immunotherapy, while the fraction of activated CD8+ T cells as well as CD56dim NK cells were notably elevated in the patients with lower risk scores. TIDE and drug sensitivity analysis further suggested that the signatures might be able to act as immunotherapy and chemosensitivity predictors. The results indicated that patients identified as being at high risk demonstrated heightened sensitivity to bleomycin, doxorubicin, and gemcitabine. This enhanced responsiveness may correlate with the expression of specific biomarkers within their tumors, suggesting that these agents could serve as viable clinical options for patients categorized as high risk. In contrast, patients classified as low risk exhibited greater sensitivity to gefitinib, imatinib, and sorafenib, possibly due to their unique molecular characteristics. These findings not only highlight the significance of personalized therapeutic strategies informed by patient risk stratification but also point to the necessity for future research to delineate the most efficacious treatment protocols for patients within varying risk strata. Based on these results, it could be inferred that the mRNAs and miRNAs might be promising indicators to predict the efficacy of immunotherapy.

To reveal potential regulatory pathways, we predicted the ceRNA network of the prognostic miRNAs and mRNAs. Currently, none of these ceRNA molecular pathways have been reported in HCC immune regulation. This ceRNA network will provide the direction and basis for investigating the molecular mechanisms by which these prognostic miRNAs and mRNAs influence HCC progression.

While this study validated immune-related prognostic signatures and HCC risk profiles, it has limitations. The sample size and diversity were insufficient to represent the heterogeneity of all HCC patients, and variations in geographical, genetic, and clinical backgrounds may limit the universality of the findings. The choice of external datasets, sample size, and quality control could affect the reliability of the validation, necessitating further validation with diverse datasets. Additionally, the ssGSEA analysis may not fully elucidate the relationship between risk scores and immune cell infiltration. Future studies should employ flow cytometry for quantitative immune cell analysis and spatial transcriptomics to explore the microenvironmental distribution and functional interactions of immune cells. Clinical trials are required to further validate predictions of drug sensitivity and immunotherapy response. Future work will aim to integrate these insights and refine the models for improved clinical utility.

Conclusion

In summary, this study discovered and validated immune-related prognostic microRNAs (miRNAs) and messenger RNAs (mRNAs) in Hepatocellular Carcinoma (HCC). Key findings include five miRNAs (hsa-miR-145-3p, hsa-miR-150-3p, hsa-miR-153-3p, hsa-miR-223-3p, hsa-miR-424-3p) and six mRNAs (KCTD17, MAFG, RAB10, SFPQ, TRMT6, UBE2D2) that predict HCC prognosis. The innovation combines bioinformatics, experimental validation, and functional assessments, enhancing our understanding of HCC immune mechanisms and paving the way for personalized prognostic markers and therapies. These results promise better management and treatment for HCC patients.

Data Sharing Statement

The sequencing information and clinical data on the TCGA-HCC cohort were acquired from TCGA database (https://xenabrowser.net) on the April of 2022. Two HCC datasets, GSE31384 and GSE76427, were acquired from GEO database (https://www.ncbi.nlm.nih.gov/).

Ethics Approval and Consent to Participate

The present study was authorized by the Research Ethics Review Board of The Nantong City No 1 People’s Hospital and Second Affiliated Hospital of Nantong University, (No. 2023KT022).

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study was funded by National Natural Science Foundation of China (81972815) and the Nantong Basic Scientific Research and Social Livelihood Science and Technology Plan Projects (MSZ2022041).

Disclosure

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

1. Forner A, Reig M, Bruix J. Hepatocellular carcinoma. Lancet. 2018;391(10127):1301–1314. doi:10.1016/S0140-6736(18)30010-2

2. Hao X, Sun G, Zhang Y, et al. Targeting immune cells in the tumor microenvironment of HCC: new opportunities and challenges. Front Cell Dev Biol. 2021;9:775462. doi:10.3389/fcell.2021.775462

3. Keenan BP, Fong L, Kelley RK. Immunotherapy in hepatocellular carcinoma: the complex interface between inflammation, fibrosis, and the immune response. J Immunother Cancer. 2019;7(1):267. doi:10.1186/s40425-019-0749-z

4. Krenkel O, Tacke F. Liver macrophages in tissue homeostasis and disease. Nat Rev Immunol. 2017;17(5):306–321. doi:10.1038/nri.2017.11

5. Chen Y, Hu H, Yuan X, Fan X, Zhang C. Advances in immune checkpoint inhibitors for advanced hepatocellular carcinoma. Front Immunol. 2022;13:896752. doi:10.3389/fimmu.2022.896752

6. Lee CK, Chan SL, Chon HJ. Could we predict the response of immune checkpoint inhibitor treatment in Hepatocellular Carcinoma? Cancers. 2022;14(13). doi:10.3390/cancers14133213

7. Yao C, Wu S, Kong J, et al. Angiogenesis in Hepatocellular Carcinoma: mechanisms and anti-angiogenic therapies. Cancer Biol Med. 2023;20(1):25–43. doi:10.20892/j.issn.2095-3941.2022.0449

8. Rupaimoole R, Slack FJ. MicroRNA therapeutics: towards a new era for the management of cancer and other diseases. Nat Rev Drug Discov. 2017;16(3):203–222. doi:10.1038/nrd.2016.246

9. Dietrich P, Koch A, Fritz V, Hartmann A, Bosserhoff AK, Hellerbrand C. Wild type Kirsten rat sarcoma is a novel microRNA-622-regulated therapeutic target for hepatocellular carcinoma and contributes to sorafenib resistance. Gut. 2018;67(7):1328–1341. doi:10.1136/gutjnl-2017-315402

10. Fritz V, Malek L, Gaza A, et al. Combined de-repression of chemoresistance associated mitogen-activated protein kinase 14 and activating transcription factor 2 by loss of microRNA-622 in Hepatocellular Carcinoma. Cancers. 2021;13(5):1183. doi:10.3390/cancers13051183

11. Khalaf K, Hana D, Chou JT, Singh C, Mackiewicz A, Kaczmarek M. Aspects of the tumor microenvironment involved in immune resistance and drug resistance. Front Immunol. 2021;12:656364. doi:10.3389/fimmu.2021.656364

12. Lou Q, Liu R, Yang X, et al. miR-448 targets IDO1 and regulates CD8(+) T cell response in human colon cancer. J Immunother Cancer. 2019;7(1):210. doi:10.1186/s40425-019-0691-0

13. Zhou C, Wei W, Ma J, et al. Cancer-secreted exosomal miR-1468-5p promotes tumor immune escape via the immunosuppressive reprogramming of lymphatic vessels. Mol Ther. 2021;29(4):1512–1528. doi:10.1016/j.ymthe.2020.12.034

14. Li Y, Ma H, Shi C, Feng F, Yang L. Mutant ACTB mRNA 3’-UTR promotes hepatocellular carcinoma development by regulating miR-1 and miR-29a. Cell Signal. 2020;67:109479. doi:10.1016/j.cellsig.2019.109479

15. Grinchuk OV, Yenamandra SP, Iyer R, et al. Tumor-adjacent tissue co-expression profile analysis reveals pro-oncogenic ribosomal gene signature for prognosis of resectable hepatocellular carcinoma. Mol Oncol. 2018;12(1):89–113. doi:10.1002/1878-0261.12153

16. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106. doi:10.1186/gb-2010-11-10-r106

17. Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9(10):1657–1659. doi:10.1046/j.1365-294x.2000.01020.x

18. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 2013;14:7. doi:10.1186/1471-2105-14-7

19. Ginestet C. ggplot2: elegant graphics for data analysis. J Royal Stat Soc Series A. 2011;174(1):245–246. doi:10.1111/j.1467-985X.2010.00676_9.x

20. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–1558. doi:10.1038/s41591-018-0136-1

21. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021;22(6). doi:10.1093/bib/bbab260

22. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–2504. doi:10.1101/gr.1239303

23. Uhlen M, Fagerberg L, Hallstrom BM, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419. doi:10.1126/science.1260419

24. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25(4):402–408. doi:10.1006/meth.2001.1262

25. Ye Y, Dai Q, Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov. 2021;7(1):71. doi:10.1038/s41420-021-00451-x

26. Fasano R, Shadbad MA, Brunetti O, et al. Immunotherapy for Hepatocellular Carcinoma: new prospects for the cancer therapy. Life. 2021;11(12). doi:10.3390/life11121355

27. Morazan-Fernandez D, Mora J, Molina-Mora JA. In silico pipeline to identify tumor-specific antigens for cancer immunotherapy using exome sequencing data. Phenomics. 2023;3(2):130–137. doi:10.1007/s43657-022-00084-9

28. Giraud J, Chalopin D, Blanc JF, Saleh M. Hepatocellular Carcinoma immune landscape and the potential of immunotherapies. Front Immunol. 2021;12:655697. doi:10.3389/fimmu.2021.655697

29. Peng B, Ge YUN, Yin G. A novel prognostic gene signature, nomogram and immune landscape based on tanshinone IIA drug targets for hepatocellular carcinoma: comprehensive bioinformatics analysis and in vitro experiments. Biocell. 2023;47(7):1519–1535. doi:10.32604/biocell.2023.027026

30. Tao YY, Shi Y, Gong XQ, et al. Radiomic analysis based on magnetic resonance imaging for predicting PD-L2 expression in Hepatocellular Carcinoma. Cancers. 2023;15(2):365. doi:10.3390/cancers15020365

31. Li YK, Wu S, Wu YS, et al. Portal venous and hepatic arterial coefficients predict post-hepatectomy overall and recurrence-free survival in patients with hepatocellular carcinoma: a retrospective study. J Hepatocell Carcinoma. 2024;11:1389–1402. doi:10.2147/JHC.S462168

32. Gong XQ, Liu N, Tao YY, et al. Radiomics models based on multisequence MRI for predicting PD-1/PD-L1 expression in hepatocellular carcinoma. Sci Rep. 2023;13(1):7710. doi:10.1038/s41598-023-34763-y

33. Essandoh K, Li Y, Huo J, Fan GC. MiRNA-mediated macrophage polarization and its potential role in the regulation of inflammatory response. Shock. 2016;46(2):122–131. doi:10.1097/SHK.0000000000000604

34. Rui X, Wu Q, Gong Y, Wu Y, Chi Q, Sun DA. A novel prognostic target-gene signature and nomogram based on an integrated bioinformatics analysis in hepatocellular carcinoma. Biocell. 2022;46(5):1261–1288. doi:10.32604/biocell.2022.018427

35. Xu G, Zhang P, Liang H, et al. Circular RNA hsa_circ_0003288 induces EMT and invasion by regulating hsa_circ_0003288/miR-145/PD-L1 axis in hepatocellular carcinoma. Cancer Cell Int. 2021;21(1):212. doi:10.1186/s12935-021-01902-2

36. Yugawa K, Yoshizumi T, Mano Y, et al. Cancer-associated fibroblasts promote hepatocellular carcinoma progression through downregulation of exosomal miR-150-3p. Eur J Surg Oncol. 2021;47(2):384–393. doi:10.1016/j.ejso.2020.08.002

37. Qi WY, Mao XB, He YB, Xiao CH. Long non-coding RNA LINC00858 promotes cells proliferation and invasion through the miR-153-3p/Rabl3 axis in hepatocellular carcinoma. Eur Rev Med Pharmacol Sci. 2020;24(18):9343–9352. doi:10.26355/eurrev_202009_23017

38. Si H, Wang H, Xiao H, Fang Y, Wu Z. Anti-tumor effect of celastrol on Hepatocellular Carcinoma by the circ_SLIT3/miR-223-3p/CXCR4 axis. Cancer Manag Res. 2021;13:1099–1111. doi:10.2147/CMAR.S278023

39. He Y, Hwang S, Cai Y, et al. MicroRNA-223 ameliorates nonalcoholic steatohepatitis and cancer by targeting multiple inflammatory and oncogenic genes in hepatocytes. Hepatology. 2019;70(4):1150–1167. doi:10.1002/hep.30645

40. Yang C, Du P, Lu W. MiR-424 acts as a novel biomarker in the diagnosis of patients with Hepatocellular Carcinoma. Cancer Biother Radiopharm. 2021;38:670–673. doi:10.1089/cbr.2020.4141

41. Teng F, Zhang JX, Chang QM, et al. LncRNA MYLK-AS1 facilitates tumor progression and angiogenesis by targeting miR-424-5p/E2F7 axis and activating VEGFR-2 signaling pathway in hepatocellular carcinoma. J Exp Clin Cancer Res. 2020;39(1):235. doi:10.1186/s13046-020-01739-z

42. Kim K, Ryu D, Dongiovanni P, et al. Degradation of PHLPP2 by KCTD17, via a glucagon-dependent pathway, promotes hepatic steatosis. Gastroenterology. 2017;153(6):1568–1580e10. doi:10.1053/j.gastro.2017.08.039

43. Fang Y, Yang Y, Zhang X, et al. A co-expression network reveals the potential regulatory mechanism of lncRNAs in relapsed Hepatocellular Carcinoma. Front Oncol. 2021;11:745166. doi:10.3389/fonc.2021.745166

44. Liu T, Yang H, Fan W, et al. Mechanisms of MAFG dysregulation in cholestatic liver injury and development of liver cancer. Gastroenterology. 2018;155(2):557–571e14. doi:10.1053/j.gastro.2018.04.032

45. Li Y, Xiong Y, Wang Z, et al. FAM49B promotes breast cancer proliferation, metastasis, and chemoresistance by stabilizing ELAVL1 protein and regulating downstream Rab10/TLR4 pathway. Cancer Cell Int. 2021;21(1):534. doi:10.1186/s12935-021-02244-9

46. Zhang YJ, Pan Q, Yu Y, Zhong XP. microRNA-519d induces autophagy and apoptosis of human Hepatocellular Carcinoma cells through activation of the AMPK signaling pathway via Rab10. Cancer Manag Res. 2020;12:2589–2602. doi:10.2147/CMAR.S207548

47. Liu Z, Xu E, Zhao HT, Cole T, West AB. LRRK2 and Rab10 coordinate macropinocytosis to mediate immunological responses in phagocytes. EMBO J. 2020;39(20):e104862. doi:10.15252/embj.2020104862

48. Pellarin I, Dall’Acqua A, Gambelli A, et al. Splicing factor proline- and glutamine-rich (SFPQ) protein regulates platinum response in ovarian cancer-modulating SRSF2 activity. Oncogene. 2020;39(22):4390–4403. doi:10.1038/s41388-020-1292-6

49. Feng M, Wang F, Liu X, et al. Neutrophils as key regulators of tumor immunity that restrict immune checkpoint blockade in liver cancer. Cancer Biol Med. 2023;20(6):421–437. doi:10.20892/j.issn.2095-3941.2023.0019

50. Hossain MA, Liu G, Dai B, et al. Reinvigorating exhausted CD8(+) cytotoxic T lymphocytes in the tumor microenvironment and current strategies in cancer immunotherapy. Med Res Rev. 2021;41(1):156–201. doi:10.1002/med.21727

51. Jenkins L, Jungwirth U, Avgustinova A, et al. Cancer-associated fibroblasts suppress CD8+ T