Glutamine Metabolism Scoring Predicts Prognosis and Therapeutic Resistance in Hepatocellular Carcinoma

Glutamine metabolism (GM) plays a critical role in hepatocellular carcinoma (HCC); however, a comprehensive methodology to quantify GM activity is still lacking. We developed a transcriptome-based GMScore to evaluate GM activity and investigated the association of GMScore with prognosis and therapeutic resistance. Two independent HCC cohorts with transcriptome data were selected from The Cancer Genome Atlas (TCGA, n = 365) and the International Cancer Genome Consortium (ICGC, n = 231). The expression of 41 GM-associated genes were used to construct and validate GMScore. Several genomic or transcriptomic biomarkers were also estimated. Tumor response to immune checkpoint inhibitors (ICIs) was predicted using the tumor immune dysfunction and exclusion algorithm. GMScore was closely correlated with patient characteristics, including stage, histology grade, alpha-fetoprotein level, and vascular invasion. High GMScore was an independent risk factor for overall survival (OS) in both cohorts (HR = 4.2 and 3.91, both p < 0.001), superior to clinical indices and other biomarkers. High GMScore presented transcriptome features to indicate cell growth advantages and genetic stability, which was associated with poor OS of patients who received transcatheter arterial chemoembolization (TACE). High GMScore was also related to high expression of immune checkpoint genes, increased infiltration of regulatory T cells, and decreased infiltration of M1 macrophages. More importantly, high GMScore indicated poor predicted responses to ICIs, which could be verified in an ICI-treated melanoma cohort. In conclusion, GMScore is a strong prognostic index that may be integrated into existing clinical algorithms. A high GMScore may indicate resistance to TACE and ICIs based on its transcriptome and immune features. Validations using other HCC cohorts, especially ICI-treated HCC cohorts, are necessary.


INTRODUCTION
Primary liver cancer is the seventh most prevalent cancer and the second leading cause of cancer-related death in the world [1]. Hepatocellular carcinoma (HCC) is the most common form of liver cancer, accounting for approximately 90% of cases. Currently, hepatitis B and hepatitis C virus infections remain the most important global risk factors for HCC. Meanwhile, metabolic risk factors, such as obesity, metabolic syndrome, diabetes mellitus, and non-alcoholic fatty liver disease are increasingly prevalent and have become important causes of HCC [2]. However, regardless of pathogenic factors, abnormal metabolism plays a critical role in hepatocellular carcinogenesis and HCC development because the liver is a metabolic organ [3].
Aberrant active glutamine metabolism (GM) is a key player in HCC. The tumorigenicity of HCC stem cells was inhibited by knocking out glutaminase 1 (GLS1), with a high expression of GLS1 predicting poor prognosis [4]. Moreover, the expression of glutamine transporter ASCT2 was significantly upregulated in HCC and was an independent prognostic risk factor [5]. In contrast, high expression of oxoglutarate dehydrogenase-like, which limited GM, was associated with favorable prognosis of HCC patients and sensitized HCC cells to sorafenib [6]. These findings indicate that the expression levels of GM-related genes reflect the molecular heterogeneity of HCC, which determines distinct clinical outcomes of HCC patients. However, there is still a lack of models to quantify GM activity based on relevant gene expressions whose acquisition has been largely promoted by nextgeneration sequencing (NGS).
In this study, we established a transcriptome-based methodology named GMScore to quantify GM activity in HCC. We showed that GMScore was a predictor not only for prognosis but also for treatment outcomes in two independent HCC cohorts.

