Identification and Characterization of Alcohol-related Hepatocellular Carcinoma Prognostic Subtypes based on an Integrative N6-methyladenosine methylation Model

Background: Alcohol consumption increases the risk of hepatocellular carcinoma (HCC), and associated with a high mortality rate and poor prognosis. N6-methyladenosine (m6A) methylations play key roles in tumorigenesis and progression. However, our current knowledge about m6A in alcohol-related HCC (A-HCC) remains elucidated. Herein, the authors construct an integrative m6A model based on A-HCC subtyping and mechanism exploration workflow. Methods: Based on the m6A expressions of A-HCC and in vivo experiment, different prognosis risk A-HCC subtypes are identified. Meanwhile, multiple interdependent indicators of prognosis including patient survival rate, clinical pathological prognosis and immunotherapy sensitivity. Results: The m6A model includes LRPPRC, YTHDF2, KIAA14219, and RBM15B, classified A-HCC patients into high/low-risk subtypes. The high-risk subtype compared to the low-risk subtype showed phenotypic malignancy, poor prognosis, immunosuppression, and activation of tumorigenesis and proliferation-related pathways, including the E2F target, DNA repair, and mTORC1 signalling pathways. The expression of Immunosuppressive cytokines DNMT1/EZH2 was up-regulated in A-HCC patients, and teniposide may be a potential therapeutic drug for A-HCC. Conclusion: Our model redefined A-HCC prognosis risk, identified potential m6As linking tumour progress and immune regulations and selected possible therapy target, thus promoting understanding and clinical applications about A-HCC.

Abnormal activity of m6A-related enzymes can induce numerous pathological responses, including tumour development, immune dysregulation, and embryonic developmental delay [4][5][6]. Therefore, understanding the m6A alterations can provide meaningful insights into the tumorigenesis, progression and potential therapeutic targets.

Ivyspring International Publisher
Liver cancer is one of the most common cancers and the fourth leading cause of cancer-related deaths, with more than 850,000 annually reported cases worldwide [7][8][9][10]. The most common risk factors for HCC are viral infection, alcoholism, and metabolic syndrome [11][12][13]. Alcohol can significantly increase the release of cytoplasmic mtDNA, thereby activating the cGAS-IRF3 signal, causing liver cell apoptosis and inflammation [14,15]. Excessive alcohol intake induces a metabolism shift from oxidative metabolism to reductive reaction, which favour the synthesis of fatty acids and promote fatty liver progress, leading to liver cirrhosis and cancer [16,17]. In addition, chronic alcohol intake can cause reactive oxygen species and DNA damage, and further promote the activation of cancer stem cell-related gene mutations, leading to a poor prognosis for A-HCC [18], which has a mortality rate that is 4 times that of the general population [19].
The specific molecular mechanisms underlying A-HCC remain to be elucidated. The two most recognised major drivers are cytochrome P450 2E1 (CYP2E1) and intestinal lipopolysaccharide (LPS) imbalance [20,21]. Alcohol could induce liver inflammation and oxidative stress cause DNA damage in hepatocytes; ultimately promote tumour initiation and progression [22]. Previously, m6A methylation was reported to play a promoting role in the occurrence and development of HCC, regulating cell proliferation, cell invasion and epithelial to mesenchymal transformation [23]. The levels and activities of m6A regulatory genes YTHDF2, ALKBH5 and FTO can inhibit the HCC malignancy [24][25][26]. For example, FTO can control liver energy homeostasis and metabolism, and it plays an anticancer role in the HCC development [27].
Here, to further explore the correlation between the level of m6A methylations and the occurrence and prognosis of A-HCC. We propose an integrative m6A model based on A-HCC subtyping and mechanism exploration workflow. Then, based on the m6A regulatory factors and multi-omics data from the cancer genome atlas (TCGA) two A-HCC subtypes and their corresponding biological and clinical characters were identified. We observed high-risk A-HCC subtypes are related to immunosuppression and some key Immunosuppressive cytokines (EZH2 and DNMT1) promote the poor prognosis of A-HCC patients. In addition, we selected possible therapy target, thereby promoting a comprehensive understanding of A-HCC and providing guidelines for its treatment.

Patients and specimens
For this study, we collected samples from 108 patients who underwent a liver biopsy at Zhujiang Hospital (Southern Medical University, Guangzhou, China) between 2018 and February 2021. After formalin fixation for 24 h, the samples were dehydrated, embedded in paraffin, and stored at 4 °C. The samples were divided into three groups: normal (no HCC history, n = 31), N-A-HCC (no history of alcohol consumption, n = 56), and A-HCC (history of alcohol consumption for more than 20 years, n = 21) (Supplementary Table 1). The data and tissue samples used in this study met the medical ethical requirements of the Southern Medical University.

