The m6A/m5C/m1A Regulated Gene Signature Predicts the Prognosis and Correlates With the Immune Status of Hepatocellular Carcinoma

RNA modification of m6A/m5C/m1A contributes to the occurrence and development of cancer. Consequently, this study aimed to investigate the functions of m6A/m5C/m1A regulated genes in the prognosis and immune microenvironment of hepatocellular carcinoma (HCC). The expression levels of 45 m6A/m5C/m1A regulated genes in HCC tissues were determined. The functional mechanisms and protein–protein interaction network of m6A/m5C/m1A regulated genes were investigated. The Cancer Genome Atlas (TCGA) HCC gene set was categorized based on 45 m6A/m5C/m1A regulated genes, and survival analysis was used to determine the relationship between the overall survival of HCC patients in subgroups. Cox and least absolute shrinkage and selection operator (LASSO) regression analyses were used to construct the risk model and nomogram for m6A/m5C/m1A regulated genes. The relationships between m6A/m5C/m1A regulated gene subsets and risk model and immune cell infiltration were analyzed using CIBERSORT. m6A/m5C/m1A regulated genes were involved in mRNA and RNA modifications, mRNA and RNA methylation, mRNA and RNA stability, and other processes. There was a statistically significant difference between cluster1 and cluster2 groups of genes regulated by m6A/m5C/m1A. The prognosis of cluster1 patients was significantly better than that of cluster2 patients. There were statistically significant differences between the two cluster groups in terms of fustat status, grade, clinical stage, and T stage of HCC patients. The risk model comprised the overexpression of YBX1, ZC3H13, YTHDF1, TRMT10C, YTHDF2, RRP8, TRMT6, LRPPRC, and IGF2BP3, which contributed to the poor prognosis of HCC patients. The high-risk score was associated with prognosis, fustat status, grade, clinical stage, T stage, and M stage and was an independent risk factor for poor prognosis in HCC patients. High-risk score mechanisms included spliceosome, RNA degradation, and DNA replication, among others, and high-risk was closely related to stromal score, CD4 memory resting T cells, M0 macrophages, M1 macrophages, resting mast cells, CD4 memory activated T cells, and follicular helper T cells. In conclusion, the cluster subgroup and risk model of m6A/m5C/m1A regulated genes were associated with the poor prognosis and immune microenvironment in HCC and are expected to be the new tools for assessing the prognosis of HCC patients.


INTRODUCTION
Liver cancer is one of the most prevalent cancers worldwide. Hepatocellular carcinoma (HCC) accounts for approximately 80% of liver cancer (1,2). Recent studies have demonstrated that targeted therapy and immunotherapy significantly affect HCC patient survival (3)(4)(5)(6). For example, Shimose et al. reported that sorafenib (SORA) improves overall survival (OS) in HCC patients. SORA can improve the prognosis of patients with unresectable HCC as first-line therapy (4). Nonetheless, the prognosis for HCC patients remains unsatisfactory. Therefore, it is extremely important to identify new targets or immunotherapies to improve the prognosis of HCC patients.
In eukaryotic messenger RNA (mRNA) regulation, N 6methyladenosine (m6A), N 1 -methyladenosine (m1A), and 5methylcytosine (m5C) modifications exist. Several studies have confirmed that m6A, m1A, and m5C regulated genes play important roles in m6A, m1A, and m5C modifications (7)(8)(9)(10)(11). Recent research has found that m6A, m1A, and m5C regulated gene expression levels are associated with tumor progression (11)(12)(13)(14)(15)(16). For example, the m6A-regulated gene methyltransferase 3 (METTL3), METTL14, and WTAP are involved in the initiation of the m6A modification process. METTL3 expression is elevated in endometrioid epithelial ovarian cancer (EEOC) tissues. The METTL3 overexpression levels correlate with the degree of malignancy and the OS of EEOC patients. Inhibiting METTL3 expression in TOV-112D and CRL-11731D cells attenuates cancer cell proliferation and migration and promotes apoptosis. METTL3 overexpression promotes EEOC progression by regulating m6A methylation (11). The m5C methyltransferase NSUN2 is significantly upregulated in gastric cancer and is predictive of a poor prognosis in gastric cancer patients. In vitro, NSUN2 promotes the proliferation, migration, and invasion of gastric cancer cells. Small ubiquitin-like modifier (SUMO)-2/3 promotes the oncogenic activity of NSUN2 by stabilizing NSUN2 (15). The risk model has been employed as a tool for determining the prognosis of cancer patients (17,18). Currently, the roles of m6A, m1A, and m5C regulated genes in HCC progression are not yet fully understood. Therefore, in this study, data from The Cancer Genome Atlas (TCGA) database were used to elucidate the biological functions and network involved in m6A/m1A/m5C regulated genes (8), and the roles of m6A/m1A/m5C regulated genes in the prognosis of HCC patients were identified. The risk model and nomogram of m6A/ m1A/m5C regulated genes were constructed to determine the prognosis of HCC patients.

