FPKM-masked sequencing data of 374 primary hepatocyte malignant tissues and 50 normal tissues were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov) database and annotated with human gene annotation files from the Ensembl official website (http://asia.ensembl.org/index.html; GRCh38.103. chr.gtf.gz). GTF files (GRCh38.103. chr.gtf.gz) were used to annotate and distinguish the mRNAs and lncRNAs for subsequent analysis. Corresponding clinicopathological data were downloaded concomitantly. Patients with pathologically confirmed HCC and complete clinical information were included in this study; those with a survival time inferior or equal to 30 days were excluded to avoid any bias originating from demise due to unrelated factors. A total of 342 HCC cases were included in the subsequent analysis.

By systematically searching and reviewing representative m6A-related literature, 23 universally acknowledged m6A regulators, namely, 8 methyltransferases, 2 demethylases, and 13 binding proteins, were selected. The details of these genes ^{[14, 31–35]} are listed in Supplementary Table S1. We calculated the correlation coefficient between the expression levels of m6A regulators and the expression level of lncRNA and selected the one with |Pearson r| > 0.4 and p < 0.001 as the m6A regulator-related lncRNA (m6ARlnc). Then, the expression matrix of the m6ARlnc genes was extracted. A flow chart of this study is illustrated in Supplementary Figure S1.

To identify prognostic m6ARlncs, univariate Cox regression analysis was conducted to screen the differentially expressed lncRNAs between the tumor and normal groups, with significance set at p ≤ 5 × 10⁻⁵, using the limma R package.

Based on the differentially expressed prognostic m6ARlncs in tumor and normal tissues, unsupervised consensus clustering was performed using the ConsensusClusterPlus R package. Based on the similarity between the expression levels of prognostic m6ARlncs and the proportion of fuzzy clustering measurements, the cluster number with the best stability was determined. In addition, CIBERSORT analysis was performed to determine the proportion of 22 types of immune cells in each sample. The ESTIMATE algorithm was used to evaluate the tumor microenvironment (TME) score (including stromal cell, immune cell, and comprehensive scores).

To fully understand the predictive value of prognostic m6ARlnc, 342 HCC cases were randomly divided into a training cohort (n = 226) and a validation cohort (n = 116) in a 2:1 ratio. In this study, all subsequent analyses were performed in the training cohort, except for model validation, hierarchical analysis of the model, and enrichment analysis of the signature, which were performed on the entire cohort. First, the glmnet, survminer, and survival R packages were used to conduct Lasso-Cox regression analysis to construct a prognosis prediction model of patients with HCC by screening the prognostic m6ARlncs. The model was optimized by cross-validation, and lncRNAs used to construct the prognostic signature were screened. A survival plot was constructed using the survival R package. The risk score of each patient was calculated using the following formula: R i s k   s c o r e = ∑ i = 1 n C o e f i × X i

Coef _i represents the coefficient of the corresponding m6ARlnc from the Cox regression model, and X _i is the expression level of m6ARlnc. The X-tile software was used to determine the optimal cutoff value of the risk score for the training cohort, and the same cutoff value was used in the validation and entire cohorts. Thus, all patients were divided into high- and low-risk groups. Second, Kaplan–Meier (K-M) survival curve analysis was conducted to compare the survival rates between groups using the survival R package.

The receiver operator characteristic (ROC) curves for 1, 3, and 5 years were constructed using the timeROC R package to evaluate the predictive accuracy of this signature, and the areas under the curve (AUC) were calculated. In addition, to confirm that the signature has consistent predictions across patients with different clinical characteristics such as sex, age, and tumor status, we conducted subgroup stratification survival analysis based on clinicopathological features using the K-M plot.

To improve the clinical application of the prognostic signature, we constructed a comprehensive prediction model combining the signature risk score and clinicopathological factors. First, we conducted a univariate Cox regression analysis, and variables with p < 0.0001 were subsequently included in the multivariate Cox regression analysis to screen independent factors in predicting the prognosis of patients with HCC. Second, we evaluated whether there were differences among clinicopathological factors, classification characteristics of prognostic m6ARlnc, stromal cell score, immune cell score, comprehensive score, and risk score between the differential expression of lncRNAs in the prognostic signature. The stromal cell, immune cell, and comprehensive scores were calculated using the ESTIMATE algorithm. Finally, a nomogram was established, according to the independent prognostic factors, to predict overall survival (OS) of patients using the rms R package.

To identify potential signaling pathways related to high- and low-risk patients determined by prognostic signature, we performed a gene set enrichment analysis (GSEA) using the GSEA software (version: 4.1.0). The gene set c2. cp.kegg.v7.4. symbols.gmt was downloaded from Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, and the gene set used for Gene Ontology (GO) enrichment analysis was c5. go.bp.v7.3. symbols.gmt.

To explore the regulatory relationship between lncRNA and m6A regulators, constructing a ceRNA network is a good choice. The miRcode database (http://www.mircode.org/) was used to predict miRNAs that interact with prognostic m6ARlncs. To identify the target mRNAs, the relationship between miRNA and target mRNA was characterized using miRDB (http://www.mirdb.org/), miRTarBase (https://mirtarbase.cuhk.edu.cn/), and TargetScan (http://www.targetscan.org/). In addition, the Cytoscape (version 3.6.0) software was used to visualize the lncRNA–miRNA–mRNA ceRNA network. Subsequently, GO and KEGG enrichment analyses of target mRNAs were performed using the clusterProfiler R package to explore the potential signaling pathways of m6ARlncs.

To predict the sensitivity of high- and low-risk patients to chemotherapeutic agents, we used the pRRophetic (https://github.com/paulgeeleher/pRRophetic) R package to predict the half-maximal inhibitory concentration (IC50) of chemotherapeutic agents for each patient. The R package was based on pre-treatment gene expression and drug sensitivity data of cancer cell lines to predict the chemotherapeutic response (Geeleher et al., 2014). Thirty-three immune checkpoints (Supplementary Table S2) were used to evaluate the performance of the model for immunotherapy and to analyze the differences in immune checkpoint expression levels between high- and low-risk groups.

The Wilcoxon rank-sum test was used to compare the differences between two groups of continuous variables, and the Kruskal–Wallis test was used to examine the difference in TEM scores between the different subtypes. The chi-square test was used to examine the expression differences between lncRNAs. The K-M curve, combined with the log-rank test, was used to evaluate the OS of the different groups. Independent predictive factors were identified using Cox proportional hazards regression analysis. Based on the results of multivariate Cox analysis, a nomogram was established and evaluated for consistency using the C-index and calibration curve. All statistical analyses were carried out using the R software (version 3.6.3). Unless indicated otherwise, p < 0.05 was considered statistically significant.

A total of 644 lncRNAs that were significantly correlated with 23 m6A regulators were identified in the TCGA dataset and defined as m6ARlnc. The co-expression network between these 23 m6A regulators and the related lncRNAs is shown in Figure 1A. We conducted a univariate Cox regression analysis, and the results indicated that out of these 644 lncRNAs, 31 were associated with the OS of patients with HCC (Figure 1B). We further compared the expression levels of these 31 lncRNAs between tumor and normal tissues and found that the tumor group presented a significantly higher expression level (Figures 1C,D).

To identify heterogeneity, unsupervised consensus clustering was performed on all samples based on the expression levels of prognostic m6ARlncs. We determined that k = 3 had the best clustering stability among the range k = 2 to 9 (Figures 2A–D). Thus, patients with HCC were divided into three subgroups: Cluster1 (n = 80), Cluster2 (n = 243), and Cluster3 (n = 19). To evaluate the prognostic value of m6ARlncs, survival analysis was performed on the different subgroups; significant differences were observed in the clinical outcomes. The OS of Cluster2 was the highest, while that of Cluster3 was the lowest (Figure 2E). The heatmap shows the differences in the expression of prognostic m6ARlncs in different subgroups, and the expression levels of most prognostic m6ARlncs were elevated in Cluster3; contrastingly, their expression levels were lowest in Cluster2. The differences in the expression of prognostic m6ARlncs in concomitance of various clinicopathological factors are also described in the above-mentioned heatmap. There was no difference in the expression of these lncRNAs among the clinicopathological factors (Figure 2F). Our results indicate a close relationship between the aggregation subgroup defined by the expression levels of prognostic m6ARlncs and the heterogeneity of patients with HCC.

Increasing evidence shows that m6ARlncs are involved in the tumor immune microenvironment. CIBERSORT analysis showed that T cells regulatory (p = 0.014) and B cells naive (p = 0.011) were significantly enriched among the three subgroups (Figure 3A). Then, the stromal cell score, immune cell score, and comprehensive score of each sample were obtained using the ESTIMATE algorithm. Results revealed that the stromal cell score of Cluster2 was higher than that of Cluster1 (p < 0.001) and Cluster3 (p = 0.028), and the comprehensive score of Cluster2 was also higher than Cluster1 (p = 0.042), indicating that Cluster2 contained more immune-related cells and tumor cells of lower purity (Figures 3B–D).

Nine lncRNAs (LINC01224, C2orf27A, AL158166.1, AL117336.2, MKLN1-AS, AL031985.3, LINC01138, POLH-AS1, and AL442125.2) were finally selected to construct a prognostic signature (m6A-9LPS) using LASSO Cox regression (Figures 4A,B). The risk score for each HCC patient was calculated based on the expression levels and coefficients of these nine lncRNAs, using the following equation (Table 1): Risk   score = ( 0.3222 × LINC 01224   expression   level ) + ( 0.0363 × C 2 orf 27 A   expression   level ) + ( 0.2281 × AL 158166 . 1   expression   level ) + ( 0.1265 × AL 117336 . 2   expression   level ) + ( 0.4807 × MKLN 1 − AS   expression   level ) + ( 0.3838 × AL 031985 . 3   expression   level ) + ( 0.0123 × LINC 01138  expression   level ) + ( 0.6103 × POLH − AS 1   expression   level ) + ( 0.5916 × AL 442125 . 2   expression   level )

Patients in the training cohort were assigned to the high- or low-risk group according to the best cutoff value of 2.9, which was determined via X-tile software analysis (Figures 4C,D). When analyzing the survival time and status of patients with HCC according to the distribution of risk scores, patients in the high-risk group showed a worse prognosis than those in the low-risk group. The nine lncRNAs were differentially expressed in the high- and low-risk groups, with a higher expression level in the high-risk group (Figures 4E–G). Patients in the validation cohort were also divided into high- and low-risk groups using the same cutoff value of the training cohort (Figures 4H–J). Furthermore, the prognosis of the two groups was compared using the K-M curve, and the results revealed that patients with a high-risk score had a poor prognosis both in the training and validation cohorts (p < 0.001, Figures 4K,L). In the training cohort, the 3-year survival rates of the high- and low-risk groups were 31.8 and 72.8%, respectively; in the validation cohort, these rates were 29.9 and 80.2%, respectively.

In addition, we evaluated the predictive ability of the established m6A-9LPS signature in the training, validation, and entire cohorts. Regardless of the ROC curve analysis in the three cohorts, m6A-9LPS showed strong predictive capabilities (Figures 5A–C). Moreover, the AUC values of m6A-9LPS were higher than those of single signatures in the training, validation, and entire cohorts, which was confirmed by time-dependent ROC analysis (Figures 5D–F). These results suggest that m6A-9LPS has a strong ability to predict the prognosis of patients with HCC.

To confirm the stability of this model in predicting the prognosis of patients with HCC under different clinicopathological conditions, we performed a stratified analysis on the entire cohort of patients divided into different subgroups based on their different clinical characteristics. Survival analysis revealed significant differences in prognosis between high-risk and low-risk groups classified by the m6A-9LPS model, regardless of sex, age, grade, or AJCC staging. In patients over 65 years of age (n = 126) and patients aged 65 or younger (n = 216), m6A-9LPS could accurately divide the high- and low-risk groups (p < 0.001) (Figures 6A,B); similarly, the male (n = 233, p < 0.001) and female (n = 109, p = 0.02) patients were also accurately divided into high- and low-risk groups (Figures 6C,D). In addition, with regard to grade and AJCC staging, the K-M curve showed that the low-risk group had a better OS (Figures 6E–P). However, for patients with grade G4 (n = 12, p = 0.79), stage N1 (n = 3), stage M1 (n = 3), and stage IV (n = 3), there was no significant difference in prognosis between the high-risk and low-risk groups (Supplementary Figure S2). These results suggest that m6A-9LPS may be a powerful prognostic signature for patients with HCC, and that it is not affected by changes in traditional clinicopathological features.

Independent prognostic analysis helped us understand whether m6A-9LPS could be independent of other clinicopathological factors that may affect the prognosis of patients with HCC. Univariate and multivariate Cox regression analyses showed that AJCC staging and m6A-9LPS risk scores were significantly correlated with the prognosis of patients (both p < 0.001). This indicates that, the m6A-9LPS risk score is an independent prognostic indicator of HCC (Figures 7A–D). The heatmap shows that the nine lncRNAs were upregulated in the high-risk group. Except for the cluster group (p < 0.001), there was no significant difference in the expression of the nine lncRNAs between the high- and low-risk groups (Figure 8A).

We then constructed a nomogram to predict the OS of patients with HCC based on all the significant prognostic factors from the multivariate Cox regression analysis. The nomogram was used to evaluate the 1-year, 3-year, and 5-year OS rates. The total score was calculated based on the scores of all the variables in the nomogram. The OS rate of each patient could be predicted by drawing a vertical line from the total score to the survival prediction axis (Figure 8B). In addition, the C-index of the nomogram was 0.726 (95% CI, 0.663–0.788). The results of the calibration curve analysis showed that, when predicting the OS of patients with HCC at 1 year, 3 years, and 5years, the nomogram had precise predictive ability, and the OS rate was close to reality in the training, validation, and entire cohorts (Figures 8C–E).

We constructed a ceRNA network based on prognostic m6ARlncs to illustrate how m6ARlnc can regulate target mRNA expression in HCC via sponging miRNA. Four prognostic m6ARlncs were extracted from the miRcode database, and 98 pairs of interactions between the four lncRNAs and 11 miRNAs were identified. Based on the results of 11 miRNAs, 352 target mRNAs were identified. Interestingly, these mRNAs contained the m6A regulator FMR1 (Edupuganti et al., 2017), which confirmed that in HCC, m6ARlnc can act as a ceRNA by targeting m6A regulators through miRNA. Finally, we constructed a ceRNA network of lncRNA–miRNA–mRNA containing four lncRNAs, 11 miRNAs, and 352 mRNAs (Figure 9A). Enrichment analysis of 352 target mRNAs showed that these genes were correlated with enriched Ras protein signal transduction, negative regulation of G1/S transition of mitotic cell cycle, negative regulation of cell cycle G1/S phase transition, cell cycle arrest, and positive regulation of cell morphogenesis involved in differentiation (the top five terms of GO biological processes). The top five KEGG signaling pathways were in hepatitis B, pancreatic cancer, endocytosis, axon guidance, and Kaposi sarcoma-associated herpesvirus infection (Figures 9B,C). This information may provide a reference for further exploration of the function of m6ARlnc in HCC.

To investigate the biological processes and pathways underlying the molecular heterogeneity between the high- and low-risk m6A-9LPS divided subgroups, we performed KEGG and GO enrichment analyses on the entire cohort. In the GO enrichment analysis, we selected the biological process (BP), combined with the results of the KEGG enrichment analysis, to explore the main biological processes involving lncRNAs in patients with HCC. In the KEGG enrichment analysis, the cell cycle, DNA replication, fatty acid metabolism, homologous recombination, and mismatch repair pathways were significantly enriched in the tumors of the m6A-9LPS-defined high-risk group. In the low-risk group, the significantly enriched pathways were mainly related to various metabolic processes (including primary bile acid biosynthesis, tryptophan metabolism, and acetone metabolism), spliceosomes, and the degradation of certain amino acids (valine, leucine, and isoleucine) (Figure 10A). The results of GO enrichment analysis showed that in the high-risk group, cell cycle DNA replication, DNA-dependent replication, mRNA export from the nucleus, ncRNA export from the nucleus, and tRNA transport played important roles, while no significant biological processes were enriched in the low-risk group (Figure 10B).

The IC50 of six commonly used chemotherapeutic agents was predicted in high- and low-risk groups for HCC, which was based on the algorithm provided in the pRRophetic R package. The IC50 values of axitinib, dasatinib, erlotinib, sorafenib, and sunitinib were higher in the high-risk group, indicating that patients in the low-risk group were more sensitive to these five drugs; at the same time, the IC50 of gemcitabine was higher in low-risk patients, which meant that high-risk patients were more sensitive to gemcitabine (Figures 10C–H). In addition, the relationship between the gene expression levels of immune checkpoints and the risk group was further analyzed. The expression levels of 33 immune checkpoints in the high-risk group were significantly higher than those in the low-risk group (p < 0.001, Figure 10I).