Mice, diets, and experimental design
C57BL/6 mice were obtained from the Guangdong Animal Experiment Center, China and they were kept in a specific pathogen-free environment at a constant temperature and light-dark cycle of 12 h. All animal handling procedures were approved by the Southern Medical University Animal Care and Use Committee.
To establish a tumour model, C57BL/6 mice were intraperitoneal injected with 25 mg/kg diethylnitrosamine (DEN; Sigma, USA) at 2 weeks of age. At 6 months, mice injected with DEN were given a non-alcoholic liquid diet for diet adaption. One week later, the experimental group was switched to an alcoholic liquid diet (the alcohol concentration was gradually increased to 4.8%), while the control group continued to receive a non-alcoholic liquid diet (maltose instead of alcohol with the same caloric content). Mice in the experimental treatment group were given teniposide (0.4 mg/d per kg bodyweight; Topscience, China) via gavage for 8 weeks at 6 months of age.

Cell culture and reagents
The human HCC-derived cell lines Huh7 and HepG2 were obtained from the American Type Culture Collection (ATCC, Manassas, VA, USA). Cells were cultured in Dulbecco's modified Eagle's medium (DMEM; Gibco, USA) supplemented with 10% foetal bovine serum (FBS; Gibco) and 1% penicillin/streptomycin (Gibco) at 37 °C in a 5% CO2 atmosphere. Mid-log phase cells were used in all experiments.
When the cells achieved the desired confluency, the cells were starved by culturing in medium without FBS for 24 h. Subsequently, the cells were incubated with medium containing 100 mM ethanol for 48 h at 37 °C in a 5% CO2 humidified environment.
These cells were used as in vitro model of A-HCC. Then, the cells were incubated for 12 h with 0.5 μM teniposide (Alexis Biochemicals, San Diego, CA, USA).

