The identification of N6-methyladenosine-related miRNAs predictive of hepatocellular carcinoma prognosis and immunotherapy efficacy

Abstract

BACKGROUND:

Hepatocellular carcinoma (HCC) has a high degree of malignancy and poor prognosis. N6-methyladenosine (m6A) modifications and microRNAs (miRNAs) play pivotal roles in tumorigenesis and development. However, the role of m6A-related miRNAs in HCC has not been clarified yet. This study aimed to identify the role of m6A-miRNAs in HCC prognosis through bioinformatics analysis.

METHODS:

The clinicopathological information and RNA sequencing data of 369 HCC tumor tissues and 49 tumor-adjacent tissues were downloaded from the TCGA database. A total of 23 m6A regulators were extracted to evaluated the m6A-related miRNAs using Pearson’s correlation analysis. Then, we selected prognosis-related m6A-miRNAs using a univariate Cox regression model and used the consensus cluster analysis to explore the characteristics of the m6A-miRNAs. The coefficient of the least absolute shrinkage and selection operator (LASSO) Cox regression was applied to construct a prognostic risk score model. The receiver operated characteristic (ROC) analysis was applied to evaluate the prognostic value of the signature. The biological functions of targeted genes were predicted by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. Then, to validate the potential predictive value for prognosis, the miRNA expression profiles from the GSE76903 and GSE6857 were used. Single sample Gene Set Enrichment Analysis (ssGSEA) and Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) were applied to assess the immune microenvironment of HCC. Additionally, a meta-analysis was used to verify the prognostic value of the m6A-microRNAs. RT-PCR was applied to validated the expression of miRNAs in HCC tissues. Cell viability, transwell assay and RNA m6A dot blot assays of HCC cells was applied to access the function of miR-17-5p.

RESULTS:

The expression of 48 m6A-related miRNAs was identified and 17 prognostic m6A-miRNAs was discovered. The expression profile of those 17 miRNAs was divided into three clusters, and these clusters were associated with the tumor microenvironment (TME) and prognosis. The nine m6A-related miRNA signature was associated with the prognosis of HCC, the AUC of the ROC was 0.771(TCGA dataset), 0.788(GSE76903) and 0.646(GSE6857). The TME and the expression of immune checkpoint molecules were associated with the risk score. The meta-analysis also validated the prognostic value of the m6A-related miRNAs (miR182-5p (HR:1.58, 95%CI:1.04-2.40) and miR-17-5p (HR:1.58, 95%CI: 1.04–2.40)). The expression of miR-17-5p was upregulated in HCC tissues and miR-17-5p showed an oncogenic role in HCC cells.

CONCLUSION:

The clinical innovation is the use of m6A-miRNAs as biomarkers for predicting prognosis regarding immunotherapy response in HCC patients.

1. Background

Liver cancer is a common and fatal malignant tumor. Globally, in 2020, 905,611 new cases of liver cancer were diagnosed and 830,180 new deaths due to liver cancer were reported [1]. Hepatocellular carcinoma (HCC) was the most commonly reported form of liver cancer [2], and there is a lack of effective treatments for HCC, especially in the advanced stages. At present, the prognosis of HCC is poor, and there is lack of effective targets for therapy.

N6-methyladenosine (m6A) is a type of the methylation modification in which the sixth nitrogen (N) atom of adenine (A) is modified. As a primary RNA modification, especially in eukaryotic mRNA, it affects the biological function of RNA [3]. There is increasing evidence indicates that the abnormal modification of m6A is involved in the progression of many cancers, including HCC [4].

MicroRNAs (miRNAs) are a class of non-coding single stranded RNA molecules of a length of 18–25 nucleotides. They are involved in many biological processes, including the development of cancer, by regulating gene expression [5]. Yan et al. demonstrated that m6A modifications lead to miRNA dysregulation and affect tumor metastasis and progression [6]. Moreover, m6A-targeted miRNAs have been demonstrated to serve as effective biomarkers for cancer detection [7]. However, the role of m6A-miRNAs in the prognosis of HCC has not been studied yet.

The tumor microenvironment (TME) consists of tumor cells, immune cells and stromal cells [8]. Accumulating studies have highlighted the important role of the TME in cancer prognoses, including the prognosis of HCC [9]. Arvanitakis et al. indicated that the TME affects the response to immunotherapies for HCC [10]. Immune checkpoint inhibitors (ICIs) have recently been a source of promising new cancer treatments when combined with anti-angiogenic agents, and and their use in the treatment of advanced HCC is effective [11]. The most extensively studied checkpoints are the programmed death protein 1(PD-1)/programmed death receptor ligand 1(PD-L1) [12].