MATERIALS AND METHODS
The Expression Levels of m6A/m1A/m5C Regulated Genes Gene expression data of HCC patients were extracted from TCGA database. A total of 50 normal liver tissue samples and 374 HCC tissue samples were included. The Writer, Reader, and Eraser genes of m6A/m1A/m5C were obtained from the literature (8), and the obtained m6A, m1A, and m5C regulated genes were merged and de-duplicated. Subsequently, the m6A/ m1A/m5C regulated gene expression data were retrieved in normal and HCC tissues, and the expression levels of m6A/ m1A/m5C regulated genes in HCC tissues were analyzed using the Limma package.
The Relationship Between m6A/m1A/m5C Regulated Genes The network of 45 m6A/m1A/m5C regulated genes in the STRING database was explored, and Cytoscape (version: 3.8.2) software was used to show the relationship between the proteinprotein interaction (PPI) network of m6A/m1A/m5C regulated genes.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analyses
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) methods are frequently used to assess the biological functions and the signaling mechanisms of polygenes (19). In this study, GO and KEGG methods were used to investigate the biological functions and signaling mechanisms associated with the 45 m6A/m1A/m5C regulated genes. In addition, the biological functions and signaling mechanisms involved in the abnormally expressed genes related to the 45 m6A/m1A/m5C regulated genes were determined using GO and KEGG methods.
Identification of m6A/m1A/m5C Regulated Gene Subgroups Using Consensus Clustering The "ConsensusClusterPlus" package was used to cluster 45 m6A/m1A/m5C regulated genes in 374 HCC tissues extracted from TCGA database into different subtypes. TCGA 374 HCC tissues were divided into cluster1 and cluster2 groups based on the optimal k value. Principal component analysis (PCA) was used to determine the differences between the two groups of patients. The Kaplan-Meier (K-M) survival analysis and logrank test were used to assess the relationship between survival and clinicopathological characteristics of HCC patients in the two cluster groups.

Cluster-Related Differentially Expressed
Genes Based on the m6A/m1A/m5C Regulated Genes Cluster1 and cluster2 groups of the m6A/m1A/m5C regulated gene subsets were matched with the tissue samples from 374 HCC patients, and the Limma package was used to determine the differentially expressed genes (DEGs) in cluster1 and cluster2 groups. The false discovery rate (FDR) <0.05 and log fold change (FC) = 2 were used to identify DEGs associated with m6A/m1A/ m5C regulated gene subsets.

