Gene Risk Signature of m6A Regulators, KIAA1429, YTHDF1, YTHDF2, METTL3 and Its Relationship with Immune Inltration in Hepatocellular Carcinoma

Background: Hepatocellular carcinoma (HCC) is a tumor of the digestive system with high mortality. N6-methyladenosine (m6A) is one of the most common forms of mRNA modication. Tumor neoantigens play an important role in anti-tumor immune response and predict the clinical response of immunotherapy. There are few studies on the relationship between m6A regulators and immune inltrating cells. The objective of this study was to determine the gene expression and prognostic value of m6A regulators in hepatocellular carcinoma. Further, we explored the relationship between m6A regulators and immune-inltrating cells. Methods: 13 m6A regulators expressions in 374 cancer tissues and 50 normal tissues were analyzed using RNA-seq data and clinical information in the TCGA database. Survival package was used to analyze the survival of patients in the two groups, and the corresponding clinical characteristics and gene expression were analyzed using univariate and multivariate Cox regression. We evaluated the expression of KIAA1429, YTHDF1, YTHDF2, and METTL3 in HCC and its correlation with TIICs using TIMER and CIBERSORT. Results: Most m6A regulators had higher expression levels in 374 cancer tissues. We conrmed two groups of HCC patients using a consensus clustering method. The prognosis of the cluster 1 group was poor compared with that of the cluster 2 group. The m6A regulators, KIAA1429, YTHDF1, YTHDF2, and METTL3 are highly expressed in the high-risk group of HCC. Furthermore, it was found that four m6A regulators (KIAA1429, YTHDF1, YTHDF2, and METTL3) are closely related to the clinicopathological characteristics and poor prognostic marker of hepatocellular carcinoma. We analyzed the correlation between KIAA1429, YTHDF1, YTHDF2, and METTL3 expression and the level of immune invasion of HCC. Conclusion: KIAA1429, YTHDF1, YTHDF2, and METTL3 as m6A regulators may regulate the tumor microenvironment through tumor immune inltration cells to exert immune anti-tumor effects. KIAA1429, YTHDF1, YTHDF2, and METTL3 as molecular markers providing a new target for therapy of HCC in the further.


Introduction
Hepatocellular carcinoma (HCC) which the most common type of primary liver cancer is a tumor of the digestive system with high mortality. The latest cancer statistics indicate that hepatocellular carcinoma is the sixth most common cancer and the fourth leading cause of death from cancer worldwide. (1).
Unfortunately, most liver cancer patients diagnosed at an advanced stage are not candidates for surgery treatment. Therefore, these patients display a dismal outcome. Novel drugs may improve the outcomes of patients with HCC, especially those in advanced stages. Exploring the molecular mechanism of liver cancer helps identify more potential therapeutic targets, providing more treatment options for HCC.
In this study, the 13 m6A regulators expression in 374 patients with HCC was systematically analyzed using the Cancer Genome Atlas (TCGA) data set. Using bioinformatics and statistical analysis, this study explored the potential application value of m6A regulators and its correlation with immune in ltration in hepatocellular carcinoma.