In this study, we try to investigate the m6A related miRNAs and to evaluate the association with HCC prognosis. The signature of m6A-miRNA for HCC prognosis was explored based on The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) dataset, and then, a meta-analysis was applied to validate the prognostic value. Furthermore, we investigated the relationship between TME and m6A-miRNAs to evaluate its potential use in predicting immunotherapy efficacy, the work flow showed in Supplementary Figure 1

2. Materials and methods

2.1 Data collection and expression profile preparation

A total of 369 HCC patients with survival information were included in our study(Supplementary Table 1). RNA sequencing data of 369 HCC tumor tissues and 49 tumor-adjacent tissues were downloaded from the TCGA database (https://portal.gdc.cancer.gov). After systematic review [13, 14, 15], 23 m6A regulators (writers: METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B; readers: YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, and RBMX; erasers: FTO and ALKBH5) were extracted. The differential expression of these regulators between tumor and tumor-adjacent samples were assessed using the “limma” package, and those with a $P$ value $<$ 0.05 were selected for further analysis.

2.2 Risk model construction

To identify the m6A-related miRNAs (m6A-miRNAs), we performed Pearson’s correlation analysis, and the m6A-miRNAs with a Pearson’s coefficient $|$ R $|$ $>$ 0.6 and $P$ value $<$ 0.05 were considered significant. The correlations between m6A and miRNAs were depicted using the “circlize” package [16] . Then, we selected prognosis-related m6A-miRNAs ( $P$ value $<$ 0.05) by univariate Cox regression model. Based on these prognostic miRNAs, we divided patients into different subgroups after 1000 resamplings using the “ConsensusClusterPlus” package [17].

To avoid overfitting, the least absolute shrinkage and selection operator (LASSO) Cox regression was conducted to select prognosis-specific m6A-miRNAs regulators for establishing the risk model using the “glmnet” package. We constructed the risk model as follows: risk score $=$ ${\sum }_{i=1}^n{\beta }_i{\text{𝑒𝑥𝑝}}_i$ , where ${\beta }_i$ represents the regression coefficient, and ${\text{𝑒𝑥𝑝}}_i$ represents the z-transformed expression level of each prognosis-specific m6A-miRNA. Then, samples were divided into low and high-risk groups according to the median value of the risk score, and the overall survival (OS) rates of these two groups were compared by univariate Cox regression.

In addition, univariate and multivariate Cox regression analyses were used to identify independent prognostic factors and construct an integrated nomogram that included predictable clinicopathological factors and risk scores.

2.3 Functional enrichment analysis

To explore the potential biological function, we search the target genes of prognosis-specific m6A-miRNAs from The Encyclopedia of RNA Interactomes (ENCORI) database [18] and the screened differentially expressed genes (DEGs) between the low- and high-risk groups based on $|$ log2FC $|⩾$ 1 and $P$ value $<$ 0.05. The Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) functional enrichment analyses were performed to annotate these genes, and the Benjamini-Hochberg method was applied for false discovery rate (FDR) estimation through “clusterProfiler” R package [19]. Pathways with FDR $<$ 0.05 were considered significant.

2.4 Risk model validation

GSE76903 and GSE6857 downloaded from the Gene Expression Omnibus (GEO) database (https://www. ncbi.nlm.nih.gov/geo), were used as the validation set [20]. The miRNA-seq data of paired adjacent normal tissues and portal vein tumor thrombosis (PVTT) from 20 hepatocellular carcinoma patients were obtained to evaluate the feasibility of the risk model (GSE76903). GSE6857, we downloaded the normalized miRNA expression data from 240 HCC patients, and the platform of the miRNA data was from OSU-CCC MicroRNA Microarray Version 2.0, miRNA-652-3p was not included.

2.5 Immune microenvironment assessment

Single sample Gene Set Enrichment Analysis (ssGSEA) was performed to explore the different infiltration degrees of the immune cell types, immune-related functions, and immune-related pathways between clustering subtypes in TCGA using the “GSVA” R package [21]. Correlations were analyzed to explore the relevance between the risk score and immune infiltration and were depicted by the “corrplot” R package. Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) [22] was applied to evaluate the tumor cellularity and tumor purity through four indicators, including the stromal score, immune score, ESTIMATE score and tumor purity, based on the gene expression profiles of the HCC samples in TCGA.