Identification of Roles m6A/m1A/m5C Regulated Gene Subgroups Associated With Differentially Expressed Genes
The relationship between the m6A/m1A/m5C regulated gene subgroups related to DEGs and the prognosis of HCC patients was using univariate Cox regression analysis (p < 0.001). Subsequently, the Limma package was used to investigate the expression levels of DEGs that were modulated by m6A/m1A/ m5C in subgroups of genes associated with prognosis in 50 normal liver tissue samples and 374 HCC tissue samples.
Cluster-Related Genes Based on the m6A/ m1A/m5C Regulated Gene Subgroups Using Consensus Clustering The expression data of 37 prognosis related DEGs in 374 HCC tissues were extracted. The optimal k value was determined by consensus clustering analysis, and the 374 HCC tissues from TCGA were grouped. PCA was used to determine the differences between subgroups of HCC patients. The K-M survival analysis and log-rank test were used to determine the relationship between prognostic and clinicopathological factors in patient subgroups.

Construction of a Risk Model and Nomograms Associated With m6A/m1A/ m5C Regulated Genes
The expression data of 45 m6A/m1A/m5C regulated genes were matched and merged with the prognosis data of HCC patients, and the relationship between m6A/m1A/m5C regulated genes and the prognosis of HCC patients was investigated using univariate Cox regression analysis (p < 0.05). The least absolute shrinkage and selection operator (LASSO) regression analysis was used to construct a risk model for m6A/m1A/m5C regulated genes, and the risk score for m6A/m1A/m5C regulated genes was calculated as follows: risk score = S (ExpmRNAn × bmRNAn) (20). Correlation analysis explored the relationship between the risk score and the expression levels of model factors, including YBX1, ZC3H13, YTHDF1, TRMT10C, YTHDF2, RRP8, TRMT6, LRPPRC, and IGF2BP3. The K-M survival analysis was performed to determine the prognosis of HCC patients in the high-and low-risk groups, as well as the construction of nomogram related to m6A/m1A/m5C regulated genes.

Identification of m6A/m1A/m5C Regulated Gene Expression in Hepatocellular Carcinoma Tissues
Cancer and normal liver tissues of 12 HCC patients diagnosed through pathology at our hospital were collected between February and April 2022. All 12 HCC patients signed the informed consent, which was reviewed and approved by the ethics committee of our hospital. With the use of a standard qRT-PCR protocol, the expression levels of YBX1, ZC3H13, YTHDF1, TRMT10C, YTHDF2, RRP8, TRMT6, LRPPRC, and IGF2BP3 in 12 HCC tissues and normal liver tissues were determined (21). The primer numbers for TRMT10C and RRP8 that were purchased from GeneCopoeia (Guangzhou, China) were HQP013870 and HQP097040, respectively. Table 1 shows the primer sequences for the m6A/m1A/m5C regulated genes YBX1, ZC3H13, YTHDF1, YTHDF2, TRMT6, LRPPRC, and IGF2BP3.

Clinical Values of a Risk Model Associated With m6A/m1A/m5C Regulated Genes
The relationship between the high-and low-risk groups and clinicopathological characteristics of HCC patients were assessed using the log-rank test. To understand the clinical value of the risk model, univariate and multivariate Cox regression analyses were performed to establish the relationship between age, gender, grade, clinical stage, T stage, M stage, N stage, and risk score value and OS of HCC patients.

Signaling Mechanisms of the Risk Model
Gene Set Enrichment Analysis (GSEA) is commonly used to investigate the signaling mechanisms of risk model (22). HCC patients were divided into the high-and low-risk groups based on the median risk score. With the use of GSEA (version 4.1.0), the effect of the risk model constructed by m6A/m1A/m5C regulated genes on each gene set was determined. GSEA was run for 1,000 cycles (Nominal p < 0.05).

The Relationship Between the Constructed Model and Immune
Infiltrating Cells Based on the m6A/m1A/ m5C Regulated Genes The expression levels of immune cells in the 374 HCC tissues were calculated using the CIBERSORT method. The immune cells of 374 HCC patients were divided into two groups: cluster1 and cluster2 groups. A t-test was used to determine whether there was a difference in the expression levels of HCC immune cells in cluster1 and cluster2 groups. The high-and low-risk models constructed based on the m6A/m1A/m5C regulated genes were combined with the immune cell data of 374 HCC patients, and the expression levels of immune cells in the highand low-risk groups were explored. In addition, correlation analysis investigated the relationship between the risk score and immune cell infiltration in HCC.
The Relationship Between the Constructed Risk Model of m6A/m1A/m5C Regulated Genes and Immune Cell Markers The expression levels of immune cell markers extracted from 374 HCC tissues in TCGA database were acquired. The risk scores were combined with the immune cell marker gene data of cancer patients, and the expression levels of immune cell markers in the high-and low-risk groups were explored using the t-test. In addition, correlation analysis explored the relationship between the risk scores and immune-infiltrating cell markers in HCC tissues.