Quantitative reverse transcription PCR (qRT-PCR)
Total RNA was extracted from clinical patient samples and cells using TRIzol reagent (Life Technologies, USA). RNA concentration and quality were measured using a spectrophotometer (Nanodrop One, Thermo Fisher Scientific, Waltham, MA, USA). Samples with a 260/280 absorbance ratio > 2 ± 0.1 were considered contaminated with protein and discarded. RNA samples were then reverse transcribed to cDNA using a reverse transcription kit (#RR037A; Takara Bio, Shiga, Japan), after which qPCR was performed using SYBR Premix Ex Taq (DRR041A; Takara Bio). Specific primers as shown in Supplementary Table 2 were used to detect the expression levels of relevant genes.
For immunofluorescence, the cells were fixed with 4% paraformaldehyde, incubated with Triton, blocked with goat serum, and incubated with primary antibodies against DNMT1 and EZH2 at 4 °C overnight and with secondary antibodies (ab150077; Abcam) at room temperature for 1 h. The nuclei were counterstained with DAPI, after which the samples were imaged using a fluorescence microscope.

Data acquisition
Tumour RNA-seq data and clinical information were obtained from two separate series of patients.  Supplementary Tables 3 and 4. Data obtained from TCGA and ICGC databases are freely available to the public, and this research strictly followed access policies and publication guidelines. Therefore, this study did not require ethical review or approval from an Ethics Committee.

Gene set enrichment analysis (GSEA)
GSEA analysis was performed using GSEA software (version 4.0.3) to detect the difference in enriched pathways between the high-risk and low-risk subtypes by the Molecular Signatures Database (MSigDB, h.all.v7.2.symbols.gmt). For each analysis, gene set permutations were performed 1,000 times.

Estimation of immune cell type
We used the single-sample GSEA (ssGSEA) algorithm to quantify the relative abundance of infiltrated immune cells. The gene set stores a variety of human immune cell subtypes, including T cells, dendritic cells, macrophages, and B cells [31,32]. The enrichment score calculated using ssGSEA analysis was used to assess infiltrated immune cells in each sample.

Statistical analysis
Relationships among the m6A regulators were calculated using Pearson's correlation based on gene expression. Continuous variables are summarised as mean ± standard deviation (SD). Differences between groups were compared using the Wilcoxon test, using the R software. Different m6A-risk subtypes were compared using the Kruskal-Wallis test. The 'ConsensusClusterPlus' package in R was used for consistent clustering to determine the subgroup of A-HCC samples from TCGA. The Euclidean squared distance metric and K-means clustering algorithm were used to divide the sample from k = 2 to k = 9. Approximately 80% of the samples were selected in each iteration, and the results were obtained after 100 iterations [33]. The optimal number of clusters was determined using a consistent cumulative distribution function graph. Thereafter, the results were depicted as heatmaps of the consistency matrix generated by the 'heatmap' R package.
We then used Kaplan-Meier analysis to compare the disease-specific survival (DSS), disease-free interval (DFI), progression-free interval (PFI) or overall survival (OS) between different subtypes using the 'survival' and 'survminer' packages in R software. The significance of differences in survival time was calculated using the log-rank test with a threshold of p < 0.05. Univariate and multivariate analyses were performed using Cox regression, followed by identification of independent risk factors for DSS, DFI, PFI, and OS in A-HCC. To evaluate the accuracy and sensibility of the model, we constructed the receiver operating characteristic (ROC) curve and calculated the area under the curve (AUC) using the 'survivalROC' package in R software.

Regulatory pattern of m6A-related genes in A-HCC
The study design is shown in Figure 1. To determine whether the clinical prognosis of A-HCC is associated with known m6A-related genes, we summarised the occurrence of 21 m6A regulatory factor mutations in A-HCC in TCGA database (n = 117). Among them, VIRMA (KIAA1429) had the highest mutation rate (20%), followed by YTHDF3, whereas four genes (YTHDF1, ELAVL1, ALKBH5, and RBM15) did not show any mutation in this sample ( Figure 2A). To systematically study all the functional interactions between proteins, we used the web site GeneMANIA to construct a network of interaction between the selected proteins and found that HNRNPA2B1 was the hub of the network ( Figure  2B-C). Furthermore, we determined the difference in the expression levels of the 21 m6A regulatory factors between A-HCC and normal liver tissue ( Figure  2D-E). Subsequently, we analysed the correlation of the m6A regulators ( Figure 2F) and found that the expression patterns of m6A-regulatory factors were highly heterogeneous between normal and A-HCC samples, suggesting that the altered expression of m6A-regulatory factors might play an important role in the occurrence and development of A-HCC.

An integrative m6A risk model
To explore the prognostic value of the expression levels of the 21 m6A methylation regulators in A-HCC, we performed univariate Cox regression analysis based on the expression levels of related factors in TCGA dataset and found seven related genes to be significantly related to OS (p < 0.05), namely YTHDF2, KIAA1429, YTHDF1, RBM15B, LRPPRC, RBM15, and YTHDF3 (Supplementary Table  5). To identify the most powerful prognostic m6A regulator, we performed LASSO Cox regression analysis. Four candidate genes (LRPPRC, KIAA1429, RBM15B, and YTHDF2) were selected to construct the m6A risk assessment model ( Figure 3A, B), the risk score was as follows: risk score = (0.0648970639115386 × KIAA1429) + (0.0370948653489106 × LRPPRC) + (0.000459715556466468 × RBM15B) + (0.0605157571421274 × YTHDF2). Based on the expression levels of these four m6A-related genes as well as k = 2, a parameter that leads to Supplementary Table clustering outcome, we identified two new clusters in TCGA dataset ( Figure 3C-E). Principal component analysis showed that cluster analysis could successfully divide A-HCC patients into two subtypes ( Figure 3F). We compared the clinical survival curves of the two subtypes and found that the survival trend of subtype C1 was significantly better than that of subtype C2 (p = 9.832e-04; Supplementary Table 6, Figure 3G, Figure S1A). The expression levels of the four selected m6A-related genes and the clinicopathological variables in the two subtypes were closely related to tumour stage and grade ( Figure 3H).
We verified the gene and protein expression of the four m6A regulators screened in the collected samples from HCC clinical patients, and the results showed that compared with normal patients, KIAA1429, LRPPCC, RBM15b and YTHDF2 were up-regulated in HCC patients, which was more significant in A-HCC patients ( Figure S1B-C). Meanwhile, to further illustrate the external applicability of the model, we conducted survival analysis of the m6A model in a variety of cancers in addition to A-HCC and found that it was predictive (p =0.003), such as liver hepatocellular carcinoma (LIHC, p =0.01), lower grade glioma (LGG, p =0.029), uterine corpus endometrial carcinoma (UCEC, p =0.033) kidney chromophobe (KICH, p =0.005) and arenal cortical carcinoma (ACC, p =0.044; Figure S1D).
To further unravel the mutation events associated with the m6A risk model, we divided the A-HCC patients into high-risk and low-risk subtypes. In the high-risk subtype, 53% of the samples had mutations in TP53 ( Figure 3I), whereas CTNNB1 mutations were  Figure 3J). TP53 is a common tumour suppressor gene, and its mutations accompany tumorigenesis [34]. The frequency of TP53 mutations in the high-risk subtype was significantly higher than in the low-risk subtype (53% vs. 23%, p = 0.001; Figure 3K). Subsequently, we divided the A-HCC patients into two subtypes according to the presence or absence of mutations in TP53 ( Figure 3L). Risk scores and model-related gene expressions were higher in the TP53-mutation group than in the non-mutated group.
To explore the function of the four identified m6A-related genes, we extracted and screened genes their co-expressed genes and performed gene ontology (GO) enrichment analysis. A total of 202 genes were co-expressed with the four m6A-related genes ( Figure 3M) and their functional categories were molecular function (MF), biological process (BP), and cellular component (CC). These terms were mainly enriched in pathways related to RNA processing, modification, and proliferation such as ncRNA metabolic processing and regulation of lipid metabolic processes ( Figure 3N). Altogether, the results suggest that TP53 mutation may be a key factor in initiating m6A methylation, which activates cancer-promoting pathways. Hence, the expression levels of KIAA1429, LRPPRC, RBM15B, and YTHDF2 could be used as a prognostic indicator in A-HCC.