The programmed death receptor 1 (PD1) and its ligand programmed death receptor ligand 1 (PDL1) are the most widely used immune checkpoints in cancer immunotherapy. To explore their potential correlation with prognosis-specific m6A-miRNAs, we compared the differences in their expression levels among the clustering subtypes and risk subgroups.

2.6 Meta-analysis of the prognostic value of the m6A-miRNAs

We conducted a meta-analysis to validate the results of the prognostic values of the related m6A-miRNAs in HCC in accordance with the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines. Two authors (Ya qian Liu and Hangju Zhu) searched the articles independently. We searched relevant studies published until April 3, 2022 on the electronic databases PubMed and Web of Science. The search keywords were as follows: “hepatocellular carcinoma” or “liver cancer” and “prognosis” or “recurrent” or “survival” and “let-7b-5p” or “miR-148a-3p” or “miR-17-5p” or “miR-182-5p” or “miR-212-3p” or “miR-22-3p” or “miR-652-3p” or “miR-9-5p” or “miR-99b-3p”. The information collected for each study was as follows: author name, publication date, country, test type, hazard ratio (HR) confident intervals (CI) calculated by survival analysis, sample size, and outcome endpoint. All studies included in the current meta-analysis met the following criteria: (a) the microRNA expression profiles were from tissues or circulation; (b) the study performed a survival analysis to examine the association between miRNAs and outcomes with a HR and 95% CI reported; and (c) studies published in English as original articles. The exclusion criteria were as follows: (a) the study tested the prognostic role of miRNA target genes instead of the miRNAs; (b) the study was a meta-analysis, review, comment, letter or duplicate publication; and (c) failure to obtain data extraction. The combined HRs and 95% CI was calculated to evaluate the association between m6A-miRNAs and OS. Moreover, the I-squared test was used to assess the heterogeneity among the studies (I ${}^2$ $<$ 50%, random effect model; I ${}^2$ $>$ 50% fixed-effects model).

2.7 HCC Samples

Seven HCC patients from Jiangsu Cancer Hospital were enrolled in the current study between January 2022 and December 2022. HCC tumor tissues and adjacent tissues without tumor were formalin-fixed paraffin-embedded tissues (FFPET). This study was performed following the procedure from the Ethics Committee of Jiangsu cancer hospital and with informed consent. 32 paired HCC tissues and adjacent normal tissues from GEO dataset (GSE227378) with non-coding RNA sequencing data.

2.8 Cell culture

Human liver cancer cells HepG2, Hep3B, Luh7 were acquired from the Cell Bank of the Chinese Academy of Sciences. Those cells were maintained in DMEM supplemented with 10% foetal bovine serum (Thermo Fisher Scientific), 80 U/mL penicillin and 0.08 mg/mL streptomycin (KeyGEN BioTECH), and incubated in a humidified 5% CO2 atmosphere at 37 ${}^{\circ }$ C.

2.9 Cell transfection

The mimic and inhibitor of miR-17-5p were synthesized from Ribobio (Guangzhou, China). HepG2, Hep3B and Luh7 cells were transfected with 20 nM mimic/inhibitor using Lipofectamine™ RNAiMAX (Thermo Fisher Scientific, Inc.) according to the manufacturer’s protocol at 37 ${}^{\circ }$ C in a 5% CO2 incubator. Following 48 h of transfection, cells were collected for cell proliferation and RT-qPCR analysis. All the transfections were repeated $>$ 3 times independently. Negative control (NC), miR-NC was used as the control.

2.10 Cell viability

Cell viability was measured using the Cell Counting Kit (CCK)-8 (Dojindo Laboratories, Inc.) assay. Briefly, Human liver cancer cells (HepG2, Hep3B and Luh7) were seeded into a 96-well plate at a density of 3,000 cells/well and cultured at 37 ${}^{\circ }$ C for 24 h prior to transfection. Then, cells were transfected with corresponding oligonucleotides. After 24, 48 or 72 h, 10 $\mu$ l CCK-8 solution was added into each well of the plate. The plates were incubated at 37 ${}^{\circ }$ C for 1 h, and the absorbance at 450 nm was measured. The experiment was repeated three times.

2.11 Transwell assay

Transwell assay was used to evaluate the invasion and migration capacities of HCC cells in vitro. Cells at a concentration of 1 $\times$ 10 ${}^5$ cells in 500 $\mu$ l of serum-free seeded in the upper chamber, coated with or without growth factor reduced Matrigel (Costar, USA) in 24-well culture plates, 10% FBS was added into the lower chamber. After incubation for 24 h, cells on the upper surface of the membrane were removed. Then the the cells were fixed with formaldehyde and stained using 0.5% crystal violet. The numbers of cells were counted in five randomly selected fields under a microscope.