Datasets
We downloaded the RNA-seq transcriptome data and clinical information of HCC patients using the TCGA database (https:// cancergenome. nih. gov/).

Immune in ltrates analysis
TIMER is a tool for systematic analysis of immune in ltration of different cancer types (https://cistroshinyapps.io/timer/) (31). We evaluated the expression of four m6A regulators in HCC and its correlation with TIICs including B cells, CD 4 + T cells, CD 8 + T cells, neutrophils, dendritic cells, and macrophages. We used CIBERSORT (https://cibersortx.stanford.edu/) to investigate the relationship between high and low levels of the four m6A regulators in the tumor immune microenvironment of HCC.

Bioinformatic Analyses
We analyzed the 13 genes expression in 374 tumor tissues and 50 normal liver tissues using Limma software package. Then, we used the vioplot software package to visualize 13 genes in 374 tumor tissues and 50 normal liver tissues. Next, in this study, 50 normal tissue samples were extracted, 374 cancer tissues were grouped using the Consensus ClusterPlus package, and the grouping results were veri ed using Principal Component Analysis (PCA). Finally, we analyzed the survival of the two clusters using the survival package. We used univariate and multivariate Cox regression to analyze the corresponding clinical features and gene expression in the TCGA data set.

Statistical Analyses
The 13 m6A regulators expression in 374 HCC tissues and 50 normal liver tissues in the TCGA dataset were contrasted using one-way ANOVA. The clinical characteristics of the primary liver cancer samples in the TCGA dataset were compared using a t-test. Kaplan-Meier method was used for the bilateral logarithmic rank test in overall survival analysis. All the statistical analyses in this study were implemented by R v3.6.3, and P < 0.05 means the difference is statistically signi cant.

Results
Identi cation of m6A regulators in hepatocellular carcinoma.

Consistent clustering of m6A regulators determined two clusters
In this study, we used the Consensus Cluster Plus package to extract 50 normal liver tissues and 374 liver cancer tissues. On the basis of the expression similarity of m6A regulators in the TCGA dataset, k = 2 seemed to have the smaller Cumulative distribution function (CDF) value (Figure 2 B, C). We divided the patients into two groups ( Figure 2A). Next, we used PCA to analyze the two clusters. The results showed that cluster 1 and cluster 2 were distinct. ( Figure 2D) Correlation between consensus cluster classi cation and clinical characteristics We analyzed the OS curves of two groups of patients. Kaplan-Meier survival curve was used to analyze OS rate of HCC patients between two clusters. In this study, we found that the OS of cluster 1 was striking shorter than cluster 2 ( Figure 3A). Most m6A regulators were highly expressed in the cluster 2. Patients of HCC in the cluster 1 had higher scores compared with the cluster 2. In this study, we investigated the positive relationship of m6A regulator and the pathological characteristics such as age, gender, grade,  Figure 4B, C). We separated the patients into a low-risk group and a high-risk group. The results indicated that the high-risk group of patients had poor survival ( Figure 4D).

Correlation between prognostic risk score and clinicopathological characteristics
We investigated the role of the ve genes in the high and low-risk patients, and the pathological characteristics including age, stage, TMN status of liver cancer ( Figure 5A). Higher expression of KIAA1429, YTHDF1, YTHDF2, and METTL3 and lower expression of ZC3H13 in hepatocellular carcinoma ( Figure 5A). Using the ROC curve, we predicted prognostic risk scores and 3-year survival rates in HCC patients. The consequences showed that risk score predicted 3-year survival in patients (AUC = 0.619) ( Figure 5B). Univariate and multivariate Cox regression analysis con rmed that the risk score, stage, T status were all related to OS ( Figure 5C,D). The clinicopathological characteristics of primary liver cancer, and KIAA1429, YTHDF1, YTHDF2, and METTL3 were risk factors for poor prognosis of liver cancer patients.
Relationship between KIAA1429, YTHDF1, YTHDF2, and METTL3 expression and tumor-in ltrating immune cells Independent tumor-in ltrating lymphocytes play a pivotal role in predicting prognosis (32). This study used TIMER to analyze the correlation between KIAA1429, YTHDF1, YTHDF2, and METTL3 expression and the level of immune invasion of liver cancer.  (Figure 6A-6D). We further investigated whether there is a difference in the tumor immune microenvironment between high and low levels of the four m6A regulators in HCC. The CIBERSORT algorithm was applied to 22 immune cell subtypes to evaluate the difference in their expression levels between the high expression group and the low expression group ( Figure 7A-7D). There was no signi cant difference in immune in ltrating cells between the KIAA1429 high expression group and the low expression group. Compared with the low expression group, the neutrophils in the YTHDF1 high expression group increased (p-value=0.039), the neutrophils in the YTHDF2 high expression group increased (p-value = 0.023), and the METTL3 memory B cells increased (p-value=0.006). The results indicate that KIAA1429, YTHDF1, YTHDF2, and METTL3 may play a crucial part in immune invasion of liver cancer.

Discussion
Globally, liver cancer is the sixth most frequent cancer and the fourth leading cause of cancer death. (1) Surgery is the preferred treatment for patients with liver cancer at an early stage (33). More evidence indicated that primary liver cancer is a multistep process including the complex interactions between genetics, epigenetics, and transcriptional changes (34).
The m6A modi cation usually rst causes the m6A methylation group to be written into the RNA through coding genes called writers (35) and then eraser genes remove the m6A methylation group on the RNA, thereby affecting the tumor biological process (30). Finally, the reader gene binds to the m6A site on the RNA in the nucleus or cytoplasm to play a speci c biological role (36,37). The expression level of m6A gene and protein may a new target of molecular targeted therapy for primary liver cancer.
In this study, we investigated that there is a close relationship between the m6A regulators and pathological characteristics of prime liver cancer. KIAA1429, YTHDF1, YTHDF2, and METTL3 were higher expressed in a high-risk group of liver cancer. The risk score is an independent prognostic indicator, which is related to the clinicopathological features of primary liver cancer.
As a downstream target of METTL3, cytokine signaling inhibitor 2 (SOCS2) is also a tumor suppressor.
Through the m6A-YTHDF2-dependent signal way, METTL3 reduced the stability of SOCS2 and promotes the occurrence and development of HCC (38). METTL3, YTHDF2, and ZC3H13 were found to be independent prognostic factors. In HCC, the expression of METTL3 is up-regulated through CNV and DNA methylation mechanisms. METTL3 may be a biomarker for the prognosis of liver cancer (39). YTHDF1 expression is related to the pathological stage of patients with HCC. Patients with low expression of YTHDF1 related to higher survival rates (40). The decreased expression of METTL14 enhanced the ability of metastasis in liver cancer cells, and the OS rate of liver cancer patients was lower. Increased METTL14 expression inhibited the ability of migration and invasion in vitro tumor cell (41). METTL14 may play a role in the progression of HCC by regulating the m6A levels of CSAD, GOT2, and SOCS2 (42). In this study, we also found that KIAA1429 protein expression was not detected in normal liver tissues and HCC tissues. Protein expression data of YTHDF1 and METTL3 were missing from the HPA data. The expression of YTHDF2 was not detected in normal tissues or moderately expressed in liver tissues, and the expression of the four genes in normal tissues and liver cancer tissues was not consistent with their mRNA expression levels. This conclusion is only based on the HPA database and needs to be con rmed by further experimental data. The study demonstrated the expression and poor prognosis of m6A regulators in HCC. The expression of m6A regulator are positively correlated with the clinicopathological features of HCC.
The composition of immune cells in the tumor microenvironment is closely related to prognosis.M2 macrophages and pro-in ammatory M1 macrophages are associated with the growth and spread of tumors and the formation of immunosuppressive microenvironments (43,44). High concentrations of CD 3+ T cells, CD 8+ cytotoxic T cells, and CD 45RO+ memory T cells suggest a better prognosis. (45,46) Compared with normal livers, there were fewer mast cells, monocytes, and plasma cells activated in HCC, and more resting mast cells, total B cells, and naive B cells, CD 4+ memory resting T cells, and CD 8+ T cells (47). Tumor neoantigens play an important role in anti-tumor immune response and predict the clinical response of immunotherapy (48,49). The potential role of m6A modi cation in the host's anti-tumor immune response is unclear. Studies have shown that long-lasting new antigen-speci c immunity can be controlled by m6A modi cation of m6A binding protein YTHDF1 to achieve anti-tumor therapy, suggesting that YTHDF1 may be a potential therapeutic target for anti-cancer immunotherapy (50). We found that KIAA1429, YTHDF1, YTHDF2, and METTL3 are positively correlated with various immune cells, suggesting that they may participate in the tumor microenvironment through the composition of immune cells and play a role in the occurrence and development of tumors.

Conclusions
KIAA1429, YTHDF1, YTHDF2, and METTL3 can not only be used for tumors by m6A modi cation mechanism, but also may regulate the tumor microenvironment through tumor immune in ltration cells to exert immune anti-tumor effects. KIAA1429, YTHDF1, YTHDF2, and METTL3 as molecular markers providing a new target for therapy of HCC in the further.