Prognostic performance of the m6A risk model in A-HCC
Since the expression levels of the four selected genes (KIAA1429, LRPPRC, RBM15B, and YTHDF2) play a crucial role in tumorigenesis and tumour development, we employed them to establish an m6A risk signature model. KIAA1429, YTHDF2, and RBM15B expression levels were not significantly different between the high/low-risk subtypes when DFI was analysed. In the LRPPRC low-expression group and low m6A risk model, patients had significantly longer DFI than those in the high-expression and high-risk model ( Figure 4A, Figure  S2A). To demonstrate the reliability of the m6A risk model, we constructed an ROC curve for DFI prediction and quantified the AUC. The AUC of the m6A risk model in 1/2/3 years was better than that of the expression of other genes (KIAA1429, YTHDF2, and RBM15B) and other factors, such as age, sex, and tumour grade ( Figure 4B-C). The clinical prognostic differences were consistent from DSS and PFI analysis ( Figure 4D-I), which indicated that the m6A risk model composed of four genes (KIAA1429, LRPPRC, RBM15B, and YTHDF2) can more accurately predict the prognosis of A-HCC. Analysis of KIAA1429, LRPPRC, RBM15B, and YTHDF2 expression levels in the high-and low-risk subtypes from TCGA database showed significantly upregulated expression in the high-risk subtype ( Figure 5A). The high-risk subtypes had a lower OS and a higher risk score than those in the low-risk subtype ( Figure 5B-C). Higher expression levels of KIAA1429, and RBM15B and higher m6A risk model scores were associated with a higher mortality rate in the high-risk subtype ( Figure 5D, Figure S2D). To further evaluate the accuracy of the m6A risk model for predicting the 1, 2, and 3-year survival rate of A-HCC patients, we performed ROC curve analysis on TCGA (n = 117) cohorts ( Figure 5E). Similarly, the performance of the m6A risk model was better than the models using the expression levels of a single gene and other factors (age, gender, tumour grade, tumour stage, and vascular invasion; Figure 5F-G).
Meanwhile, the same verification was performed in the ICGC database ( Figure S3). The above data show that the m6A risk model predicts the OS of A-HCC patients with more accuracy and reliability than any of the other models analysed.

m6A risk model to evaluate the occurrence and development of A-HCC
Considering that m6A methylation is closely related to the occurrence and development of tumours, we explored the relationship between the m6A risk model and clinicopathological characteristics. In TCGA cohort, the expression levels of LRPPRC and RBM15B and the m6A risk score were significantly correlated with tumour grade and T stage. Increases in tumour grade and T stage were associated with higher m6A risk score or gene expression levels. Additionally, KIAA1429, LRPPRC, and RBM15B and the m6A risk scores were significantly different between different tumour stages ( Figure 6A). Subsequently, we validated our conclusions again using the ICGC dataset. In the ICGC cohort, KIAA1429, LRPPRC, RBM15B, and YTHDF2 expression levels and the m6A risk score were significantly correlated with tumour grade. Moreover, the increase in tumour grade was associated with a gradual increase in the m6A risk model score. Only RBM15B expression levels and the m6A risk model score were associated with tumour stage and T stage. We also evaluated the relationship between the m6A risk model and vascular invasion and found that KIAA1429 and LRPPRC expression levels and the m6A risk model score were significantly correlated with vascular invasion ( Figure  6B). This indicates that tumour vascular invasion is highly correlated with the model score and that patients with higher scores are more likely to exhibit vascular invasion.
Next, we used the Cox regression model to perform univariate and multivariate survival analyses on the m6A risk model. In TCGA dataset, both univariate and multivariate analyses showed that tumour stage and the m6A risk model score were strongly associated with OS ( Figure 6C), which was replicated in the ICGC database ( Figure 6D). Thus, we concluded that the m6A risk model can used to evaluate the occurrence and development of A-HCC.