2.12 RNA isolation and Reverse transcription-quantitative PCR (RT-PCR)

Total RNA of the cultured cells was extracted using TRIzol reagent (Thermo Fisher Scientific, Inc.) and FFPE tissues was extracted by GIAGEN miRNeasy FFPE Kit (No.217504). The TaqMan Reverse Transcription kit (Takara Biotechnology Co., Ltd.) was used to obtain cDNA for mRNA detection (42 ${}^{\circ }$ C for 1 h), while the TaqMan MicroRNA Reverse Transcription kit (Takara Biotechnology Co., Ltd.) was used for miRNA detection (42 ${}^{\circ }$ C for 1 h). qPCR was performed using SYBR Premix Ex Taq (Takara Bio, Inc.) on an ABI PRISM 7500 Sequence Detection System (Applied Biosystems; Thermo Fisher Scientific, Inc.). Primers for cDNA amplification were as follows:

miR-17-5p forward, 5’-AGCTGGACCAAAGTG CTTACAGT-3’, and

Figure 1.

The profile of m6A-miRNA regulators. (A)Differential expression of 23 m6A regulators; (B) the correlations between m6A-miRNAs and different expressed m6A regulators. * $P<$ 0.05, ** $P<$ 0.01,*** $P<$ 0.001.

universal reverse primer, 5’-ATCCAGTGCAGGGT CCGAGG-3’; RT primer, 5’-GTCGTATCCAGTGCAGGGTCCG AGGTATTCGCACTGGATACGACCTACCT-3’; miR-182-5p forward, 5’-CGGCGGAGACAAGA GGAGC-3’, reverse, 5’-CAGGGTGTCAGGAGGTGGGAG-3’. RT primer, 5’-GTCGTATCCAGTGCAGGGTCCG AGGTATTCGCACTGGATACGACAGTGTG-3’; U6 forward, 5’-GCTTCGGCAGCACATATACTAA AAT-3’, and reverse, 5’-CGCTTCACGAATTTGC GTGTCAT-3’;

The reaction mixtures were denatured at 95 ${}^{\circ }$ C for 3 min, followed by 40 cycles of 95 ${}^{\circ }$ C for 10 sec and 60 ${}^{\circ }$ C for 30 sec. The relative expression levels of miRNA and mRNA were normalized to that of U6. The relative expression were calculated and normalized using the 2- $\mathrm{\Delta }\mathrm{\Delta }$ Cq method. All reactions were conducted in triplicate.

Figure 2.

Survival analysis of m6A-miRNAs. (A) Forest map of cox regression based on the m6A-miRNAs and HCC overall survival; (B) Consensus cluster analysis of the m6A-miRNAs; (C) Kaplan-Meier curve of HCC overall survival and clusters.

2.13 RNA isaolation and RNA m6A dot blot assay

Total RNA was isolated as described above and the poly (A) RNA (600, 300, and 150 ng) dots were subjected to RNA m6A dot hybridization analysis on nylon membranes (GE Healthcare, USA) and then probed with m6A antibody. The membranes were then UV cross-linked, blocked, incubated with m6A antibody (1:1000, Abcam, USA) overnight at 4 ${}^{\circ }$ C. Then, incubated with HRP-conjugated goat anti-mouse IgG. Finally, they were dyed with methylene blue and photographed (Bio-Rad, USA).

2.14 Statistical analysis

Student’s $t$ -tests were performed to compare the differences between two groups, and the Kruskal-Wallis test was used for comparisons between two or more groups. OS refers to the interval from the date of diagnosis to the date of death and were compared between different groups by Kaplan-Meier analysis and the log-rank test. The prognostic ability of the risk model, clinicopathological factors and nomogram for 5-year OS was assessed by time-dependent receiver operating characteristic (ROC) curve and the area under the curve (AUC) values through the “timeROC” package. All statistical analyses were conducted with R software version 4.1.2.

3. Results

3.1 Identification of m6A-miRNAs

A total of 20 m6A regulators were differentially expressed (Fig. 1A). Of these, 15 were upregulated in the HCC tumor tissues, and 5 were downregulated in the HCC tumor tissues. Further correlation analysis identified 48 m6A-miRNAs (Fig. 1B, Supplementary Table 2).

3.2 Characteristics of the prognostic value of m6A-miRNAs

