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

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 = ∑ i = 1 n β i ⁢ 𝑒𝑥𝑝 i , where β i represents the regression coefficient, and 𝑒𝑥𝑝 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 ∘ 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 ∘ 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 ∘ C for 24 h prior to transfection. Then, cells were transfected with corresponding oligonucleotides. After 24, 48 or 72 h, 10 μ l CCK-8 solution was added into each well of the plate. The plates were incubated at 37 ∘ 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 × 10 5 cells in 500 μ 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 ∘ C for 1 h), while the TaqMan MicroRNA Reverse Transcription kit (Takara Biotechnology Co., Ltd.) was used for miRNA detection (42 ∘ 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 OPEN IN VIEWER

The reaction mixtures were denatured at 95 ∘ C for 3 min, followed by 40 cycles of 95 ∘ C for 10 sec and 60 ∘ 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- Δ ⁢ Δ Cq method. All reactions were conducted in triplicate.

OPEN IN VIEWER

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 ∘ 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.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).

OPEN IN VIEWER

OPEN IN VIEWER

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, 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.

OPEN IN VIEWER

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, 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).

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

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, γ ⁢ δ 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).

OPEN IN VIEWER

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).