GSEA signalling pathways
To further explore the pathways potentially involved in the development of A-HCC, we divided the patients from TCGA and ICGC databases into high-risk and low-risk subtypes based on risk scores and performed GSEA enrichment analysis (Supplementary Table 7). Pathways enriched in the high-risk subtype were mainly related to tumour formation and proliferation, such as E2F targets, DNA repair, and MTORC1 signalling pathways ( Figure  7A). Interestingly, the enriched pathways were shown to be closely related to tumour development and anti-apoptosis. For example, the E2F pathway plays a key role in cell proliferation by regulating the cell cycle [35].

Utility of the m6A risk model in diagnosing and assessing the disease status of A-HCC
To explore the potential role of the m6A risk model in the diagnosis of A-HCC as well as its reliability and accuracy, we compared it with known A-HCC-related genes and diagnostic markers. Alpha-fetoprotein (AFP) is the most commonly used clinical HCC marker [36]. Other proteins closely related to A-HCC include patatin-like phospholipase domain-containing protein 3 (PNPLA3), hydroxysteroid 17-beta dehydrogenase 13 (HSD17B13), serpin family A member 1 (SERPINA1), and transmembrane 6 superfamily member 2 (TM6SF2) [37][38][39][40]. We found that the m6A risk model (AUC = 0.888) had a better predictive accuracy for A-HCC diagnosis compared with that of AFP, SERPINA1, TM6SF2, and PNPLA3 expression levels ( Figure 7B). We validated these results using the ICGC database ( Figure 7C).
We next evaluated the specificity of the m6A model in distinguishing A-HCC from alcoholassociated non-malignant changes. Surprisingly, the m6A risk model score was significantly increased in the A-HCC samples compared with A-HCC paracarcinoma and normal tissue samples in both TCGA and ICGC databases; additionally, the m6A model showed a marked sensitivity in A-HCC diagnosis ( Figure 7D-E). We also verified that this model was superior to other related factors in distinguishing cancer and paracarcinoma tissue samples ( Figure S4), demonstrating that dysregulated expression levels of m6A genes are highly specific in the tumorigenesis of A-HCC. Based on TCGA database, we established a nomogram to predict the OS of patients according to various possible influencing factors ( Figure 7F).

The immune landscape of A-HCC patients
To explore the immune landscape of A-HCC patients, we performed ssGSEA using TCGA and ICGC databases. The resulting heatmap is shown in Figure S5A-B and the infiltration levels of various immune cell types are shown in Figure 8A. The infiltration levels of most immune-activated cells, such as activated CD8 + cells, activated CD8 + T cells, effector memory CD8 + T cells, gammadelta T cells, and immature B cells, were reduced in the high-risk subtype. However, the proportion of activated CD4 + T cells and CD56 dim natural killer cells in the high-risk subtype was higher than in the low-risk subtype. We then found a positive correlation between these immune cells, and the proportion of myeloid-derived suppressor cells were closely correlated with that of regulatory T cells (R = 0.91; Figure 8B). Subsequently, we downloaded the immunosuppressive cytokines related to the cancer-immunity cycle from the Tracking Tumour Immunophenotype website [41] and compared the relationship between the m6A risk model and immunosuppressive cytokines using box plots. The results showed that levels of most of the immunosuppressive cytokines, such as Arg2, CCL28, DNMT1, and EZH2, were upregulated in the highrisk subtype (Figure 8C), suggesting that high-risk A-HCC patients have reduced cancer-immunity cycle activity. Similarly, we analysed the correlation between these immunosuppressive cytokines and found that DNMT1 and EZH2 were highly correlated (R = 0.71; Figure 8D). Kaplan-Meier analysis of DNMT1 and EZH2 showed that patients with higher DNMT1/EZH2 expression have poorer OS ( Figure  8E-F). To verify the above conclusions, we generated a Venn diagram of immune cells and immunosuppressive cytokines of TCGA/ICGC databases, which resulted in an overlap of 19 immune cells and 26 immunosuppressive cytokines ( Figure  S4C). Subsequently, we explored the correlation between immunosuppressive cytokines (DNMT1 and EZH2) and all immune cells (Supplementary Table 8) and found that three types of immune cells (activated CD4 + T cells, monocytes, and neutrophils) were closely related to DNMT1 and EZH2 levels ( Figure  8G-H). These results were further validated in the ICGC database ( Figure S6A-H). Therefore, we performed logistic regression analysis of the model risk score and immune cells/immunosuppressive cytokines levels and found that they were closely correlated ( Figure S7). Altogether, these results indicate that an increase in activated CD4 + T cell infiltration is associated with higher expression levels of DNMT1 and EZH2, whereas the opposite was observed for monocyte and neutrophil infiltration. Therefore, immunosuppressive cytokines, such as DNMT1 and EZH2, and immune cells, such as activated CD4 + T cells, monocytes, and neutrophils, may form a TIM regulatory system, representing a new target for A-HCC therapy.