A total of 17 m6A-miRNAs were associated with HCC prognosis and screened by Cox regression analysis (Fig. 2A). Based on the expression level of the 17 m6A-miRNAs, the $k=$ 3 was selected as the most stable clusters when $k$ was set from 2 to 9 (Fig. 2B). Cluster 3 had the highest survival probability (Fig. 2C).

Figure 3.

Construction of m6A-miRNAs signature for prognosis. (A) and (B) M6a-miRNAs screening for model construcyopm by LASSO regression; (C) Risk score correlated with OS; (D) Distribution of the risk score between the HCC patients; (E) Kaplan – Meier curve of OS for HCC patients based on the risk score; (F) Time-dependent ROC curve for evaluating the predictive efficacy of the risk score.

Figure 4.

Prognostic performance of clinical features and m6A-miRNAs signature. (A) Kaplan-Meier curve of HCC overall survival and clinical stages; (B) Time-dependent ROC curve for evaluating the predictive efficacy of the clinical stage; (C) A nomogram was constructed based on age, tumor stage, risk score; (D) Time-dependent ROC curve for evaluating the predictive efficacy of the clinical stage conbined risk socre.

3.3 Construction of the m6A-miRNA signature for prognosis

Then, a signature of nine m6A-miRNAs was identified by LASSO regression analysis to predicted the prognosis of HCC. (Fig. 3A and 3B, Supplementary Table 3). The risk score was calculated based on the coefficients and expression levels, and we divided the risk scores into high-risk and low-risk based on the median (Fig. 3C). According to the follow-up outcomes, we divided patients into the deceased group and the survival group. The distribution of the risk scores showed that the deceased group had a higher risk score (Fig. 3D, $P<$ 0.0001). Survival analysis indicated that patients in the high-risk score group demonstrated worse survival (Fig. 3E). The AUC of the m6A-miRNAs signature for predicting prognosis was 0.772 (Fig. 3F), suggesting that this risk model could predict the prognosis.

Figure 5.

The tumor immune microenvironment characteristics of the prognostic related m6A-miRNAs clusters and risk score. (A and B) Distribution of m6A-miRNA clusters in Tumor immune microenvironment A was immune cell infiltration; B was scores; (C and D) Distribution of risk scores in tumor immune microenvironment immune cells infiltration and scores. * $P<$ 0.05, ** $P<$ 0.01,*** $P<$ 0.001.

3.4 The prognostic performance of the clinical features and m6A-miRNAs signature

To explore the prognostic value of the clinical information from HCC patients, we collected clinical information, including sex, age, tumor stage, BMI, neoplasm histological grade and viral hepatitis serology. Further survival analysis showed that the tumor stage was associated with the prognosis of HCC (Fig. 4A, $P<$ 0.0001) and the AUC of the tumor stage was 0.631 (Fig. 4B). No associations were found between prognosis and the other clinical information(Supplementary Table 4). The nomogram showed that an older patient with an advanced tumor stage and a higher risk score had a higher HCC mortality risk (Fig. 4C). Moreover, the risk score model could improve the prognosis performance of the clinical features, as the AUC of tumor stage was 0.771 (Fig. 4D).

Figure 6.

Association of the checkpoints (PD1 and PD-L1) with prognostic related m6A-miRNAs clusters and risk score. (A) The expression of PD1 in the three cluster; (B) The expression of PD-L1 in the three cluster; (C) The expression of PD1 in the two risk score group; (D) The expression of PD-L1 in the two risk score group; * $P<$ 0.05, ** $P<$ 0.01, *** $P<$ 0.001, NS $>$ 0.05.

3.5 The effect of the tumor immune microenvironment characteristics on the prognostic value of the related m6A-miRNAs clusters and risk scores

Figure 7.

Functional analysis of m6A-miRNAs signature. (A) KEGG; (B) GO.

Figure 8.

Validation of m6A-miRNAs signature for HCC prognostic prediction in GEO dataset. (A) ROC curve for evaluating the predictive efficacy of the m6A-miRNAs signature in HCC portal vein tumor thrombosis from GSE76903 and GSE6857; (B) Distribution of the m6A-miRNAs signature risk score in liver normal tissee and portal vein tumor thrombosis from GSE76903; (C) ROC curve for evaluating the predictive efficacy of the m6A-miRNAs signature in HCC survival from GSE6857; (D) ROC curve for evaluating the predictive efficacy of the m6A-miRNAs signature in HCC metastasis from GSE6857.

