Abstract Given the resistance to conventional treatments and limitations of immune checkpoint blockade therapy in hepatocellular carcinoma (HCC), it is imperative to explore novel prognostic models and biomarkers. The dependence of cancer cell on exogenous methionine, known as Hoffman effect, is a hallmark of HCC, with numerous studies reporting a strong correlation between methionine metabolism and tumor development. Betaine-homocysteine S-methyltransferase (BHMT), a critical component of methionine metabolism pathway, has polymorphisms linking to poor prognosis in multiple cancers. Nevertheless, there is little literature regarding the relationship between methionine metabolism and incidence, mortality of HCC, as well as the function of BHMT in HCC progression. In this study, by analyzing multiple datasets, we constructed a methionine metabolism-related prognostic model and thoroughly investigated the influence of BHMT on the prognosis of HCC. Bioinformatics analysis revealed a marked decrease in BHMT expression in HCC, which was linked to adverse clinical outcomes. CIBERSORT results suggest that BHMT promotes infiltration of M1 macrophages. Our results suggest its potential as an ideal prognostic biomarker for anti PD-L1 immunotherapy. In summary, this study innovatively provides first methionine metabolism-related prognostic model and unveils the tumor suppressive function of BHMT in HCC, providing potential mechanism by which BHMT exert its function. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-93650-w. Keywords: Methionine metabolism, Lipid metabolism, Hepatocellular carcinoma, BHMT, Tumor microenvironment, Immunotherapy Subject terms: Cancer metabolism, Cancer microenvironment, Cancer therapy, Gastrointestinal cancer, Tumour immunology Introduction The global incidence of hepatocellular carcinoma (HCC) has seen a gradual increase recently^[50]1, positioning it as the world’s sixth most frequently diagnosed neoplasm and the third primary contributor to cancer-induced fatalities^[51]2. In 2022, there were more than 860,000 new cases of primary liver cancer and over 750,000 deaths, with hepatocellular carcinoma accounting for 80% of all primary liver cancer cases^[52]2. The risk factors for HCC include HBV infection, obesity, patterns and duration of alcohol consumption, genetics, and dysbiosis of gut microbiota^[53]3. The application of high-definition imaging technology and specific molecular biomarkers has benefited the screening of HCC^[54]4. However, due to the covert onset, rapid progression, and ease of metastasis and recurrence characterizing HCC^[55]5, most patients are diagnosed at a late stage of the disease^[56]6. Early detection of HCC, followed by timely treatment, can result in a five-year survival rate exceeding 70%^[57]7. Over the past five years, with an enhanced understanding of the tumor microenvironment and the development of emerging treatments, immunotherapy, especially immune checkpoint inhibitors, has been increasingly applied in the therapeutic approach to advanced HCC^[58]8. However, treatment outcome varies widely among patients, which may be related to the heterogeneity of the tumor microenvironment in different individuals. The aim of this study is to work towards the development of more sensitive and accurate prognostic models, and to search suitable biomarker which can also predict the outcome of immunotherapy of HCC. Tumor metabolic reprogramming is a key feature of HCC and involves alterations in a variety of metabolic pathways including glycolysis, pentose phosphate pathway and lipid metabolism^[59]9. Methionine metabolism is one of the liver-specific metabolic pathways, with the key enzyme methionine adenosyltransferase (MAT) being almost exclusively expressed in the liver^[60]10. The liver is crucial for the synthesis of S-adenosylmethionine (SAM), with approximately 85% of SAM being synthesized in this organ. It has been reported that in models of liver cirrhosis and primary liver cancer, SAM levels are significantly reduced^[61]11. Hoffman effect, the dependency of cancer cells on exogenous methionine, suggests that even when provided with homocysteine, a direct precursor of methionine, cancer cells are unable to synthesize methionine^[62]12. This indicates the deficiency of SAM caused by a disruption in the methionine synthesis pathway. Previous studies using mouse models and cohort studies have confirmed that reducing methionine intake in animals or humans can effectively inhibit tumor development and enhance the efficacy of antimetabolite drugs in chemotherapy and radiotherapy^[63]13. A prospective study has shown that limiting dietary methionine intake can effectively reduce the all-cause mortality rate of invasive breast cancer^[64]14. Despite numerous investigations into the process of methionine metabolism and its association with tumor development, there is a scarcity of research exploring the correlation between methionine metabolism and the occurrence and mortality of HCC. Betaine-homocysteine S-methyltransferase (BHMT), an enzyme predominantly expressed in the kidney and liver, regulates the activity of the betaine pathway, the alternative synthesis route of methionine^[65]15. BHMT plays a crucial role in preventing homocysteine accumulation, regulating renal osmolality, and promoting adipocyte differentiation and insulin resistance^[66]15–[67]18. Altered expression of BHMT has been proved with multiple cancer prognosis. The BHMT G742A polymorphism increases the risk of head and neck squamous cell carcinoma^[68]19, while the BHMT Arg239Gln polymorphism is associated with an increased risk of colorectal cancer^[69]20. Interestingly, studies using BHMT knockout (BHMT-/-) mice have shown BHMT deficiency leads to a higher incidence of fatty liver and HCC^[70]21. Besides, in HCC, downregulation of BHMT expression has been reported and is associated with poor prognosis^[71]22. Nevertheless, there are still limited reports on the potential pathways through which BHMT exerts its effects. In this study, as manifested in Fig. [72]1, combing clinically relevant prognostic indicators and key genes from the cysteine and methionine metabolic pathways, a prognostic computing model for HCC was conducted through the analysis of multiple datasets from the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). The model demonstrated excellent predictive capabilities for patient prognosis. Moreover, BHMT, a critical gene within liver methionine synthesis, was subjected to independent analysis. Our bioinformatic results indicate a significant downregulation of BHMT in HCC, and its expression level was associated with poor prognosis. Through comprehensive analysis of BHMT expression levels in relation to staging and prognosis, we found that BHMT expression is negatively correlated with the progression of HCC. Further analysis revealed that BHMT is implicated in pathways affecting fatty acid metabolism and in the reshaping of the tumor microenvironment. Intriguingly, it also exhibited remarkable predictive accuracy for the outcomes of immunotherapy in HCC. Besides, BHMT expression was positively associated with Peroxisome Proliferator-Activated Receptor (PPAR) signaling and negatively associated with cytokine-cytokine receptor interaction. In summary, we constructed first methionine metabolism related HCC prognostic model and unveiled the pivotal role of the key gene BHMT in tumor development and sculpting of tumor microenvironment. Fig. 1. [73]Fig. 1 [74]Open in a new tab Overall design of the study. Materials and methods Patient data This study compiled data from 96 HCC patients, pathologically confirmed at Nanchang Ninth Hospital between 2008 and 2018. The cohort comprised 82 males and 14 females, ranging in age from 27 to 70 years, with a mean age of 47 ± 9.8 years. Preoperative serum levels of Aspartate Aminotransferase (AST), Alanine Aminotransferase (ALT), Cholesterol (CHO), High-Density Lipoprotein (HDL), Low-Density Lipoprotein (LDL), Triglyceride (TG), Apolipoprotein A1 (APOA1), Apolipoprotein B (APOB), and lipoproteins were collected for these HCC patients. Ethics statement This research was reviewed and approved by Medical Ethics Committee of Infectious Diseases Hospital of Nanchang University. All experiments attached to the regulations and guidelines of Nanchang University. The participants provided their written informed consent to participate in this study. Data acquisition and preprocessing The GDC-TCGA-LIHC dataset from Xena ([75]https://xena.ucsc.edu/) was downloaded to provide hepatocellular carcinoma mRNA matrices in FPKM and count formats from TCGA. Survival data for the samples in the mRNA matrices, containing overall survival (OS) and disease-specific survival (DSS), were obtained from Xena Pan-Cancer Atlas, which contains de-batching data from TCGA and Genotype-Tissue Expression (GTEx). Additionally, the R package ‘TCGAbiolinks’ was utilized to download corresponding clinical data, encompassing detailed tumor staging and tumor size. Only samples with available mRNA expression matrices, OS, and DSS were included in subsequent analyses. In survival analyses, if multiple samples from the same patient were present, only one sample was retained for analysis, prioritizing frozen samples. GEO ([76]https://www.ncbi.nlm.nih.gov/geo/) provided mRNA expression matrices and clinical treatment outcomes in [77]GSE202069, a HCC cohort utilizing anti-PD1 immunotherapy, as well as mRNA expression matrices and associated clinical data in [78]GSE14520 and [79]GSE116174. mRNA expression matrix and clinical data of ICGC-LICA-FR cohort was obtained from International Cancer Genome Consortium (ICGC) ([80]https://dcc.icgc.org/). Subsequently, Transcripts Per Kilobase Million (FPKM) formatted matrices were converted to Transcripts Per Million (TPM) format, and all mRNA matrices were processed into log(x + 1) format to ensure successful data visualization. Identification of candidate genes for prognostic modeling To identify candidate genes for prognostic modeling, differential analysis was conducted on HCC tumor and normal liver tissues from TCGA and GTEx data utilizing R package ‘limma’. Only genes expressed above 30 counts in 50% of samples were considered, with selected thresholds of log[2] |Fold Change| > 1 and p value < 0.05. An intersection of the differential genes and the cysteine and methionine metabolism pathway geneset from Kyoto Encyclopedia of Genes and Genomes (KEGG) yielded candidate genes for the prognostic model^[81]23. Establishment and verification of the prognostic model The advancement of machine learning has laid the groundwork for the steady improvement in the accuracy of prognostic models. To obtain the most accurate model possible, previously reported methods^[82]24 were employed. This approach integrates multiple models and trains them across various datasets to achieve optimal performance, using TPM data that were log2-transformed and Z-score normalized. Specifically, we utilized 30 different machine learning algorithms and 13 candidate genes, with OS as the outcome event. The optimal model was selected by comparing the average C-indices in the TCGA, [83]GSE116174, [84]GSE14520-[85]GPL3921, and ICGC-LICA-FR datasets. To ensure robustness, model training employed 10-fold cross-validation to determine the parameters for regularized models and boosting methods, while the random forest model was trained by constructing 1,000 decision trees with default parameters. The C-indices were calculated using the Cox proportional hazards model based on risk scores generated by each model. Notably, integrating machine learning models with widely used clinical indicators and gene expression-based models proved essential for building high-quality prognostic models^[86]25. To further enhance the accuracy of our model, based on TCGA-LIHC dataset, we conducted univariate Cox regression analysis with threshold of p value < 0.05 to select clinically relevant prognostic indicators, identifying stage and HBV infection status. Multivariate Cox regression analysis was utilized to combine the optimized model with clinically relevant prognostic indicators. In summary, we first used an integrated machine learning approach to construct an optimal model based on transcriptomic features, which is well-suited for interpreting high-dimensional and complex data. Subsequently, to further improve model accuracy and incorporate clinical information, we employed multivariate Cox regression analysis to construct a secondary model by integrating clinical factors. To further investigate the superiority of the combined model, a clinical model incorporating only the selected clinical factors was constructed using the multivariate Cox regression method. The combined model was exhibited in form of nomogram by R package ‘regplot’. For primary and combined model, as well as clinical model, functional validation was performed using TCGA and [87]GSE14520-[88]GPL3921 datasets. Model’s ROC curves were plotted at 1, 3, and 5-year cutoff times. Expression analysis To elucidate the role of BHMT in tumors, the expression of BHMT was initially explored. ‘TCGAplot’^[89]26, a highly integrated R package based on TCGA data, was utilized with methods ‘pan_boxplot’ for analyzing and visualizing BHMT expression levels across cancers. BHMT expression boxplots in HCC tumor/adjacent normal tissues and paired hepatocellular carcinoma samples were generated using TCGA datasets downloaded from Xena. Furthermore, BHMT mRNA expression among different Gender, Grade and Stage were also plotted based on TCGA dataset. Survival analysis The predictive potential of BHMT expression was analyzed using Kaplan-Meier (KM) curve analysis in the TCGA-LIHC, GSE14520_GPL3921, [90]GSE116174 and ICGC-LICA-FR dataset with OS as outcome event. The ‘surv_cutpoint’ function from the R package ‘survminer’ was used to identify the optimal cutoff point based on BHMT expression levels. Tumor microenvironment analysis The immune infiltration levels in TCGA samples were inferred based on transcriptomic expression profiles using the R package ‘CIBERSORT’ with threshold p value < 0.05 for high-confidence inferences for subsequent analyses. Samples were classified into high and low BHMT expression groups using median BHMT expression as threshold, and violin plots were drawn to compare the infiltration levels of various immune cells between the groups. Furthermore, dot plots correlating BHMT with significantly different immune cell infiltrations were drawn to explore the association between BHMT mRNA level and the infiltration levels of specific immune cells. The “ESTIMATE” package was utilized to calculate the correlation between BHMT expression and stromal as well as immune cell infiltration in HCC. Previous research reported key immune checkpoints^[91]27,[92]28. To further investigate the correlation between BHMT and immune landscape of tumors, TCGA data were divided based on the median BHMT mRNA level, and boxplots of various key immune checkpoints were drawn. Pan-cancer cytokine and cytokine receptor expression level was analyzed with R package ‘TCGAplot’. In [93]GSE202069, samples from HCC immunotherapy were split into high and low BHMT groups according to the median BHMT expression and combined with clinical immunotherapy outcomes to evaluate the BHMT’s predictive functionality for immunotherapy outcomes. Pathway enrichment analysis Downregulated genes identified in differential expression analysis between normal and tumor samples were utilized to perform KEGG analysis. To focus on metabolism-associated pathways, only those pathways associated with three major metabolic processes were included in this analysis. Genesets for BHMT pathway enrichment were generated by differential analysis based on the R package ‘limma’ The analysis included samples expressing the highest and lowest 20% of BHMT in TCGA HCC samples, considering only genes with expression levels above 30 counts in at least 50% of samples. The analysis thresholds were log[2] |Fold Change| > 1 and p value < 0.05, identifying 1171 positively associated genes and 1065 negatively associated genes. Above associated genes were analyzed for Gene Ontology (GO) and KEGG pathways using the R package ‘clusterProfiler’. Besides, method ‘gseKEGG’ was also utilized to perform Gene Set Enrichment Analysis (GSEA) in TCGA-LIHC with differential expression results mentioned above. Furthermore, BHMT expression association analysis was performed with Enoyl-CoA Hydratase and 3-Hydroxyacyl CoA Dehydrogenase (EHHADH) and Acyl-CoA Synthetase Long Chain Family Member 1 (ACSL1). Preparation of liver paraffin tissue arrays Liver paraffin tissues of 96 HCC patients diagnosed pathologically at Nanchang Ninth Hospital from 2008 to 2018 were collected. After reviewing the HE staining results of all wax blocks, those containing HCC and adjacent normal tissues were selected, marked, and arranged in order of pathology numbers. Tissue-containing wax columns from the cancer and adjacent areas were extracted using a 2 mm diameter puncher, placed face-down in a 96-well plate in sequence. A cleaned embedding iron box was prepared, and the tissue-containing end of each wax column was gently dabbed with glue, moistened, tapped on paper to remove excess glue, and affixed to the iron box. Twelve columns per row were arranged in 12 rows. The embedding box with wax columns was placed on an embedding machine until the remaining wax melted, followed by embedding with additional wax. Immunohistochemistry Immunohistochemistry was performed on liver tissue arrays to verify the expression of BHMT, EHHADH, and ACSL1 proteins in HCC and adjacent normal tissues. Expression of Glypican 3 (GPC3), Hepatocyte, Alpha-Fetoprotein (AFP), Ki67, and other antibodies were tested on wax blocks containing both cancerous and adjacent normal tissues from HCC patients. All wax block samples underwent 4 μm sectioning, were baked at 68℃ for 2 h, and sequentially processed for antigen retrieval, blocking, primary antibody application (50–100 µl), overnight incubation at 4 °C, and secondary antibody application (50–100 µl) at 37 °C in PBS for 30 min. Staining, counterstaining, dehydration, and coverslipping were followed by microscopic observation and statistical analysis. Primary antibodies included anti-BHMT polyclonal antibody ab243698 (1:500 Abcam UK), anti-EHHADH polyclonal antibody 26570-1-AP (1:350, Proteintech USA), and others. The antigen retrieval solution used was EDTA (pH8.0, ZLI-9079, 50× dilution). Assessment of immunohistochemical results Immunohistochemical assessment scores were based on staining intensity and range for each sample. Staining intensity was graded as 0 (no staining), 1 (light yellow), 2 (brown-yellow), and 3 (brown). Staining range was rated as 0 (no staining), 1 (5-25% range), 2 (26-50% range), 3 (51-75% range), and 4 (> 75% range). The final assessment score for each sample was the product of staining intensity and range. A score below 1 indicated negative staining; otherwise, it was considered positive. Statistical analysis Statistical analysis was performed using R (version 4.2.2). Intergroup significance was generated by the t-test, while paired-sample intergroup significance was calculated by the paired t-test. The Kaplan-Meier method was used for survival analysis. Besides, Chi-square test, Fisher’s exact test, and Kruskal Wallis test were utilized to assess correlations between BHMT mRNA and HCC clinical pathology, serum biochemistry results, and other antibody relevancies in liver arrays. The significance of BHMT for predicting immunotherapy outcomes was calculated by the Chi-square test. Multiple comparison corrections were applied as follows: Benjamini-Hochberg (BH) correction for differential expression analysis conducted using the limma package, false discovery rate (FDR) adjustment for pathway enrichment and GSEA results, and Holm-Bonferroni adjustment for multiple pairwise expression analyses and Cox regression analyses. Results Identifying prognostically relevant differential genes in the methionine and cysteine pathway Differential analysis of 160 normal tissue samples and 369 HCC samples from the TCGA-LIHC dataset, revealed 2745 differential genes (Fig. [94]2a, Table [95]S1), with 1239 being upregulated and 1506 downregulated. The results of metabolic associated pathway enrichment analysis for the downregulated genes indicated that the methionine and cysteine metabolism pathway is specifically downregulated in HCC (Fig. [96]2b). Considering the crucial role of methionine and cysteine metabolism in HCC and the high dependence of HCC on methionine^[97]12,[98]29, our subsequent analysis will focus on this pathway. Exploration of the KEGG database yielded 52 genes related to the methionine and cysteine pathway. By intersecting genesets derived from both methods, 13 Prognostic Methionine and Cysteine Related Differential Expression Genes was identified, namely IL4I1, BCAT1, AGXT2, MAT1A, PSAT1, CTH, SDS, TAT, PHGDH, GNMT, BHMT, CBS, and CDO1, which further served as candidate genes to construct methionine metabolism associated prognostic model. Fig. 2. [99]Fig. 2 [100]Open in a new tab Construction of a Prognostic Model Related to Methionine Metabolism. Volcano plot (a) illustrate differential analysis of tumor and normal samples from TCGA and GTEx, generating 1239 regulated genes and 1506 downregulated genes. Downregulated genes were utilized to perform metabolism associated KEGG pathway enrichment analysis (b). The ratio is defined as the proportion of pathway genes present among the differentially expressed genes relative to the total number of genes in the pathway geneset. With thirteen candidate genes generated by intersection of genes from methionine and cysteine pathway geneset and differentially expressed genes between tumor and normal samples, thirty machine learning algorithms were employed to construct the model, which was ranked based on the average C-index across datasets (c). Forest plot shows the results of univariant (d) and multivariant (e) Cox regression analysis with predict scores of prognostic model based on transcriptomics generated by algorithm ‘RSF + SuperPC’, and clinical indicators. In panel d, clinical indicators include stage, gender, age, and HBV infectious status, while in panel e, clinical indicators include stage and HBV infectious status. The nomogram (f) predicts 1-year, 2-year, and 5-year survival outcomes in TCGA-LIHC. Establishment of the prognostic model Thirty machine learning algorithms were employed to build the model, which was then ranked according to the average C-index across four datasets (Fig. [101]2c, Table [102]S2-S3). The ‘RSF + SuperPC’ algorithm, consisting of 11 genes and achieving an average C-index of 0.586, was determined to have the best prognostic capability. The risk score of primary model = [BCAT1 * 1.9597079] + [PSAT1* 0.9240857] + [CTH * (-1.0413683)] + [PHGDH* 0.8376703] + [TAT * (-1.6079071)] + [GNMT * (-1.2796806)] + [BHMT * (-1.2599247)] + [MAT1A* (-1.0592199)] + [CDO1 * (-1.3322304)] + [SDS * (-0.7058260)]. Expression of each gene are in log[2](x + 1) format. To further develop a model with high prognostic potential, we incorporated clinical information including age, tumor stage, gender, and HBV infection into our analysis. Utilizing datasets from TCGA-LIHC and GSE14520_GPL3921 that encompass these clinical factors, we initially conducted univariate Cox regression analysis with a significance threshold of p < 0.05, which identified tumor stage, HBV infection as having prognostic predictive potential (Fig. [103]2d). Subsequently, multivariant cox analysis was utilized to constructed clinical combined prognostic model, integrating the initial RNA expression-based model scores with the identified high-predictive-value clinical factors (Fig. [104]2e). The nomogram exhibits good prognosis stratification capability (Fig. [105]2f). Model predictions on HCC patient survival Primary and combined models’ efficacy in forecasting survival outcomes in HCC patients were validated using the TCGA-LIHC and [106]GSE14520-[107]GPL571 datasets. HCC patients were split into high and low-risk categories based on the median risk score of the model. KM curve analysis demonstrated that both primary and combined models effectively distinguished patient prognoses in both the training and validation sets. High-risk patients exhibited significantly worse OS outcomes compared to low-risk patients (p < 0.001, Fig. [108]3a-b, Fig [109]S1a-b). The prognostic model constructed using only clinical indicators also demonstrates strong predictive potential (p < 0.001, Fig [110]S2a-b). Fig. 3. [111]Fig. 3 [112]Open in a new tab Functional Validation of the Prognostic Model. Samples were divided into high and low-risk groups based on the median model prediction score, and Kaplan-Meier survival curve analysis was performed for TCGA-LIHC (a) and [113]GSE14520-[114]GPL571 (b). To further evaluate the prognostic value of the model, time-dependent ROC curves were plotted for the model scores in the aforementioned datasets (c-d). To evaluate the prognostic prediction capability of the constructed models over time, Time-dependent ROC curve analysis was administrated in both cohorts to exam constructed model’s prognostic predicting capability in different time nodes. For the primary model, the AUC values at 1, 3, and 5 years in the TCGA-LIHC cohort were 0.624, 0.631, and 0.602, respectively (Fig [115]S1c), while in the validation cohort, these values were 0.701, 0.653, and 0.662, respectively (Fig [116]S1d). In contrast, the combined model exhibited superior performance. in TCGA-LIHC, AUC values for one, three, and five years (Fig. [117]3c) are 0.689, 0.723, 0.715 in sequence. Intriguingly, AUC values in validation cohort (Fig. [118]3d) are successively 0.733, 0.700, 0.719, which indicates perfect potential for model generalization. Moreover, in both datasets, the combined model consistently outperformed the clinical-factors-only model. In TCGA-LIHC, the AUC values for one-, three-, and five-year survival predictions were 0.660, 0.695, and 0.690, respectively (Fig [119]S2c). Similarly, in the validation cohort, the clinical-factors-only model achieved AUC values of 0.677, 0.635, and 0.664 (Fig [120]S2d), all of which were lower than those of the combined model. Reduced expression of BHMT in HCC Earlier research underscored the indispensable function of the gene BHMT in methionine synthesis, a process intimately linked to tumor progression^[121]30. To further investigate function of BHMT in tumor development, we conducted an initial analysis of BHMT mRNA expression in tumor compared to normal tissues. Pan-cancer analysis (Fig. [122]4a) revealed that BHMT is significantly underexpressed in multiple cancers including Breast Cancer (BRCA), Bile Duct Cancer (CHOL), Colon Cancer (COAD), Esophageal Cancer (ESCA), Glioblastoma (GBM), Kidney Chromophobe (KICH), Kidney Clear Cell Carcinoma (KIRC), Kidney Papillary Cell Carcinoma (KIRP), Liver Cancer (LIHC), Rectal Cancer (READ), Sarcoma (SARC), and Stomach Cancer (STAD), while it is significantly overexpressed in Head and Neck Cancer (HNSC), Lung Adenocarcinoma (LUAD), and Lung Squamous Cell Carcinoma (LUSC). Within the TCGA-LIHC dataset, BHMT expression significantly decreased in tumor samples (Fig. [123]4b) and matched tumor tissues (Fig. [124]3c) relative to normal samples. Using a tissue microarray of paraffin-embedded samples from 96 clinical HCC patients, we verified the expression of BHMT protein by immunohistochemical methods. In HCC tissue, there was a notable reduction in BHMT expression contrasted with adjacent non-tumor tissues (Fig. [125]4g-h, Fig S3a, p < 0.001). These results indicate that alterations in methionine metabolism occur in HCC, with BHMT being underexpressed. Additionally, our results suggest that BHMT expression is significantly higher in male HCC tissues compared to female HCC tissues (Fig. [126]4d, p = 0.022). Fig. 4. [127]Fig. 4 [128]Open in a new tab Expression of BHMT in HCC. The ‘TCGAplot’ R package was used to depict the expression levels of BHMT in pan-cancer tissues, contrasting tumor with normal tissues (a). Within TCGA-LIHC, expression levels of BHMT were analyzed in all (b) and paired (c) tumor versus normal samples. Besides, BHMT expression among different Gender (d), Stage (e) and Grade (f) patients were analyzed. Immunohistochemical detection of BHMT expression in cancerous and adjacent normal tissues (g) and box plot comparing BHMT scores between the normal group and the tumor group is presented in penal h. The scale bar represents 100 μm. Low BHMT expression associated with onset and progression of HCC Decreased tendency of BHMT expression was observed with elevated Stage (Fig. [129]4e) and Grade (Fig. [130]4f) of HCC, indicating an association between BHMT expression and the progression of HCC. Given the aberrant expression of BHMT in HCC tumor tissues, survival analysis of BHMT was conducted. KM curves were employed to visually assess the predictive capability of BHMT expression on the prognostic outcomes of HCC patients. Survival curve analysis based on four dataset indicates that patients with low expression of BHMT exhibit worsen overall survival outcomes (Fig. [131]5a-d). Analysis of immunohistochemically assessed BHMT protein expression in relation to clinically relevant pathological information of HCC patients (Table S4) revealed that lower BHMT protein expression is associated with a younger age of onset of HCC (p = 0.0478) and a rising number of solitary tumors (p = 0.0387). Additionally, we utilized immunohistochemical methods to detect the expression of HCC-related antibodies (Ki67, CK7, CK19, Hepatocyte, GPC-3, AFP) in a panel of 96 HCC and adjacent non-tumor tissues (Fig S4). Analyzing the correlation between BHMT and the expression of Ki67, CK7, CK19, Hepatocyte, GPC-3, AFP, we found that BHMT protein expression is related to the expression of the HCC-associated antibody GPC3 (p = 0.0349) (Table S5). Fig. 5. [132]Fig. 5 [133]Open in a new tab The prognostic implication of BHMT in HCC. Kaplan-Meier survival curve analyses highlighted differences in OS between patients with high and low BHMT expression in the TCGA-LIHC (a), GSE14520_GPL3921 (b), [134]GSE116174 (c) and ICGC-LICA-FR (d) databases. BHMT reduction shapes an immunosuppressive microenvironment To further investigate the role of BHMT in HCC, TCGA-LIHC database was utilized to assess whether BHMT expression influences the tumor microenvironment in HCC. Violin plots (Fig. [135]6a) indicate that lower expression of BHMT is primarily associated with reduced infiltration of plasma cells and M1 macrophages, but increased infiltration of memory B cells (Fig. [136]6b), based on deconvolution of 22 immune cell types. The ESTIMATE results suggest that BHMT is overall negatively correlated with stromal cell infiltration and immune cell infiltration (Fig S4a). Additionally, a pan-cancer correlation analysis indicated a general negative correlation between BHMT expression and the expression levels of chemokines (Fig. [137]6c) and chemokine receptors (Fig. [138]6d) in HCC. Analysis of mRNA expression levels of immune checkpoints between high and low BHMT expression groups reveals an overall increase in immune checkpoint expression in the low BHMT group (Fig. [139]6e). Cytokines play a critical role in the recruitment of immune cells; however, the high expression of immune checkpoints inhibits the anti-tumor activity of these immune cells. An analysis using the HCC anti-PD1 immunotherapy dataset [140]GSE202069 indicated that patients with lower BHMT from the model had significantly better outcomes from immunotherapy compared to those with higher expressions (Fig. [141]6f), suggesting BHMT is meaningful for guiding anti-PD1 immunotherapy decisions. In addition, GSEA results revealed that BHMT is negatively correlated with epithelial-mesenchymal transition (EMT), TNFα signaling, and MYC Targets V1 (Fig S4b-d). These findings suggest that downregulation of BHMT may be linked to the formation of an immunosuppressive microenvironment and immune escape in HCC. Fig. 6. [142]Fig. 6 [143]Open in a new tab Low BHMT Expression Induces an Immunosuppressive Tumor Microenvironment. Through CIBERSORT deconvolution of TCGA-LIHC expression profiles, the infiltration levels of immune cells in each sample were quantified. (a) shows the infiltration levels of different cells in groups with high versus low BHMT expression, with white dots indicating median cell infiltration levels. Scatter plots reveal a positive correlation between M1 macrophage infiltration levels and BHMT expression (b). Correlation of BHMT expression at the pan-cancer level with cytokine (c) and cytokine receptor (d) expression. Expression of immune checkpoints in high versus low BHMT expression groups (e). HCC immunotherapy dataset [144]GSE202069 was divided into high and low-BHMT groups based on the median BHMT expression, with Panel F demonstrating significant correlation between immunotherapy outcomes (response and non-response) and BHMT expression (f). *p value < 0.05, **p value < 0.01, ***p value < 0.001. BHMT is associated with fatty acid metabolism reprogramming The top and bottom 20% of HCC patients based on BHMT expression were used for differential analysis, employing the thresholds mentioned in the methods, which yielded 1171 upregulated and 1076 downregulated genes (Fig S5a). Heat maps display the 50 genes most positively (Fig S5b) and negatively (Fig S5c) correlated with BHMT expression. GO (Fig. [145]7a) and KEGG (Fig. [146]7b) pathway analyses suggest that BHMT is involved in various lipid metabolism processes, such as fatty acid metabolism and steroid metabolism, which are essential process involved in tumor reprogramming^[147]31,[148]32 and might be the potential mechanism of BHMT to exert its HCC suppression function. GSEA verified the association results (Fig S5d-e). Increasing evidence suggests that imbalances in fatty acid metabolism are associated with various diseases, including neurodegenerative disorders, metabolic syndrome, and cancer^[149]33, EHHADH and ACSL1 are two key enzymes involved in fatty acid metabolism^[150]34,[151]35. Association analysis demonstrated strong positive association between BHMT expression with EHHADH (Fig. [152]7c, p < 2.2e − 16), as well as with ACSL1 (Fig. [153]7d, p < 2.2e − 16) in HCC. Immunohistochemical assessment of liver tissue microarrays for EHHADH and ACSL1 expression in HCC revealed that both EHHADH (Fig. [154]7e, Fig S3b, p < 0.001) and ACSL1 (Fig. [155]7f, Fig S3c, p < 0.001) were downregulated in HCC tissues. Additionally, by analyzing lipid metabolism and liver injury markers in serum from a cohort of 96 HCC patients, it was found that patients in the high BHMT protein expression group had lower AST levels than those in the low BHMT expression group (p = 0.0112) (Table S6). Nevertheless, the comparison of serum markers, including CHO, HDL, LDL, TG, and APOA1, APOB, did not yield statistically significant differences (all p > 0.05) across the evaluated cohorts (Table S6). These clinical analyses suggest that aberrations in lipid metabolism in HCC could provide a direction for analyzing the mechanisms by which BHMT participates in hepatocyte function. Fig. 7. [156]Fig. 7 [157]Open in a new tab Analysis of Potential Mechanisms by Which BHMT Affects Tumor Prognosis. Utilizing BHMT associated genes, GO (a) and KEGG (b) enrichment analyses were conducted using BHMT-associated genes. In panel A, GO enrichment includes three categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Furthermore, association analysis between BHMT and EHHADH (c), as well as ACSL1 (d) were applied. Based on tissue microarray of paraffin-embedded samples from 96 clinical HCC patients, immunohistochemistry about EHHADH (e) and ACSL1 (f) was performed to investigate their expression in BHMT positive and negative HCC tissues. The scale bar represents 100 μm. BHMT is positive correlation with PPAR signaling pathway and negative correlation with Cytokine-Cytokine receptor interactions Besides, KEGG analysis (Fig. [158]7b) also shows a strong positive correlation between BHMT-related differential genes and the PPAR signaling pathway, a pathway involved in tumor growth arrest, cellular senescence, and clearance, which inhibits tumor progression^[159]36. Additionally, BHMT-related differential genes exhibit a negative correlation with cytokine-cytokine receptor interactions (Fig. [160]7b), which have been proven to play a significant role in remodeling the tumor microenvironment^[161]37. GSEA results verified positive association between BHMT and PPAR signaling pathway, as well as the negative association between BHMT and cytokine-cytokine receptor interaction (Fig S5f-g). Discussion Alterations in the methionine cycle and associated polyamine changes play a fundamental pathogenic role in HCC and are closely linked to poor prognosis in HCC^[162]29. Li et al. found that methionine deprivation induced irreversible cell cycle arrest in hepatocellular carcinoma cells, confirming the potential applicability of methionine-restricted dietary therapies in HCC^[163]38. Rosa M. Pascale and colleagues highlighted how disruptions in methionine metabolism promote HCC through reduced expression of SAM. SAM acts as a tumor suppressor by providing antioxidant effects, inhibiting NF-κB activation, and inducing the overexpression of tumor suppressor gene PP2A^[164]29. Another notable study emphasized the high dependence of tumor-initiating cells (TICs) on methionine, suggesting that transient pharmacological inhibition of the methionine cycle can significantly reduce the tumor-initiating capacity of TICs, highlighting methionine metabolism as a potential therapeutic target^[165]39. In this study, our enrichment analysis of metabolic pathways among downregulated genes in HCC identified significant enrichment of the methionine and cysteine metabolism pathway. Building on this finding—and considering the established role of methionine metabolism in HCC progression—we focused our investigation on this pathway, emphasizing its potential for effective patient stratification and prognostic prediction. In recent years, with the increasing research on the pivotal role of the liver in metabolism and tumor immune reprogramming, metabolic-related prognostic models for HCC have emerged rapidly, focusing on critical processes such as arachidonic acid metabolism, cholesterol metabolism, and selenium metabolism^[166]40–[167]42. Furthermore, HCC has been shown to exhibit dependence on exogenous methionine, which further supports the metabolic reprogramming of methionine metabolism and its importance in HCC^[168]12. In this work, utilizing multiple databases and through the selection of optimal machine learning algorithms, we identified a prognostic model related to methionine metabolism involving genes such as DNMT3B, DNMT3A, GNMT, DNMT1, CTH, BCAT1, MAT1A, BHMT, CDO1 and TAT. Among these, BCAT1, PSAT1, and PHGDH were identified as risk factors for HCC, whereas CTH, TAT, GNMT, BHMT, MAT1A, CDO1, and SDS were identified as protective factors. However, this model, entirely based on transcriptomic information, captures transcriptional differences among patients but does not account for clinical factors that have been proven to significantly predict prognosis. To further enhance the model’s predictive accuracy and fully integrate clinical indicators, we incorporated tumor stage and HBV infection status using multivariate Cox regression analysis to construct a clinical feature-enhanced methionine metabolism-associated model. This represents the first prognostic model based on methionine metabolism, addressing a significant gap in the field. The model underscores the critical role of methionine metabolism pathways in HCC progression, laying the groundwork for future research into methionine metabolism and HCC therapeutic strategies. Notably, in our model, HBV functions as a protective factor against HCC mortality. Existing research primarily focuses on the mechanisms by which HBV contributes to HCC development, including the induction of genomic instability, chronic inflammation, and alterations in the tumor microenvironment^[169]43,[170]44. However, our findings suggest that HBV may play a more diverse role in the progression of HCC, warranting further investigation. Constructed model exhibits excellent prognostic capability among various examinations, with a suitable fit for generalization. The primary pathway for methionine synthesis in the liver is the betaine synthesis pathway^[171]45, while research has reported that the liver-specific distribution of BHMT determines the confinement of the betaine synthesis pathway to the liver, highlighting BHMT as a key enzyme in this pathway^[172]46. We hypothesize that dysregulated BHMT expression may induce tumor exogenous methionine dependence. Given that previous studies focused on the prognostic impact of BHMT SNPs^[173]19,[174]47–[175]49, we expanded the analysis of BHMT at the RNA level. Analysis of pan-cancer data indicated significant downregulation of BHMT in various cancers, including BRCA, CHOL, COAD, ESCA, GBM, KICH, KIRC, KIRP, LIHC, READ, SARC, and STAD, while significantly upregulated in HNSC, LUAD, and LUSC. TCGA-LIHC data analysis and IHC results also indicated significant downregulation of BHMT in HCC, with lower BHMT expression correlating with HCC prognosis and worsen patient outcome, consistent with previous studies^[176]22. Interestingly, we found that low expression levels of BHMT are associated with an earlier onset age and prognosis of HCC, providing guidance for early screening in high-risk populations. BHMT expression is significantly upregulated in male HCC tissues, suggesting that BHMT may be regulated by gonadal hormones and may perform distinct functions between sexes. This gap warrants further investigation. Additionally, BHMT expression positively correlates with GPC3, a biomarker used for diagnosing HCC and predicting poor prognosis. Our results suggest that BHMT holds significant value in the diagnosis and prognosis of HCC and is intricately linked to disease progression. The development of solid tumors depends on intrinsic signals within cancer cells and external factors from the tumor microenvironment, where an immunosuppressive environment promotes tumor proliferation, dissemination, and immune evasion^[177]50. Nevertheless, no literature has reported the regulatory role of BHMT in the HCC immune microenvironment. In this work, we found BHMT expression positively correlates with plasma cell and M1 macrophage infiltration, while generally negatively correlating with memory B cells and stromal cells infiltration, chemokines and receptors, suggesting a novel link between BHMT and tumor microenvironment remodeling. Moreover, in samples with low BHMT expression, we observed a widespread high expression of immune checkpoints, which further impairs the anti-tumor activity of immune cells, leading to the accumulation of non-cytotoxic immune cells within the tumor. Further analysis revealed that BHMT is an effective predictor of response to anti-PD-L1 therapy, with the potential to benefit a wide range of patients. Interestingly, our results indicate that BHMT is negatively correlated with the EMT, the TNFα signaling pathway, and the MYC targets V1 pathway, all of which are known to promote tumor progression as documented in previous studies^[178]51–[179]53. This suggests that the tumor-suppressive role of BHMT may be partially mediated through its negative regulation of these pathways. Studies have reported that memory B cells within tumors exhibit high heterogeneity, with different subsets playing either anti-tumor or pro-tumor roles^[180]54. In our study, we observed a negative correlation between BHMT expression and memory B cell infiltration. Given the diversity of memory B cell infiltration, further subclassification of memory B cells to identify the specific subsets influenced by BHMT is essential, providing direction for future research. Similarly, existing studies have highlighted the diverse roles of stromal cells in shaping the tumor microenvironment^[181]55, underscoring the need for further research to identify the specific cellular subpopulations influenced by BHMT. Plasma cells actively communicate with other immune cells, enhancing the cytotoxic effects of CD8^+ T cells and NK cells in HCC^[182]56,[183]57. Besides, studies have reported that M1 polarized macrophages activate anti-tumor immune responses by directly phagocytosing tumor cell and presenting associated antigen, secreting multiple pro-inflammatory cytokines^[184]58. This potential linkage warrants further exploration using advanced technologies such as spatial transcriptomics and single-cell techniques. Hepatic fat accumulation induces steatohepatitis, which can subsequently lead to the development of HCC. Chul et al. reported that betaine, involved in the methionine cycle, can reduce homocysteine levels by upregulating BHMT, thereby inhibiting hepatic lipid accumulation and reducing HCC incidence^[185]59. Teng YW proposed that BHMT is involved in the regulation of lipid synthesis, suggesting a potential interaction between BHMT and lipid metabolism^[186]60. However, no studies have specifically pinpointed the processes affected by BHMT. Through KEGG and GO enrichment analyses, we found that BHMT is closely related to fatty acid metabolic pathway. Further analysis revealed a strong positive correlation between BHMT and the gene expression of two key enzymes in fatty acid metabolism, EHHADH and ACSL1. EHHADH, a member of the 3-hydroxyacyl-CoA dehydrogenase family, plays a critical role in the peroxisomal fatty acid β-oxidation pathway in hepatocytes. Previous studies have demonstrated that alterations in fatty acid synthesis, β-oxidation, and cellular lipid composition contribute to the initiation and progression of hepatocellular carcinoma (HCC)^[187]61. ACSL1, a member of the long-chain acyl-CoA synthetase (ACSL) family, plays a pivotal role in fatty acid metabolism by catalyzing the activation of fatty acids for energy production or their conversion into phospholipids, cholesterol esters, and triglycerides. Studies have shown that both EHHADH and ACSL1 regulate fatty acid metabolism through PPAR family member-dependent signaling pathways^[188]62,[189]63. Serum lipid metabolites and liver function indicators showed that the AST levels were significantly reduced in the group with low BHMT expression, whereas ALT levels did not exhibit a significant difference. This may be related to the increased activity of anti-apoptotic genes or the downregulation of pro-apoptotic genes^[190]64. Fatty acids serve as substrates for various major lipids in the body, and their metabolism is a crucial downstream process in lipid metabolism. It has been reported that fatty acid metabolism maintains lipid metabolic balance and regulates the overall lipid metabolic rate through various mechanisms^[191]65,[192]66. Our study provides direction for future studies to uncover the specific sites where BHMT influences these pathways. Further exploration of the mechanisms and potential functions of BHMT in inhibiting HCC development revealed also through enrichment analysis. Interestingly, BHMT expression positively correlates with PPAR signaling pathway and negatively correlates with cytokine-cytokine receptor interaction, in addition to its association with fatty acid metabolism. PPARs are members of the nuclear hormone receptor superfamily, known to induce growth arrest and apoptosis in various malignancies^[193]67, while cytokine-cytokine receptor interactions play an irreplaceable role in shaping the tumor microenvironment^[194]68. These pathways may represent the potential mechanisms through which BHMT exerts its tumor-suppressive function, thereby expanding the understanding of BHMT’s role in HCC, providing foundational basis for future research into its multifaceted functions in the context of HCC. In this study, we constructed the first prognostic model related to methionine metabolism, which demonstrated superior performance. Additionally, transcriptome analysis of BHMT indicated its downregulation in HCC and association with HCC prognosis and poor prognosis. CIBERSORT was employed to estimate immune cell infiltration, revealing that samples with high BHMT expression had increased infiltration of plasma cells and M1 macrophages. BHMT expression is negatively associated with multiple immune checkpoints, and we demonstrate that BHMT serves as an excellent biomarker for predicting the outcomes of anti-PD-L1 immunotherapy. Combining enrichment results with experimental validation, we show BHMT is associated with fatty acid metabolism reprogramming, providing potential targets in this pathway. Besides, enrichment analysis also found that BHMT positively correlates with PPAR signaling pathway and negatively correlates with cytokine-cytokine receptor interaction. These suggest potential mechanisms by which BHMT exerts its tumor-suppressive functions in HCC. Electronic supplementary material Below is the link to the electronic supplementary material. [195]Supplementary Material 1^ (865KB, xlsx) [196]Supplementary Material 2^ (18.4MB, docx) Acknowledgements