m6A model predicts A-HCC treatment efficacy
In TCGA database, patients in the m6A high-risk subtype had lower immune and stroma scores as well as lower ration immune score -stroma Score/microenvironment score than patients in the m6A low-risk subtype ( Figure 8I). Thus, our model could predict the TIM state and the therapeutic responses of A-HCC. Recently, an ImmuCellAI estimation method was developed to predict the response of HCC patients to immunotherapy [42]. We evaluated whether the m6A risk model can make similar predictions and analysed the difference in KIAA1429, LRPPRC, RBM15B, and YTHDF2 expression levels and the risk score between the responder and non-responder group.
Significantly upregulated expression of KIAA1429, LRPPRC, and RBM15B and high-risk scores were observed in the non-responder group compared with the responder group ( Figure 8J). This was further verified using the ICGC database ( Figure  S5I-J). As shown in Figure 8A, most high-risk subtypes lacked immune cells; immunoreactive cell deficiency is known to cause immunotherapy tolerance [43,44], which indicates that high-risk may be related to non-response to immunotherapy (immune tolerance).

Validation of A-HCC core genes (DNMT1/EZH2) and potential drugs
We collected pathological samples from normal, N-A-HCC, and A-HCC patients and performed immunohistochemical staining and qRT-PCR. DNMT1 and EZH2 levels in the liver tissues of normal individuals and N-A-HCC patients were barely detecSupplementary Table, while they were diffusely expressed in A-HCC patients ( Figure 9A-C), indicating that DNMT1 and EZH2 expression in A-HCC patients is increased in comparison with normal and N-A-HCC individuals.
We then evaluated the role of DNMT1 and EZH2 in guiding A-HCC treatment. As the therapeutic effects of PX-12 [45], LRRK2-IN-1 [46], and GSK-J4 [47] in A-HCC have been already described, we decided to explore the therapeutic effect of teniposide on A-HCC. We employed two HCC cell lines, Huh7 and HepG2, and treated them with 100 mM alcohol, as a cellular model of A-HCC. DNMT1 and EZH2 gene expression and protein levels, evaluated by qRT-PCR, western blotting and immunofluorescence staining, were significantly higher in the alcohol-treated group (100 mM) than in the control group. Administration of teniposide (0.5 μM) to alcohol-treated cells abolished these effects ( Figure  9D-F). Given that DNMT1 and EZH2 are barely expressed in the control group but are significantly up-regulated by alcohol-treatment and significantly down-regulated after teniposide treatment, the results suggest that DNMT1 and EZH2 may be core proteins in the aetiology of A-HCC and highlight teniposide as a potential therapeutic drug.

The potential therapeutic effect of teniposide against A-HCC in vivo
We evaluated the role of teniposide in the occurrence and development of HCC in mice; an overview of the experimental procedure is provided in Figure 10A. Mice begin to form HCC 7-10 months after injection of DEN solvent [48]; hence we administered alcohol and drugs (teniposide, TEN) at 6 months of age, and divided the mice into five groups: Control+NC (without TEN and alcohol), Alcohol+NC (without TEN), Alcohol+TEN, Control+TEN (without alcohol) under DEN stress, and Control without DEN stress (10 mice per group). MRI imaging (AVANCE IIITM HD 600MHz) of the mouse liver was obtained at 8 months and representative liver images of each group are shown in Figure 10B. Tumour number analysis showed that teniposide significantly reduced tumours numbers in A-HCC ( Figure 10C). Haematoxylin and eosin staining of liver sections demonstrated that Alcohol+NC group had the most obvious liver lesions and that teniposide was more effective in treating A-HCC than HCC ( Figure 10D). We then determined the expression of two A-HCC core genes (DNMT1 and EZH2) using qRT-PCR and IHC, and found significantly higher expression in the A-HCC group than in the HCC group, which significantly decreased following teniposide treatment ( Figure 10E-G). Taken together, these data suggest that teniposide has a potential therapeutic effect on the occurrence and progression of A-HCC by acting on the A-HCC core genes, DNMT1 and EZH2.