To investigate the effect of the tumor immune microenvironment on the prognostic related m6A-miRNAs, we evaluated the immune cells of the infiltration and immune scores of HCC tissues. The distribution of the three clusters were different among all the immune infiltration cell subsets (Fig. 5A). Then, to confirm the distribution of the three clusters, we calculated the immune scores using ESTIMATE. The distribution of the three clusters were different in tumor microenvironment (Fig. 5B) as well.

Furthermore, we conducted a correlation analysis to evaluate the relationship between risk scores and immune infiltration. The results showed that CD8+ T cells, CD56+ cells NK cells, $\gamma \delta$ T cells, and mast cells were not associated with risk scores. However, other immune infiltration cells were positively correlated with risk scores (Fig. 5C). The distribution of risk scores was different in the tumor microenvironment (Fig. 5D).

Figure 9.

Meta-analysis of the prognosis m6A-miRNAs. (A) forest plot of miR-17-5p; (B) forest plot of miR-182-5p; (C) forest plot of miR-652-3p.

3.6 Association between the checkpoints (PD-1/PD-L1) and the prognosis-related m6A-miRNAs clusters and risk scores

We next explored the potential association between the m6A-miRNAs and the checkpoints. We found that the expression level of PD-1 was different in the different m6A-miRNAs expression clusters (Fig. 6A and 6B). Patients with high PD-1/PD-L1 expression level had a high risk score (Fig. 6C and 6D).

3.7 Functional analysis of the m6A-miRNAs signature and validation of the m6A-miRNAs signature for prognosis, recurrence and metastasis prediction

To explore the biological behaviors regulated by these miRNAs, we performed the KEGG and GO enrichment analyses. Those miRNAs were involved in pathways involved in cancer, the PI3K-Akt signaling pathway and DNA-binding transcription activity (Fig. 7A and 7B).

Next, We conducted logistical regression analysis of the 9-m6A-miRNAs and the risk score was summation of the product of nine miRNA expression levels and corresponding regression coefficients. we applied the nine-m6A-miRNA model for the analysis of an additional 20 HCC patients with portal vein tumor thrombosis (PVTT), which is a hallmark event of tumor poor prognosis from the GSE76903 dataset, to validate the predictive ability. Normal tissues were regarded as no recurrence events, and the PVTT tissues were regarded as recurrence events. The AUC of the validation cohort was 0.788, and the distribution of the risk score indicated that PVTT tissues were associated with higher risk scores (Fig. 8A and 8B).

Then we applied the nine-m6A-miRNA model for a larger dataset with 240 HCC patients, only eight-m6A-miRNAs were included for miR-652-3p was not included in the platform (GSE6857). The AUC of the validation cohort was 0.646 for prognosis prediction, The AUC was 0.774 for metastasis prediction. (Fig. 8C and 8D).

3.8 Meta-analysis of the prognostic m6A-miRNAs

To validate the prognostic value of the nine-m6A-miRNA model, we conducted a meta-analysis. Four hospital-based studies [23, 24, 25, 26] were included, and only three miRNAs were investigated, miR-17-5p, miR-182-5p and miR-652-3p (Table 1). The results of the current studies from TCGA were included. The pooled results revealed that miR-17-5p (HR: 1.58, 95%CI: 1.04–2.40, Fig. 9A) and miR-182-5p were associated with worse OS for HCC patients (HR: 1.17, 95% CI: 1.03–1.32; Fig. 9B).

Table 1

Characteristics of the studies included in the meta analysis

Prognostic miR	First author, date	Country	Cut-off	Sample size	Sample type	PMID	Prognostic test	HR	95%CIL	95%CIH	$P$	Resource
miR-17-5p	Chen(2012)	China	Median	120	Tissue	22583011	Multivariate	1.790	1.140	2.980	0.042	Hospital
miR-17-5p	Zheng(2013)	China	Median	96	Serum	23108086	Multivariate	2.192	1.024	4.691	0.043	Hospital
miR-182-5p	Matsui(2021)	Japan	Median	111	Tissue	34038612	Multivariate	2.415	1.268	4.599	0.007	Hospital
miR-652-3p	Chi(2021)	China	Median	126	Tissue	34608826	Multivariate	3.249	1.323	7.979	0.010	Hospital

Figure 10.

The expression of miR-17-5p and miR-183-5p. (A) and (B) The expression of miR-17-5p and miR-183-5p in HCC and tumor adjacent were determined by RT-PCR in seven HCC patients. (C) The expression of miR-17-5p in 32 paired HCC and adjacent tissues from GES227383.

3.9 The expression of miR-17-5p is up-regulated in HCC tissues