Statistical Analysis
The expression levels of m6A/m1A/m5C regulated genes and m6A/m1A/m5C regulated gene subgroup-related genes in HCC tissues were identified using the Limma package. Cox, LASSO, and K-M survival analyses were used to identify m6A/m1A/m5C regulated genes associated with the prognosis of HCC patients. The relationship between models based on m6A/m1A/m5C regulated genes and immune cell infiltration was determined using correlation analysis (p < 0.05).
The Functional Mechanisms and Roles of m6A/m1A/m5C Regulated Genes GO annotation revealed that the m6A/m1A/m5C regulated genes participate in RNA modification and methylation, macromolecule methylation, mRNA methylation and modification, RNA regulation and mRNA stability, ncRNA processing, regulation of mRNA catabolic process, RNA splicing, mRNA catabolic process, mRNA transport, and other processes ( Table S1). KEGG analysis involved the microRNAs in cancer, spliceosome, and RNA transport. Figure S2 showed the PPI network between m6A/m1A/m5C regulated genes.

Evaluation of Clinical Values of m6A/m1A/ m5C Regulated Gene Subgroups in Hepatocellular Carcinoma Patients Using Consensus Clustering
Based on 45 m6A/m1A/m5C regulated gene expression data, and the optimal k value of 2, 374 HCC patients are divided into cluster1 and cluster2 groups (Figures 2A-C). PCA revealed significant differences between cluster1 and cluster2 groups ( Figure 2D). The K-M survival analysis demonstrated that the OS of patients in cluster1 was significantly better than that of patients in cluster2 ( Figure 2E). In addition, there was a statistically significant difference between HCC patients in the two cluster groups in terms of clinical stage, T stage, tumor grade, and fustat status ( Figure 2F).
The Functions and Signaling Pathways of Cluster-Related Differentially Expressed Genes Based on the m6A/m1A/m5C Regulated Genes Compared to cluster1 group, cluster2 group HCC tissues displayed 371 cluster-related DEGs (Table S2), including 315 overexpressed genes and 56 downregulated genes. Figure S3 shows the top 20 cluster-related DEGs in HCC tissues. The m6A/ m1A/m5C regulated gene subsets related DEGs were found to be involved in signaling release, neurotransmitter transport, drug catabolic process, positive regulation of protein secretion, regulation of secretion, calcium ion regulated exocytosis, hormone transport, regulation of protein secretion, and others using GO analysis (

Clinical Values of the Risk Model
Associated With the m6A/m1A/m5C Regulated Genes The risk model was related to the tumor grade, clinical stage, T stage, distant metastasis, and fustat status of HCC patients, according to an analysis of clinicopathological characteristics of HCC patients in the high-and low-risk groups ( Figure 8A). Univariate Cox regression analysis revealed that the clinical stage, T stage, and risk score were risk factors for the dismal prognosis in HCC patients ( Figure 8B). Multivariate Cox regression analysis revealed that distant metastasis and risk score were risk factors for poor prognosis in HCC patients ( Figure 8C).

Signaling Mechanisms Associated With a High-Risk Score
A high-risk score is associated with spliceosome regulation, cell cycle, base excision repair, RNA degradation, oocyte meiosis, ubiquitin-mediated proteolysis, pyrimidine metabolism, DNA  Figure S5 and Table 5).