Discussion
Increasing evidence has demonstrated that the interaction between multiple m6A regulators plays an essential role in the development and progression of several types of cancers. Here, we summarised the m6A-regulatory genes involved in the pathways associated with tumorigenesis (Supplementary Table  9). To clarify the relationship between m6A-related genes and the prognoses of patients with A-HCC, we selected 21 m6A regulators and mapped the m6A modifications mediated by these regulators and their potential biological functions in disease occurrence ( Figure 11). Demethylases (FTO and ALKBH5) and methyltransferases (such as Metl3 and Metl14) have been reported to regulate the progression of several types of cancers, including liver, lung, and breast cancers [49][50][51][52]. For example, silent information regulator 1 (SIRT1) can deregulate FTO to guide the m6A methylation of downstream molecules [53], and ALKBH5 can act as a tumour suppressor by reducing the expression of LYPD1 in HCC [54].  Currently, the mortality rate of HCC remains high, especially for A-HCC patients, owing to the lack of early diagnostic markers, treatment resistance, and poor immunotherapy sensitivity [55]. The TIM plays a decisive role in the effectivity of treatment methods and can lead to tumour development and recurrence after treatment [56,57]. In the progression of liver cancer, immune cells and immunosuppressive cytokines can adjust the balance of the immune system by, for example, modifying the proportion of T-cell subsets [58]. Importantly, the immunosuppressive state in the tumour can be transformed into an immune-activated state by adjusting the TIM and inducing an effective immune response, which can lead to an enhanced sensitivity to immunotherapy. Therefore, evaluating the TIM, and its relationship with prognosis in A-HCC patients is urgent. In the present study, we clarified the role of m6A methylation in regulating the TIM and developed a risk model based on the expression levels of m6A regulators, which may help predict the prognosis of A-HCC patients and more effective immunotherapy targets.
The risk model developed herein integrated four m6A methylation regulators, LRPPRC, YTHDF2, KIAA14219, and RBM15B. LRPPRC encodes a nuclear protein that is significantly negatively correlated with the immune response in a variety of tumour tissues [59][60][61]. YTHDF2 encodes a member of the YT521-B homology (YTH) domain family that has been reported to regulate mRNA stability [62,63], and Yu et al. found that YTHDF2 plays a negative regulatory role in the inflammatory response induced by LPS [64]. KIAA14219 encodes a methyltransferase that acts as a m6A 'writer' [65] and has been reported to be involved in tumour cell proliferation and metabolism pathways in liver and breast cancers [66,67]. Finally, RBM15B, which encodes another m6A 'writer', modifies m6A by binding to m6A methylation complexes (METTL3 and WTAP) [68]. RBM15B has been reported to be associated with the immune landscape in various diseases [69].
In this study, we used the four m6A regulators to divide A-HCC patients into two subtypes and predicted their prognosis, and the model was validated in clinical patient samples we collected. We notably found that m6A high-risk subtypes had a high frequency of mutations in TP53. As TP53 is a tumour suppressor gene, this indicates that TP53 mutations may cause changes in m6A methylation levels. In addition, the pathways associated with the high-risk subtype were mainly related to RNA processing modification, and tumour development, suggesting that these four m6A regulators can be used as indicators of the occurrence and prognosis of A-HCC. In analysing different survival interval (DFI, DSS, PFI and OS), we found that the prognosis of the m6A high-risk subtype was significantly worse and that the m6A risk model was more reliable and accurate than single genes in prediction efficiency, which could be used as an independent predictor. Meanwhile, the model was more reliable than the common clinical indicators AFP, PNPLA3, HSD17B13, SERPINA1, and TM6SF2 in predicting patient outcome. Finally, we constructed a nomogram based on various confounding factors, with the aim of applying this model to clinical guidance in the future.
GSEA indicated that the pathways enriched in the high-risk subtype were related to tumour formation and proliferation, which included the common E2F pathway and the PI3K/Akt/mTOR pathway [70,71]. E2F is a transcription factor that controls the expression of all cell division genes, of which E2F8 is significantly increased in HCC and ovarian cancer [72]. It can transcriptionally inhibit CDK1-induced hepatocyte polyploidy, interact with HIF1 to form a complex, improve VEGFA level, promote angiogenesis, and induce tumour metastasis [72,73]. In addition, the PI3K/Akt/mTOR pathway is essential for tumour survival and growth, and induces resistance to radio-therapy, chemo-therapy, and cytostatic drugs [74].
A large amount of data from various disease conditions have indicated a correlation between m6A modifications and TIM [75][76][77]. Although several studies have investigated the role of single regulatory factors or a single immune-infiltrating cell type in the immune response [78,79], the comprehensive role of multiple m6A regulators in the immune response has not been studied to date. In this study, we describe the relationship between m6A regulators and the A-HCC immune response. In our model, there were clear differences in the TIM cell infiltration characteristics, higher m6A risk scores were associated with a higher infiltration of activated CD4 + T cells, higher levels of immunosuppressive cytokines (DNMT1 and EZH2) and reduced levels of monocytes and neutrophils infiltration.
These features indicate an immunosuppressive TIM in the high-risk subtype, corresponding to the so-called 'immune desert type'. In contrast, the low-risk subtype had an immuneactivated state. Therefore, the immunosuppressive cytokines DNMT1 and EZH2; and the immune cells activated CD4 + T cells, monocytes, and neutrophils appear to form a TIM regulatory system that greatly impacts the prognosis of A-HCC. DNMT1, a common DNA methyltransferase, is involved in DNA methylation in eukaryotes [80]. DNMT1 is closely related to the occurrence and development of various diseases, including several types of cancers, as it affects the methylation levels of CD4+T cell-related genes, thereby inhibiting the immune response [81][82][83][84]. EZH2 acts as a catalyst for polycomb repressive complex 2 (PRC2) formation, catalysing the trimethylation of lysine 27 on histone H3 (H3K27me3) and mediating gene silencing [85]. Several studies have reported that EZH2 can regulate the development and function of B cells and neutrophil migration and change the plasticity of CD4+T cells, highlighting the important role of EZH2 in the immune regulation of many diseases [86][87][88].
CD4 + T cells act as central orchestrators of immune regulation. Depending on the specific TIM, activated CD4 + T cells can differentiate into CD4+ T helper (Th) cells, which collaborate with B cells and CD8+ T cells promote immune response [89,90]. Monocytes are an important part of innate immunity and have been reported to be key regulators of cancer development [91]. During tumorigenesis, monocytes perform several antitumor immunity functions, including phagocytosis and recruitment of lymphocytes, and can even differentiate into tumour-related immune cells [92,93]. Neutrophils exhibit powerful antimicrobial functions, including phagocytosis and formation of neutrophil extracellular traps [94,95]. Under pathological conditions, neutrophils are activated and infiltrate lesions, thereby changing the tissue microenvironment [96][97][98].
We evaluated the performance of the m6A risk model in assessing the sensitivity of immunotherapy and found that high score models were associated with reduced sensitivity to treatment. This may be because activated CD4 + T cells, monocytes, and neutrophils in the m6A high-risk subtype interact with DNMT1 and EZH2, resulting in an immunosuppressive, desert type microenvironment. DNMT1 and EZH2 expression levels were then compared between normal, N-A-HCC and A-HCC samples, while activating activated CD4T cells and inhibiting monocyte and neutrophil. DNMT1 and EZH2 expression levels were revealed to be correlated with changes in immune cells in the TIM and may improve the TIM state by inhibiting its expression. Through drug sensitivity analysis, we found that A-HCC patients were generally sensitive to teniposide, PX−12, LRRK2−IN−1, and GSK−J4 drugs, which can help clinicians better select treatment strategies. Among these four drugs, teniposide has not been reported in HCC studies. In our study, we found that teniposide has a potential therapeutic effect on A-HCC by down-regulating the expression of A-HCC core genes (DNMT1 and EZH2), thereby reversing the malignant degree of A-HCC and improving the prognosis.
In conclusion, we employed the expression levels of m6A regulators to construct a risk model that can accurately predict the prognosis of A-HCC patients and aid further understanding of the TIM state in A-HCC. The model can also predict the sensitivity of A-HCC patients to immunotherapy and drug therapy, which can greatly help guide future clinical selection of A-HCC targeted therapy and immunotherapy. Our finding also demonstrated that DNMT1 and EZH2 can be exploited as core genes of A-HCC and that teniposide can be used for the treatment of A-HCC.