The expression of miR-17-5p and miR-183-5p in HCC and tumor adjacent were determined by RT-PCR in seven HCC patients (Fig. 11A and B), miR-17-5p was upregulated in HCC tissues ( $P<$ 0.05, Fig. 10B) and validated the result in 32 paired HCC tissues from GEO ( $P<$ 0.05, Fig. 10C).

3.10 miR-17-5p promotes HCC cell growth and invasion in vitro

We applied three HCC cell lines (Hep3B, HepG2 and Luh7) to evaluate the efficiency of overexpression and silencing, the expression were of miR-17-5p were detected by RT-PCR after the miR-mimics and miR-inhibitors were transfected (Fig. 11A). The cell viability promoted after the miR-mimics were transfected and repressed by miR-inhibits (Fig. 11B).Transwell assay was applied to evaluate invasive and migrate capabilities, the results showed that miR-17-5p-mimic/inhibitor affected cell migration and invasion ( $P<$ 0.05, Fig. 11C). Then we assayed the m6A dot blot assays of Luh7 cells with inhibitor and mimic of miR-17-5p. It was obvious that the methylation of RNA increased significantly after miR-17-5p inhibitor, while decreased after miR-17-5p mimic (Fig. 12).

Figure 11.

miR-17-5p promotes HCC tumorigenesis in vivo. A.RT-qPCR analysis of the expression levels of miR-17-5p-inhibitor/mimic/control groups. B. HepG2, Hep3B and Luh7 cells were cultured after miR-17-5p-inhibitor/mimic/control transient transfection, CCK8 assay was used to access the viability of HCC cells. C. The effect of miR-17-5p on cell migration and invasion in LUH7 cells.

Figure 12.

m6A dot blot assays the detection of Luh7 and HepG2 cells with inhibitor of miR-17-5p and NC control. It was obvious that the methylation of RNA increased significantly after miR-17-5p inhibitor.

4. Discussion

In this study, we identified a nine-m6A-miRNA signature that is associated with HCC prognosis and validated this signature associated with HCC recurrence. This study is the first to systematically evaluate the prognostic performance of m6A-miRNAs. Our results showed that patients with high risk scores according to the m6A-miRNAs signature have decreased cancer survival rates and increased HCC recurrence rates. These m6A-miRNAs may affect tumor progression and metastasis. Then, we further investigate the immune-related characteristics of the prognostic m6A-miRNAs and found that the three clusters were differentially distributed in the TME. Moreover, functional analysis revealed that these miRNAs are involved in cancer pathways. The meta-analysis showed that miR-17-5p and miR-182-5p were associated with HCC prognosis. We further verified that expression of miR-17-5p was upregulated in HCC tissues. Based on the transient transfection of miR-17-5p-mimic or miR-17-5p-inhibitor in HCC cell lines, we observed that miR-17-5p increased HCC growth and invasion in vitro and in vivo. Moreover, the expression level of miR-17-5p was negative correlations with the methylation of total RNA.

Many studies have reported that miRNAs have good diagnostic value for cancers [27, 28], including HCC [29]. Recently, Zhang et al, found that m6A-targeted microRNAs in serum had potential use in cancer diagnosis [7], and 18 m6A-miRNAs exhibited a good efficacy in the diagnosis of HCC. Among those m6A-miRNAs, miR-17-5p was identified in the current study as a prognostic marker for HCC patients. It has reported that miR-17-5p could suppress HCC metastasis by blocking the HGF/ERBB3-NF- $\kappa$ B positive feedback loop. This provides a possible new strategy for preventing metastasis after HCC resection [30]. Moreover, an increasing number of studies have shown that miRNAs could be used as prognostic markers. Xi et al. constructed a three-miRNA signature for use as a prognostic biomarker for HCC patients based on TCGA database [31]. Different from the current study, the three miRNAs were identified by a differential expression strategy, and that study did not focus on their role in HCC recurrence.