Immune Cell Infiltration in Hepatocellular Carcinoma Correlates With the Cluster Groups and Risk Model of m6A/m1A/m5C Regulated Genes
There were significant differences in the expression levels of B cells naive, B cells memory, CD4 memory resting T cells, CD4 memory activated T cells, M0 macrophages, M1 macrophages, resting mast cells, and eosinophils among the 45 m6A/m1A/ m5C regulated gene subgroups ( Figure 9A). The expression levels of B cells naive, B cells memory, T cells CD8, CD4 memory resting T cells, M0 macrophages, and eosinophils were significantly different between the 37 cluster-related differentially expressed gene subgroups ( Figure 9B). There were significant differences in the expression levels of the stromal score, CD4 memory resting T cells, follicular helper T cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, resting mast cells, eosinophils, and neutrophils in high-and low-risk groups ( Figure 9C). The risk score was significantly correlated with the levels of the stromal score, CD4 memory resting T cells, CD4 memory activated T cells, follicular helper T cells, M0 macrophages, M1 macrophages, and resting mast cells (Figures 9D-I and Figure S6).

DISCUSSION
Numerous studies indicate that the expression levels of many genes change during HCC progression (23)(24)(25). Wang et al. for example discovered that the expression level of MAGEH1 is significantly downregulated in HCC tissues, which is associated with poor prognosis in HCC patients. MAGEH1 can inhibit proliferation, migration, and invasion and delay HCC progression (24). Yang et al. reported that elevated SHOX2 expression is associated with tumor recurrence and TNM stage in HCC patients. Increased SHOX2 expression is observed in HCC cells. SHOX2 expression can be suppressed to inhibit HCC cell proliferation and invasion (25). Current research demonstrates that RNA modification is associated with the progression of malignant tumors (12-16, 26, 27). For example, the level of IGF2BP3 overexpression correlates with cancer progression and survival. IGF2BP3 knockdown inhibits DNA replication in the S phase of the cell cycle, cell proliferation, and angiogenesis by reading the m6A modification of CCND1 and VEGF (12). WTAP is highly expressed and is a risk factor for poor prognosis in HCC patients. WTAP can promote cell proliferation and tumor growth in HCC patients. ETS1 is a downstream effector of WTAP. WTAP can induce post-transcriptional repression of ETS1 by regulating m6A modification, which in turn mediates  the p21/p27-dependent mechanism involved in the G2/M phase regulation of HCC cells (26). The expression level of METTL3 is significantly elevated in HCC tissues and cells. Elevated METTL3 expression is associated with poor OS. When METTL3 expression is inhibited, the ability of HCC cells to invade, migrate, and proliferate is significantly decreased. Mechanistically, METTL3 may regulate the expression level of USP7 via m6A methylation modification, thereby promoting the growth and migration of HCC cells (27). This suggests that m6A/ m1A/m5C regulated genes play important roles in HCC progression. In this study, the roles of 45 m6A/m1A/m5C regulated genes were explored in HCC progression, and significant differences in OS, clinical stage, T stage, tumor grade, and fustat status between subgroups of m6A/m1A/m5C regulated genes were observed. The identification of the expression of the risk model genes YBX1, ZC3H13, YTHDF1, TRMT10C, YTHDF2, RRP8, TRMT6, LRPPRC, and IGF2BP3 in HCC tissues was consistent with that in TCGA database. In addition, a high-risk score constructed using YBX1, ZC3H13, YTHDF1, TRMT10C, YTHDF2, RRP8, TRMT6, LRPPRC, and IGF2BP3 was associated with poor prognosis, tumor grade, clinical stage, T stage, distant metastasis, and fustat status and is an independent risk factor for poor prognosis in HCC patients. Preliminary evidence suggests that m6A/m1A/m5C regulated genes play important biological roles in the progression of HCC, which can be used to determine the prognosis and survival of HCC patients. The signaling pathways, including the cell cycle, DNA replication, and WNT are closely associated with cancer progression (28)(29)(30)(31)(32)(33)(34). Previous studies have found that ASF1B is highly expressed in HCC, and elevated ASF1B expression level is associated with poor prognosis in HCC patients. GSEA revealed that ASF1B may regulate the cell cycle, DNA replication, and oocyte meiotic signaling pathway. ASF1B silencing inhibits the growth and cell cycle arrest, induces apoptosis, and reduces the expression levels of PCNA, cyclinB1, cyclinE2, and CDK9 in HCC cells (28). The level of TCF3 expression is significantly increased in HCC tissues. TCF3 expression level correlates with clinical stage, tumor size, TNM stage, grade, OS, disease-specific survival, and progression-free survival (PFS) in HCC patients. TCF3 knockdown inhibits cancer cell proliferation and cell cycle, which is associated with dysregulation of the WNT signaling mechanism (31). m6A/m1A/m5C regulated genes are involved in multiple functions, including RNA modification and methylation, mRNA methylation and modification, regulation of RNA, and mRNA stability. The cell cycle, base excision repair, RNA degradation, oocyte meiosis, DNA replication, homologous recombination, basal transcription factors, pathways in cancer, WNT signaling pathway, and long-term potentiation processes are involved in the construction of a high-risk model for m6A/ m1A/m5C regulated genes. Previous studies indicate that m6A/ m1A/m5C regulated genes are closely associated with HCC progression. However, the mechanism of m6A/m1A/m5C regulated genes in HCC progression remains to be confirmed in future studies. Currently, targeted therapy and immunotherapy are among the treatment options available to HCC patients (2,(35)(36)(37)(38). For example, PD-1 inhibitors are well tolerated by HCC patients. HCC patients have achieved good clinical outcomes. The median OS for PD-1 inhibitor-treated HCC patients is 6.6 months, the median PFS is 5.3 months, and the overall response rate is 30.8%. A patient achieves complete remission (38). The relationship between immunotherapy and the immune microenvironment is well established (39,40). In this study, the expression levels of CD4 memory resting T cells, follicular helper T cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, resting mast cells, eosinophils, and neutrophils were found to be significantly different in the risk model constructed by m6A/ m1A/m5C regulated genes YBX1, ZC3H13, YTHDF1, TRMT10C, YTHDF2, RRP8, TRMT6, LRPPRC, and IGF2BP3. The risk score correlated with levels of the stromal score, CD4 memory resting T cells, CD4 memory activated T cells, follicular helper T cells, M0 macrophages, M1 macrophages, and resting mast cells. In addition, the risks core correlated with the HCC immune cell markers, including IL10, ITGAM, STAT5B, CD68,  HLA-DPB1, KIR2DL4, IRF5, CSF1R, CD274, HLA-DRA, CD8B,  STAT1, NOS2, ITGAX, CD86, CD8A, BCL6, TGFB1, CD163,  CCR8, TBX21, CCL2, CD3E, TNF, CD1C, CD2, HAVCR2,  NRP1, STAT5A, CD3D, LAG3, HLA-DPA1, PDCD1, VSIG4,  STAT3, GZMB, MS4A4A, GATA3, IFNG, and HLA-DQB1. This suggests that the risk model constructed using m6A/m1A/m5C regulated genes is associated with HCC immune cell infiltration.
This study identified the involvement of m6A/m5C/m1A regulated genes and constructed a risk model for HCC using HCC data from TCGA database and our hospital data. In the prognosis in HCC patients. The risk score is associated with HCC stromal score and levels of CD4 memory resting T cells, M0 macrophages, M1 macrophages, resting mast cells, CD4 memory activated T cells, and follicular helper T cells.

DATA AVAILABILITY STATEMENT
The data used in database are available from the TCGA website and the data of qRT-PCR by contacting the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the ethics committee of Affiliated Hospital of Zunyi Medical University. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
DK, R-SS, and DL designed the study and supervised the manuscript writing process. DL, KL, WZ, and K-WY analyzed the data and wrote the manuscript. DL, D-AM, and G-JJ performed data visualization and edited the language of the manuscript. All authors approve the submission and publication of the manuscript.