Development of N6-methyladenosine Related Signature as New Biomarker for Prognosis of Hepatocellular Carcinoma and Correlates with Sorafenib and anti-PD-1 Immunotherapy Treatment Response

Background: N6-methyladenosine (m6A) modication plays an essential role in diverse key biological processes and may take part in the development and progression of hepatocellular carcinoma (HCC). Here, we systematically analyzed expression proles and prognostic values of 13 widely reported m6A modication related genes in HCC. Methods: mRNA expression of 13 m6A modication related genes and clinical parameters of HCC patients were downloaded from TCGA, ICGC, GSE109211 and GSE78220. Univariate and LASSO analysis was used to develop risk signature. Time-dependent ROC was performed to assess the predictive accuracy and sensitivity of risk signature. Results: FTO, YTHDC1, YTHDC2, ALKBH5, KIAA1429, HNRNPC, METTL3, RBM15, YTHDF2, YTHDF1 and WTAP were signicantly overexpressed in HCC patients. YTHDF1, HNRNPC, RBM15, METTL3, YTHDF2 were independent prognostic factors for OS and DFS in HCC patients. Next, a risk signature was also developed and validated with ve m6A modication related genes in TCGA and ICGC HCC cohort. It could effectively stratify HCC patients into high risk patients with shorter OS and DFS and low risk patients with longer OS and DFS and showed good predictive eciency in predicting OS and DFS. Moreover, signicantly higher proportions of macrophages M0 cells, neutrophils and Tregs were found to be enriched in HCC patients with high risk score, while signicantly higher proportions of memory CD4 T cells, gamma delta T cells and naive B cells were found to be enriched in HCC patients with high low score. Finally, signicantly lower risk scores were found at sorafenib treatment responders and anti-PD-1 immunotherapy responders compared to that in non-responders, and anti-PD-1 immunotherapy treated patients with lower risk score had better OS than patients with higher risk score. Conclusion: A risk signature developed with the expression of 5 m6A related genes could improve the prediction of prognosis of HCC and correlate with sorafenib treatment and anti-PD-1 immunotherapy response. we Univariate analysis showed that higher expression of YTHDF1, WTAP, HNRNPC, RBM15, METTL3, KIAA1429, YTHDC1, YTHDF2 and lower expression of ZC3H13 were statistically related to poorer OS of HCC patients (all p<0.05, supplementary gure 1A); multivariate analysis showed that expression of YTHDF1, WTAP, HNRNPC, RBM15, METTL3, KIAA1429 and YTHDF2 still remained signicantly related with OS after adjusting for gender, age, histologic grade, T stage, N stage, M stage and TNM stage (all p<0.05, supplementary gure 1B-1J).Then, the prognostic values of m6A modication related genes for recurrence of HCC patients were also analyzed. Univariate analysis indicated that overexpression of YTHDF1, WTAP, HNRNPC, RBM15, METTL3, YTHDC1 and YTHDF2 were statistically related with shorter DFS (all p<0.05, supplementary gure 2A); and, multivariate analysis showed that expression of YTHDF1, HNRNPC, RBM15, METTL3 and YTHDF2 were still statistically related with DFS after adjusting for gender, age, histologic grade, T stage, N stage, M stage and TNM stage (all p<0.05, supplementary gure 2B-2H). These results strongly conrmed the important roles played by m6A modication related genes in HCC. multivariate prognostic T stage and risk of patients (all p<0.05, gure 4A). risk signature still remained statistically related with OS after adjusting for T stage, M stage and TNM stage by multivariate analysis. In multivariate analysis, after adjusting for TNM stage, the risk signature was still signicantly related with OS (p<0.01, gure 4B). univariate analysis T stage, TNM stage and risk signature DFS patients. In univariate analysis, T stage, TNM stage and risk signature signicantly associated DFS in patients (all p<0.001, gure 4C). these into multivariate analysis, the only the risk signature was statistically related with DFS (p<0.001, gure 4D). To conclude, these results indicated that risk signature was independent prognostic factor for OS and DFS of HCC patients.