Among the miRNAs in the current nine-m6A-miRNA signature, most miRNAs were involved in cancer pathways, as the functional enrichment analysis showed that the potential target genes of these miRNAs were involved in cancer-related pathways. MiR-148-3p is a prognostic marker for many types of tumors, such as prostate cancer [32], acute myeloid leukemia [33] and colorectal cancer [34]. Further mechanistic study indicated that miR-148a-3p was sponged by circANKS1B and then increased the expression level of the transcription factor USF1. This ultimately activates TGF- $\beta$ 1/Smad signaling, which could promote EMT [35]. HCC was inhibited by miR-182-5p, as it regulates liver cancer cell invasion, adhesion, migration and proliferation through the Wnt/ $\beta$ -catenin signaling pathway [36]. It was reported that TRERNA1 acted as a ceRNA and could regulate NRAS proto-oncogene (NRAS) through sponging miR-22-3p. Then, it could increase the expression of TRERNA1 through the activation of the RAS/Raf/MEK/ERK signaling pathway [37]. Chi et al. reported that miR-652-3p could take part in the tumor progression of HCC [26]. Cell viability and anchorage were promoted by miR-9-5p in cervical cancer cell lines, and it acted as a tumor suppressor [38]. miR-99b-3p directly targets PCDH19, resulting in HCC metastasis and proliferation [39]. Li et al. reported that miR-212-3p was directly upregulated by ZLM-7 in breast cancer. Additionally, VEGFA expression was inhibited by miR-212-3p through Sp1, and then, the angiogenesis and progression of breast cancer were inhibited [40]. Many studies have reported that the small non-coding RNA let-7b-5p is a biomarker for cancer prognosis in cancers such as pancreatic cancer [41], breast cancer [42] and colon cancer [43]. A systematic regulatory network analysis showed that let-7b-5p might promote the pathogenesis of non-small cell lung cancer [44].

Subsequent analysis on the TME and m6A-miRNAs revealed that the three consensus clusters of the prognostic m6A-miRNAs were differentially distributed in the immune infiltration cells. Then, the patients with high risk scores may have higher stromal scores and immune scores and lower tumor purity. Previous studies have shown that immune infiltration is associated with the prognosis of tumors [45, 46]. Interestingly, we found that the consensus clusters were associated with the expression of PD1, and high risk scores were associated with high expression levels of PD1/PD-L1, as well as actived CD4+ T cell. Tfh cells were the promoter of anti-tumor immune responses and it was suffered immunosuppression, when the expression of PD-1 was upregulated [47]. The findings could reveal the poor prognosis of HCC patients in the high-risk group.

Many studies have reported that non coding RNAs may affect tumor progression and metastasis such as long non-coding RNAs(lncRNAs), Zhong et al. found that m6A related lncRNAs could forecast the survival of breast cancer, colorectal cancer and HCC [48, 49, 50]. In the current study, we firstly explored the m6A-miRNAs in HCC prognosis and TME. In conclusion, we first found a nine-m6A-miRNA signature with predictive value for HCC prognosis and recurrence. Furthermore, the immune related characteristics of the prognostic m6A-miRNAs was determined. MiR-17-5p may play an oncogene role in HCC by interactions with methylation modification. The limitations existed in this study. First, we validated the model of m6A-miRNAs signature in two cohort from GEO, further larger sample size, multi-ethnic and different etiologies cohort is needed to establish. Second, further mechanistic investigation for the m6A-miRNAs are needed in subsequent studies.

Author contribution

Conception: Jing Han, Longbiao Zhu and Jifeng Feng.

Interpretation or analysis of data: Yaqian Liu and Jing Han.

Preparation of the manuscript: Renrui Zou, Yaqian Liu and Sangsang Qiu.

Revision for important intellectual content: Longbiao Zhu, Ya Lu, Yan Chen and Hui Yu.

Supervision: Hangju Zhu and Wenbo Zhu.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-230263.

sj-pdf-1-cbm-10.3233_CBM-230263.pdf - Supplemental material

Supplemental material, sj-pdf-1-cbm-10.3233_CBM-230263.pdf

sj-xlsx-1-cbm-10.3233_CBM-230263.xlsx - Supplemental material

Supplemental material, sj-xlsx-1-cbm-10.3233_CBM-230263.xlsx

sj-xlsx-2-cbm-10.3233_CBM-230263.xlsx - Supplemental material

Supplemental material, sj-xlsx-2-cbm-10.3233_CBM-230263.xlsx

sj-xlsx-3-cbm-10.3233_CBM-230263.xlsx - Supplemental material

Supplemental material, sj-xlsx-3-cbm-10.3233_CBM-230263.xlsx

sj-xlsx-4-cbm-10.3233_CBM-230263.xlsx - Supplemental material

Supplemental material, sj-xlsx-4-cbm-10.3233_CBM-230263.xlsx

Footnotes

Acknowledgments

This work was supported by the National Natural Science Foundation of China (81702686, 81602920), Natural Science Foundation of Jiangsu Province (BK20211383) and the Young Talents Program of Jiangsu Cancer Hospital, Jiangsu.

Conflict of interest

The authors declare that they have no conflicts of interest.