Patients and Data Collection
We screened patients from The Cancer Genome Atlas (TCGA) as the discovery set and patients from the International Cancer Genome Consortium (ICGC) as the validation set. The following enrollment criteria were used: 1) available sequencing data for GMScore calculation, 2) pathological diagnosis of HCC, and 3) no prior history of radiation therapy, chemotherapy, target therapy, immunotherapy, or other anticancer medications (including neoadjuvant therapy). Sequencing data and corresponding clinical information, recorded up to August 1, 2021, were downloaded from TCGA (https://portal.gdc.cancer.gov/ repository) and ICGC (https://dcc.icgc.org/projects/LIRI-JP) data portals. An additional melanoma cohort (GSE78220) treated with anti-PD-1 was used to verify the association of GMScore with immunotherapy outcomes, and data were downloaded from Gene Expression Omnibus (GEO) [7]. Genetic expression data presented as the fragments per kilobase per million (FPKM) values were transformed into transcripts per kilobase million (TPM) values to improve comparability between samples [8]. The American Joint Committee on Cancer criteria was used for clinical and clinicopathological classification and staging.

GMScore Construction
The 41 GM-related genes (Supplementary Table S1) were extracted from the Gene Ontology (GO) initiative and a published study [9]. Of them, differentially expressed genes between cancerous and the adjacent normal tissues were identified by the limma R package in HCC patients of TCGA, with a false discovery rate of <0.05. Subsequently, genes associated with overall survival (OS) were selected using univariate Cox regression models. The optimal cut-off values to define high and low expression subgroups were determined based on the association of genetic expression with OS using the Survminer R package. Subsequently, gene expression level was evaluated as 0 or 1; a value of 0 was assigned when the gene expression was less than the corresponding cut-off value, and a value of 1 otherwise. Furthermore, the least absolute shrinkage and selection operator (LASSO) Cox regression model was used to screen the most useful prognostic genes. A GMScore model was then constructed based on the fraction of selected genes using Cox regression coefficients. The formula was established as follows: GMScore sum (each gene's expression × corresponding coefficient).

Immune, Stromal, and ESTIMATE Scoring
The ESTIMATE algorithm was used to estimate the relative fraction of stromal and immune cells in the tumor microenvironment (TME) and was exhibited in the form of the StromalScore, ImmuneScore, and ESTIMATEScore. Of these, ESTIMATEScore correlated with DNA copy numberbased tumor purity [10].

Tumor Mutation Burden Calculation
The total counts of somatic nonsynonymous variations, including missense, nonsense, splice-site, and frameshift mutations, in coding regions were defined as TMB [11].

Gene Set Enrichment Analysis
Differentially expressed genes between the low and high GMScore subgroups were determined, with the criteria of adjusted p-value <0.05 and log2 (fold change) >1, and GSEA based on GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using NetworkAnalyst 3.0 [12].

Immune Infiltration Estimation
The CIBERSORT algorithm and the LM22 gene signature were utilized to quantify the abundance of infiltrating immune cells in the TME, based on transcriptome data [13].

Tumor Immune Dysfunction and Exclusion Scoring
TIDE, a computational method based on transcriptome to model the induction of T cell dysfunction and the prevention of T cell infiltration in tumor immune evasion, was used to calculate T cell dysfunction and exclusion scores, which were further merged as the TIDE score to predict tumor response to immune checkpoint inhibitors (ICIs) [14].

Statistical Analysis
R software (version 4.0.3, http://www.r-project.org) or IBM SPSS Statistics ver.20 (IBM Corp., Armonk, NY, United States) were used to analyze the related data and plot graphs. The Student's t-test, chi-square test, Fisher's exact probability test, or Mann-Whitney U test were used to compare the differences between groups. The Kaplan-Meier method with the log-rank test was used to compare OS between different parts. Hazard ratios (HRs) and their 95% confidence intervals (CIs) for prognostic factors were calculated using univariate and multivariate Cox proportional hazard models. Patients were separated into high and low GMScore groups based on the optimal value associated with OS calculated using the Survminer R package. Receiver operating characteristic curve (ROC), time-dependent ROC curves, and the areas under the ROC curves (AUC) depicted or calculated by the timeROC R package, were used to evaluate the predictive power of the models. Statistical significance was set at p < 0.05, and all p-values were two-tailed.

Patient Characteristics
A total of 365 patients from TCGA and 231 patients from ICGC were included (Supplementary Table S2). Patients in ICGC were older and had more stage III/IV diseases than those in TCGA (p < 0.05). Histology grade, alpha-fetoprotein (AFP), and vascular invasion data were only available in TCGA.

Derivation of the GMScore and its Association With Clinical Features
A total of 30 GM-related genes were differentially expressed between HCC and normal tissues ( Figure 1A). Of these, 20 genes were significantly associated with OS ( Figure 1B Patients with high GMScores in both TCGA and ICGC included a significantly higher proportion of patients with stage III/IV disease (p < 0.001 and p 0.009, respectively) than those with low GMScores. Histology grade 3/4 (p < 0.001), AFP larger than 200 ng/ml (p < 0.001), and vascular invasion (p 0.001) were also more frequent in patients with high than in those with low GMScores ( Table 1).

GMScore is an Independent Prognostic Factor
Patients were divided into high and low GMScore subgroups to establish a prognosis-predictive model in TCGA, which was validated using ICGC ( Figure 2A). Kaplan-Meier survival analysis confirmed the survival discrepancy between high and low GMScore groups ( Figure 2B), and the results of the ROC curve analysis verified the predictive value of the established risk model in both TCGA and ICGC ( Figure 2C). After univariate selection for prognostic significance of variables (Supplementary Figure S1), multivariate analysis in both TCGA and ICGC showed that high GMScore was an independent predictor of poor OS (HR 4.2, 95% CI 2.38-7.4, p < 0.001, and HR 3.91, 95% CI 1.92-7.97, p < 0.001, respectively), even superior to staging in terms of HR ( Figure 2D).

GMScore is Superior Than Other Biomarkers for Prognostic Prediction
ImmuneScore, StromalScore, ESTIMATEScore, and TMB have been reported as prognostic predictors in many cancers [10,11,15]. We showed that they were also associated with or tended to be associated with OS of HCC (Supplementary Figure S2). However, prognostic stratification according to these factors could be further optimized using GMScore. The OS of both patients in the high and low subgroups of these biomarker values could be further stratified by the GMScore level ( Figures 3A,B). Moreover, the ROC curve analysis revealed that GMScore was better than other biomarkers in terms of the prediction of 3-years OS (only two patients were followed up for over 5 years in ICGC) ( Figures 3C,D).

Transcriptome Features of High GMScore Indicate Therapeutic Resistance
Differentially expressed genes between the high and low GMScore subgroups were determined ( Figure 4A and Supplementary Table S3). GSEA based on the GO biological  Table S4). We revealed that genes associated with proliferation and cell cycle were significantly enriched in the high GMScore subgroup, indicating that tumors with high GMScores might have cell growth advantages ( Figure 4B). Meanwhile, high GMScore also enriched genes associated with DNA repair and chromatin regulation, suggesting that tumors with high GMScores may be characterized by genomic stability ( Figure 4C). It is known that genomically stable or elevated DNA repair correlates with therapeutic resistance; thus, we examined the impact of high GMScore on outcomes of transcatheter arterial chemoembolization (TACE). In the TCGA HCC cohort, high GMScore significantly decreased OS in patients who received adjuvant TACE after surgery (p 0.003; Figure 4D) or salvage TACE after recurrence (p < 0.001; Figure 4E) than those with low GMScore. A similar result was found in the ICGC HCC cohort ( Figure 4F), although the timing of TACE was unknown.

GMScore is Associated With Immune Profiles and May Predict Response to ICIs
TMB has been identified as a predictor of ICI efficacy in multiple cancers. In this study, we found that GMScore level did not affect TMB ( Figure 5A). However, high GMScore subgroup had significantly higher gene expressions than low GMScore subgroup in terms of several immune checkpoints, including PD-1, CTLA-4, TIM-3, and TIGIT in both TCGA and ICGC ( Figure 5B). For immune cell infiltration in tumors, high GMScore subgroup had a significantly higher abundance of regulatory T cells (Tregs) but a significantly lower abundance of M1 macrophages than low GMScore subgroup ( Figure 5C). These findings indicate that GMScore may partly reflect TME.
Because no HCC cohorts treated by ICIs with transcriptome data were available, we used the TIDE algorithm to calculate T cell dysfunction and exclusion scores, and predict ICI response by the merged TIDE score. We observed that high GMScore subgroup had significantly higher exclusion score and TIDE score compared to the low GMScore subgroup ( Figure 5D). Consequently, for the predicted response rate, high GMScore HCC was significantly inferior to low GMScore HCC in both TCGA and ICGC ( Figure 5E). Furthermore, we used an ICItreated melanoma cohort to verify the impact of GMScore on immunotherapy outcomes and found that the high GMScore subgroup had a significantly lower objective response rate than the low GMscore subgroup (14.3 vs. 66.7%, p 0.016; Figure 5F). More importantly, high GMScore significantly reduced OS compared to low GMScore (p 0.039; Figure 5G).

DISCUSSION
This is the first study to evaluate GM activity by quantifying GMassociated gene expression as the GMScore in HCC. We revealed that GMScore was an independent OS predictor superior to clinical indices and other biomarkers. HCC with high GCS scores may be genomically stable and resistant to TACE. In addition, GMScore was associated with immune profiles and may be helpful in the prediction of response to ICIs. Together, the methodology of GM evaluation, such as ours, may improve the understanding and clinical outcomes of HCC.
In the recent decade, improvements in therapeutic approaches have been achieved in HCC. In contrast, there has not been any significant clinical improvement in HCC biomarkers in predicting prognosis and treatment response. Some existing clinical algorithms that have been proposed for HCC do not include biomarkers, such as the Barcelona Clinic Liver Cancer (BCLC) system, the Japan Integrated Scoring (JIS) system, and the Hong Kong Liver Cancer (HKLC) system. On the other hand, other existing clinical algorithms, such as the albumin-bilirubin (ALBI) grade and the BALAD score, include only traditional biochemical parameters such as AFP, serum bilirubin, and albumin [16]. Novel biomarkers, especially those developed by NGS, such as genomic and transcriptome features, are still far from clinical application, although increasing evidence indicates their critical impacts on clinical outcomes of HCC [17,18]. In this study, GMScore had a stronger prognostic effect on HCC than staging, which was independent of clinical indices. GMScore also has potential influences on treatment outcomes. Further studies are needed to test whether GMScore can be integrated with existing clinical algorithms to further improve clinical practice.
GMScore allowed us to explore GM-associated transcriptome features. It is reasonable to observe an association between high GMScore and cell growth advantages indicated by GSEA, because glutamine is an essential nutrient for cancer cell proliferation. Interestingly, high GMScore may improve DNA repair and genomic stability, which may induce TACE-resistance in both the TCGA and ICGC HCC cohorts. Many pieces of evidence support our findings. Glutamine deficiency induces DNA alkylation damage and DNA damage accumulation to trigger genomic instability and hyposensitize cancer cells to alkylating agents [19]. Glutamine synthetase (GS), an enzyme catalyzing glutamate and ammonia to glutamine, improves DNA repair and causes radiation resistance [20]. Moreover, high GS expression decreased OS and increased early phase recurrence in 554 HCC patients from two independent cohorts who underwent adjuvant TACE [21].
The deleterious role of GM in anticancer immunity is intriguing. Glutamine blockade suppressed cancer cells but induced a long-lived, highly activated phenotype for effector T cells [22]. Targeting GM rendered ICI-resistant tumors susceptible to immunotherapy by modulating myeloid-derived suppressor cells (MDSCs) in a breast cancer model [23]. In combination with anti-PD-L1, glutamine depletion in mice strongly promoted antitumor efficacy of T cells by increasing Fas/CD95 levels [24]. Our study revealed that high GMScore may promote the genetic expression of immune checkpoints and the infiltration of immunosuppressive Tregs, but impede the infiltration of immunoactivated macrophages M1. Similarly, a previous study on kidney cancer also identified a glutamine signature (GlnS) to show that high-GlnS tumors had higher levels of Tregs and impaired T-cell cytotoxic function [9]. However, the impact of GM on immunotherapy efficacy has not been investigated in cancer patients. In this study, we predicted the tumor response to ICIs in patients with HCC using the well-established TIDE algorithm. We showed that GMScore may be associated with such a response, which could be validated in an ICI-treated melanoma cohort. These findings further support the development of GM inhibitors in clinical trials and GM evaluation in the efficacy prediction of immunotherapy.
This study has several limitations. First, more HCC patients are necessary to validate the accuracy of prognosis prediction based on GMScore. Second, the findings of this study have not been verified experimentally. Moreover, predictions rather than true tumor responses to ICIs in HCC patients were used, although the impact of GMScore on immunotherapy was confirmed in melanoma.
In conclusion, we showed that GMScore, a multigene model related to GM, could predict both prognosis and therapeutic resistance in HCC. In particular, GMScore may impact immunotherapy response, which is highly interesting considering the recent arrival of the immunotherapy era in HCC. These results will help in understanding of the cancer biology, stratifying the prognosis, developing precise medical strategies, and improving the survival of HCC patients. Our study stressed that GM evaluation and GM blockade could be an interesting area for further investigation.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: HCC data in our study is available in TCGA (https://portal.gdc.cancer.gov/repository) and ICGC (https://dcc. icgc.org/projects/LIRI-JP) data portals. The data of ICI-treated melanoma cohort is available in Gene Expression Omnibus under accession number GSE78220 (https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc GSE78220).

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.