Introduction
Hepatocellular carcinoma (HCC) is a common type of cancer and represents the leading cause of cancerrelated death worldwide. HCC is still a serious burden to public health (1). There were about 841,000 patients developed HCC and 782,000 patients died from HCC alone in 2018 because of late diagnosis and limited treatment options (1,2). Moreover, incidence of HCC is increasing rapidly with 50% recurrence rate after surgical treatment (3,4). It is well recognized that development and progression of HCC is the result of multistep process, where interactions between genetics and epigenetics have played important roles (5)(6)(7)(8). By understanding the pathogenesis of HCC is key to discover new diagnostic biomarkers and therapeutic targets.
RNA modi cation, discovered in the 1970s, has recently been recognized as a third layer of epigenetics that could modify a plethora of native cellular RNAs. (9)(10)(11). N6-methyladenosine (m6A) modi cation is the most abundant form of internal mRNA methylation among the kinds of RNA modi cations in eukaryotes (12). m6A modi cations in mammalian cells are dynamic and reversible, and are commonly regulated by binding proteins ('readers'), methyl-transferases ('writers'), demethylases ('erasers') (13). Among m6A modi cation related genes, 13 genes, including ZC3H13, WTAP, KIAA1429, METTL3, METTL14, RBM15, YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC, ALKBH5 and FTO, are the most prominent m6A modi cation related genes (14)(15)(16). These m6A modi cation related genes are primarily involved in modulation of alternative mRNA splicing, precession of pre-miRNA, stability of mRNA and enhance of translation e ciency of mRNA [13]. Not only do these 13 m6A modi cation related genes play essential roles in many important biological processes, such as development of embryonic and neural cell, differentiation of stem cell, and stress responses (17)(18)(19), they also take park in the development, progression and radio resistance of various kinds of cancers (20)(21)(22)(23). For example, overexpression of YTHDF1 is found to be related with poorer survival of HCC patients, and KIAA1429, METTL3 are found to regulate migration and invasion of HCC, indicating an important roles of m6A modi cations related genes played in HCC (24)(25)(26).
Recently, Zhou et al explored expression pattern and prognostic values of m6A modi cation related genes of HCC patients, but they mainly focused on the role of METTL3and YTHDF1 (27). In the present study, we comprehensively analyzed the expression pattern and prognosis of thirteen widely reported m6A modi cation related genes in TCGA HCC cohort. Besides, we also developed and validated a risk signature with expression of 5 selected m6A modi cation related genes, and analyzed its prognostic value for HCC patients and its relation with tumor-in ltrating immune cells in TCGA and ICGC HCC cohort.
Moreover, the prediction values of risk signature in sorafenib treatment and anti-PD-1 immunotherapy response was also evaluated.

Ethics statement
All the data analyzed in the present study were received from TCGA, ICGC and GEO dataset, written consents were already obtained before our study.
Data collection mRNA expression of TCGA HCC cohorts, which included 374 HCC cases and 50 normal controls, were got from GDC Data portal (https://cancergenome.nih.gov/). Meanwhile, corresponding clinical-pathological data, including gender, age, histologic grade, tumor T stage, N stage, M stage (M), TNM stage, overall survival (OS) time and disease-free survival (DFS) time were also downloaded. It was of note that 9 of 374 HCC patients were excluded because of absence of corresponding clinical-pathological data, and basic characteristics of 365 HCC patients were summarized in table 1. In addition, a total of 232 HCC patients with available OS information and mRNA expression were got from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). mRNA expression of 67 sora nib-treated HCC patients of GSE109211was downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/) and there were 21 sora nib treatment responders and 46 non-responder in GSE109211. Moreover, mRNA expression of 27 melanomas patients with anti-PD-1 checkpoint inhibition therapy of GSE78220 was also downloaded from GEO database. 4 patients achieved complete response, 10 patients achieved partial response and 13 patients achieved no response.

Data analysis ow chart
To make the study to be better understood, a work ow of the study was depicted and was shown at gure 1.

Expression of m6A modi cation related genes of HCC patients and their associations with clinical-pathologic parameters
First, mRNA expression of 13 m6A modi cation related genes were download from TCGA and compared between HCC patients and normal controls. As was shown at gure 2A and 2B, signi cantly higher expression of FTO, YTHDC1, YTHDC2, ALKBH5, KIAA1429, HNRNPC, METTL3, RBM15, YTHDF2, YTHDF1 and WTAP were found at HCC compared to normal tissues (all p<0.001). Interestingly, we also found that expression of most of the 13 m6A modi cation related genes seemed to be lower than that of other 32 kinds of tumors. Besides, most of the 13 m6A modi cation related genes were positively correlated with each other ( gure 2C). Moreover, genetic change, such as missense mutation, truncating mutation, ampli cation, deep deletion, diploid, gain were observed in near 80% of the HCC patients ( gure 2D), Speci cally, each HCC patient may have one or more kinds of genetic changes. The genetic rate of WTAP, KIAA1429, RBM15, METTL3, METTL14, ALKBH5, YTHDC1, YTHDC2, HNRNPC, YTHDF1, YTHDF2, FTO and ZC3H13 were 7%, 4%, 17%, 40%, 5%, 5%, 7%, 8%, 18%, 11%, 9%, 13%, 17%, respectively, suggesting that higher expression of m6A modi cation related genes may be the result of genetic changes in related genes. Taken together, these results indicated that m6A modi cation related genes played important roles in HCC.
Prognostic value of m6A modi cation related genes in HCC cases These results strongly con rmed the important roles played by m6A modi cation related genes in HCC.
Development of risk signature with 5 m6A modi cation related genes and its association with clinical-pathologic parameters To better explore the prognostic value of m6A modi cation related genes, a risk signature was developed. Based on the results of univariate analysis ( gure 3A), ZC3H13, YTHDF1, WTAP, HNRNPC, RBM15, METTL3, KIAA1429, YTHDC1 and YTHDF2 were associated OS and were considered as prognosticrelated genes. Then, LASSO analysis was used to further screen the prognostic-related genes. In the end, 5 genes, including YTHDF2, YTHDF1, METTL3, KIAA1429 and ZC3H13, were used to develop the risk signature ( gure 3A, 3B). The risk score was then constructed based on the coe cients weighted by LASSO analysis and calculated as follows: risk score= (0.07*YTHDF2) + (0.02*YTHDF1) + (0.11*METTL3) + (0.04*KIAA1429) -(0.1*ZC3H13). We calculated risk score for every HCC cases and assigned them into high risk group and low risk group on the basis of the median risk score. Expression of YTHDF2, YTHDF1, METTL3, KIAA1429 tend to be higher in patients with high-risk score; and expression of ZC3H13 seemed to be higher in patients with low-risk scores ( gure 3C Prognostic value of risk signature for OS and DFS of HCC cases The risk signature was found to be associated with clinical-pathologic parameters. We next performed univariate and multivariate analysis to analyze its prognostic value. Based on the univariate analysis, T stage, M stage, TNM stage and risk signature were statistically related with OS of HCC patients (all p<0.05, gure 4A). The risk signature still remained statistically related with OS after adjusting for T stage, M stage and TNM stage by multivariate analysis. In multivariate analysis, after adjusting for TNM stage, the risk signature was still signi cantly related with OS (p<0.01, gure 4B). Similarly, univariate analysis also showed that T stage, TNM stage and risk signature were statistically related with DFS of HCC patients. In univariate analysis, T stage, TNM stage and the risk signature were also signi cantly associated with DFS in HCC patients (all p<0.001, gure 4C). By incorporating these factors into multivariate analysis, the result suggested only the risk signature was statistically related with DFS (p<0.001, gure 4D). To conclude, these results indicated that risk signature was independent prognostic factor for OS and DFS of HCC patients.

Validation of risk signature
To independently test the applicability of the signature, 232 HCC patients with available OS information from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP) was further used to examine the applicability of the signature. Risk score for every patient was computed. Similarly, the signature could effectively stratify high risk HCC patients with poorer OS and low risk patients with better OS (HR=2.309, 95%CI: 1.302-4.369, P=0.006, gure 6A). Moreover, the AUC of risk signature for predicting 1-, 3-and 5-year OS were 0.7, 0.74 and 0.714 ( gure 6B), respectively, which convincingly suggested the good discrimination and prediction of our signature.
Correlation of risk signature with tumor-in ltrating immune cells in TCGA and ICGC HCC cohort CIBERSOR was used to calculate 22 kinds of in ltrating immune cells in patients with different risk scores. In TCGA HCC cohort, signi cantly higher proportions of macrophages M0 cells, memory B cells, follicular helper T cells and neutrophils were found to be enriched in HCC patients with high risk score, while signi cantly higher proportions of resting memory CD4 T cells and monocytes were found to be enriched in HCC patients with low risk score (all p<0.05, gure 7A). In ICGC HCC cohort, signi cantly higher proportions of macrophages M0 cells and Treg cells were found to be enriched in HCC patients with high risk score, while signi cantly higher proportions of naive B cells and gamma delta T cells were found to be enriched in HCC patients with low risk score (all p<0.05, gure 7B). These results suggested that the risk signature was signi cantly associated with tumor-in ltrating immune cells, and different kinds of in ltrating immune cells in patients with different risk score may contribute to their different prognosis.
Risk signature as indicator in sora nib treatment response for HCC patients To investigate the association between risk signature and sorafenib treatment response, we calculated risk score for each HCC patients treated with sorafenib of GSE109211, which contained 21 sorafenib treatment responders and 46 non-responders. Signi cantly lower risk scores were found at sorafenib treatment responders compared to that in non-responders (P < 0.001, gure 8A). Moreover, the AUC for predicting sora nib treatment response was 0.794 ( gure 8B). Taken together, the risk signature may be served as an indicator for sora nib treatment response in HCC patients.
Correlation of risk signature with anti-PD-1 immunotherapy As a major breakthrough in cancer therapy, immunotherapies represented by immunological checkpoint blockade (PD-1/L1 and CTLA-4) has proved promising clinical e cacy and previous study have proved that combination treatment with anti-PD-1 antibodies and Sorafenib exhibited a more potent antitumor effect, but only a small number of patients could achieved durable responses (31,32), so in the present study, we also explored whether the risk signature could predict patients' response to immune checkpoint blockade therapy in an anti-PD-1 cohort of GSE78220. Encouragingly, patients with lower risk score had better OS than patients with higher risk score (HR=3.81, 95%CI: 1.13-11.08, p=0.03, gure 9A). Besides, despite there was no statistical difference, lower risk score was found at patients with complete immunotherapeutic response compared to that in patients with partial response and patients with no response, and lower risk score was also found in alive patients treated with anti-PD-1 than that in patients of death, which may due to the limitation number of patients in the cohort ( gure 9B-9C). Moreover, the AUC of the risk signature for predicting 1 year-, 1.5-year and 2-year OS of patients with anti-PD-1 immunotherapies was 0.669, 0.725 and 0.639 ( gure 9D). In a word, the aboved results strongly indicated that risk signature was signi cantly correlated with response to anti-PD-1 immunotherapy, which may be used as a new biomarker for predicting the response to anti-PD-1/L1 immunotherapy. Discussion m6A modi cations are mainly controlled by methyl-transferases, binding proteins and (13). Studies has reported conservative role and mechanism of m6A modi cation related genes in regulating RNA modi cation, but only a few literatures have studied the role of m6A modi cation related genes in HCC patients. Zhou et al found that YTHDF1 was signi cantly up-regulated in HCC and positively correlated with pathology stage (24). Cheng et al also reported that expression of KIAA1429 was higher in HCC and HCC cell lines, and KIAA1429 could regulate the progression of HCC by regulating ID2 m6A modi cation (26). Chen et al discovered that METTL3 was signi cantly up-regulated in HCC. Knockdown of METTL3 was also found to suppress tumorigenicity and progression of HCC through YTHDF2-dependent post transcriptional silencing of SOCS2 (25). Moreover, Yang et al found YTHDF2 was signi cantly related to malignancy of HCC, and miR-145 could inhibit tumorigenicity of HCC by decreasing YTHDF2 (33). Collectively, these results indicated that m6A modi cation related genes promoted the tumorigenesis of HCC.
Whether expressions of m6A modi cation related genes could be considered as prognostic biomarker is one of the trending research topics in m6A modi cation research (20). Up-regulation of YTHDF1 and METTL3 expression were found to be related to poorer OS of HCC patients (24,25,27). Similarly, in our study, THDF1, HNRNPC, RBM15, METTL3, YTHDF2 are independent prognostic factors for OS and DFS in HCC patients. Next, a risk signature based on the expression ve genes could differentiate HCC patients into high risk patients with poorer OS and DFS and low risk patients with better OS and DFS. Interestingly, this risk signature together showed better predictive e ciency in predicting OS and DFS than TNM stage, or any single gene estimation alone. Therefore, this risk signature may be an advantageous method for individualized therapeutic strategies in HCC patients. In addition, we also found that the risk signature was signi cantly associated with tumor-in ltrating immune cells, which may in uence prognosis of patients with different risk scores. Signi cantly higher proportions of macrophages M0 cells, neutrophils and Tregs cells were found to be enriched in HCC patients with high risk score. Previous studies have showed that macrophages could be recruited to tumor tissues and become pro-angiogenic cells, which was signi cantly associated with micro-vessel density and poor OS and DFS of HCC (34,35); Zhou et al also found that tumor-associated neutrophils could promote the progression of HCC and resistance to sorafenib by recruiting macrophages and Tregs cells (36). These results may partly explain the reason for poorer OS and DFS in HCC patients with high risk score. Moreover, signi cantly higher proportions of memory CD4 T cells, gamma delta T cells and naive B cells were found to be enriched in HCC patients with low risk score, suggesting higher proportions in ltrated T cells and B cells. Garnelo et al. found that the degree of in ltrated T cells and B cells of tumor tissues signi cantly related with improved prognosis of HCC patients (37), which may also partly explain the reason for longer OS and DFS in HCC patients with low risk score.
As an oral multi kinase inhibitor, sorafenib is one of the standard care therapy for advanced stage HCC patients approved by FDA. It can prolong the survival time of HCC patients by inhibiting cell proliferation and angiogenesis, and promoting cell apoptosis through inhibiting a variety of intracellular and cell surface kinases (such as c-raf, BRAF, RET) vascular endothelial growth factor receptor (VEGFR) and platelet-derived growth factor receptor (PDGFR) (38,39). However, some studies have also found that HCC rapidly become sorafenib-resistant and only about 30% of patients can bene t from sorafenib treatment, which may greatly limit the wide clinical application of sorafenib (40,41). Besides, as a major breakthrough in cancer therapy, immunotherapies represented by immunological checkpoint blockade (PD-1/L1 and CTLA-4) has proved promising clinical e cacy and previous study have proved that combination treatment with anti-PD-1 antibodies and Sorafenib exhibited a more potent antitumor effect, but only a small number of patients could achieved durable responses (31,32), so identifying the HCC patients suitable for sorafenib treatment or anti-PD-1 immunotherapy or their combination therapy may be urgent and clinical signi cance. Encouragingly, in the present study, we found the m6A-related risk signature was signi cantly correlated with response to sorafenib treatment and anti-PD-1 immunotherapy. Signi cantly lower risk scores were found at sorafenib treatment responders or anti-PD-1 immunotherapy responders and anti-PD-1 immunotherapy treated patients with lower risk score had better OS than patients with higher risk score, which strongly indicated that the risk signature may be used as a new biomarker for predicting the response to sorafenib treatment and anti-PD-1immunotherapy, and even the combination of them, but independent prospective studies with a larger sample size are still needed to con rm our ndings.
Though the risk signature exhibited good performance for the prognosis of HCC, several limitations should be addressed. First of all, although the prognostic value of the risk signature has been validated in external cohort, independent cohorts consist of more HCC patients are required to further verify the model. Secondly, we did not explore the potential biological functions and pathways of risk signature. In vitro and in vivo experiment should be carried out to uncover the relevant mechanisms. Finally, previously, Huang et al. suggested signi cant expression of m6A modi cation related genes were found in circulating tumor cells (CTCs) (42). Further studies are needed to examine whether these m6A modi cation related genes could be detected in peripheral blood in HCC patients and whether the risk signature in blood could still have good prognostic value.
In conclusion, THDF1, HNRNPC, RBM15, METTL3, YTHDF2 are independent prognostic factors for OS and DFS in HCC patients. A risk signature developed with the expression of YTHDF2, YTHDF1, METTL3, KIAA1429 and ZC3H1 could improve the prediction of prognosis and correlate with sorafenib treatment and anti-PD-1 immunotherapy response.

Declarations
Ethics approval Figure 1 The work ow chart of the present study.        Association of risk signature with anti PD-1 immunotherapy treatment response of GSE78220 cohort. Kaplan-Meier analysis of OS of anti PD-1 immunotherapy treated patients with different risk score (A); Difference of risk score among complete anti PD-1 immunotherapy response, partial anti PD-1 immunotherapy response and no anti PD-1 immunotherapy response (B); Difference of risk score between alive patient with anti PD-1 immunotherapy and dead patients with anti PD-1 immunotherapy (C); AUC of risk signature in predicting1-year, 1.5-year and 2-year OS in patients 1-year, 3-year and 5-year OS in patients with anti PD-1 immunotherapy response (D).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.