Bulk RNA-seq data from a prior Renji Hospital study (14), including paired primary tumor and adjacent normal tissues from three FH-deficient RCC patients, were analyzed to assess differential gene expression patterns and their functional roles in glycolysis. Additionally, another FH-deficient RCC bulk RNA-seq cohort [GSE157256 dataset (15)] was adopted to investigate the different expression patterns across metastasis, primary tumors, and adjacent normal tissues.

To explore the associations between glycolysis and ICI response, four scRNA-seq datasets were evaluated: an FH-deficient RCC cohort [from a published article in the journal Clin Cancer Res ( 16 ).], a clear cell renal carcinoma (ccRCC) cohort [SRP308561 (17)], a skin cutaneous melanoma (SKCM) cohort [GSE115978 (18)], and a basal cell carcinoma (BCC) cohort (GSE123813 (19)). The obtained gene expression matrices were converted into Seurat objects, and all subsequent analyses were conducted using R software. Probable doublets were first removed using the DoubletFinder package. After integrating Seurat objects across all samples, strict quality control (QC) filters were applied to exclude cells meeting any of the following criteria: genes detected >7,500 or <200, total reads >75,000, or mitochondrial RNA content >20%. Qualifying cells underwent normalization and scaling via the LogNormalize method in Seurat, which concurrently identified highly variable genes. Principal component analysis (PCA) was performed on the highly variable gene matrix, with significant principal components (PCs) selected based on elbow plot inflection points and heatmap-driven evaluation of variance contributions. Using these PCs, a k-nearest neighbor graph was constructed (FindNeighbors), and graph-based clustering (FindClusters) was executed at a resolution parameter of 1. Cluster identities were visualized via UMAP projection, and cell annotations for all clusters were provided by the literature from the data sources.

Glycolysis-related genes from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) were used to compute glycolysis pathway enrichment scores via gene set variation analysis (GSVA). Cohort details are summarized in Supplementary Table S1 .

A pan-cancer glycolysis signature (Glyc.Sig) was developed using 41 scRNA-seq datasets (405 patients, 1,220,365 cells, 21 cancer types, Supplementary Table S2 ) encompassing malignant, stromal, and immune cells from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and TISCH portal (http://tisch.comp-genomics.org/) (20), as well as directly from the published articles. These included breast cancer (BRCA), basal cell carcinoma (BCC), clear cell renal cell carcinoma (ccRCC), colorectal cancer (CRC), cholangiocarcinoma (CHOL), glioma, fumarate hydratase-deficient renal cell carcinoma (FH-deficient RCC), head and neck cancer (HNSC), liver hepatocellular carcinoma (LIHC), multiple myeloma (MM), Merkel cell carcinoma (MCC), medulloblastoma (MB), non-small cell lung cancer (NSCLC), neuroendocrine tumor (NET), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), and uveal melanoma (UVM) (16–19, 21–53).

TCGA pan-cancer transcriptomic data (10,154 patients; 30 cancer types) from the UCSC Xena data portal (https://xenabrowser.net) (54) were analyzed to evaluate the associations between Glyc.Sig and immune suppression. Acute myeloid leukemia (AML), diffuse large B-cell lymphoma (DLBC), and thymoma (THYM) were excluded due to immunocyte predominance (55). Tumor mutation burden (TMB) data from cBioPortal (https://www.cbioportal.org) (56, 57) and intratumor heterogeneity (ITH) data from a published article [Thorsson et al. (58)] were introduced to assess their correlations with Glyc.Sig scores calculated by GSVA.

Pretreatment transcriptomic data with clinical information from 10 ICI RNA-seq cohorts were used to develop and validate a Glyc.Sig-driven predictive model. These include five SKCM cohorts, provided by Van Allen in 2015 (59), Riaz in 2017 (60), Hugo in 2016 (61), Liu in 2019 (62), and Gide in 2019 (63); one renal cell carcinoma (RCC) cohort provided by Braun in 2020 (64); two urothelial carcinoma (UC) cohort from Mariathasan in 2018 (65) and Synder in 2017 (66); one glioblastoma multiforme (GBM) cohort from Zhao in 2019 (67); and one gastric cancer (GC) cohort from Kim in 2018 (68). Detailed cohort characteristics are demonstrated in Supplementary Table S3 .

Seven CRISPR/Cas9 screening datasets derived from previous studies by Freeman (69), Kearney (70), Manguso (71), Pan (72), Patel (73), Vredevoogd (74), and Lawson (75) across multiple cancer types (BRCA, CRC, RCC, and SKCM) were analyzed to identify immune resistance-related genes. Following Fu et al. (76), who curated the first six datasets (except Lawson’s cohort), we divided the data from these 7 datasets into 17 distinct groups ( Supplementary Table S4 ). The CRISPR screening methodological framework involved genome-wide CRISPR-Cas9 knockout in cancer cell lines subjected to cytotoxic lymphocyte (CTL) co-culture/not subjected to CTL (in vitro), or xenograft models in immune-deficient or immune-competent (in vivo) mice, followed by sgRNA abundance quantification through RNA sequencing (RNA-seq). Immunomodulating effects were quantified by calculating log-fold changes in sgRNA reads between experimental pairs (CTL-treated vs. untreated; immune-deficient vs. immune-competent mice) (75). Subsequent Z-score normalization was conducted for cross-dataset comparisons. Gene rankings were determined through an average Z-score across all datasets, and top-ranked genes with low Z-scores were regarded as potential mediators of immune resistance.

GSVA (R package “GSVA”) was employed to calculate scores of HALLMARK pathways, Glyc.Sig, and glycolysis intensity (using glycolysis-related genes obtained from MsigDB). Pathway enrichment analysis was conducted using data from the Reactome Knowledgebase (https://reactome.org) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/) via the R package “clusterProfiler.” Immune infiltration was quantified using the R package MCP-counter v1.1.0.

Immunohistochemical (IHC) staining was carried out on five tumor samples with adjacent normal tissues from FH-deficient RCC patients at Renji Hospital, with ethical committee approval. The experimental procedure included the following steps:

The final IHC score was determined by multiplying the scores for the percentage positivity of target protein-expressing cells.

FH-deficient RCC patients’ tumors and adjacent normal tissues were lysed in 2% SDS, followed by thermal denaturation at 99°C for 30 min. Proteins were separated by SDS-PAGE and then transferred onto nitrocellulose membranes. After blocking with 3% BSA in TBST for 1 h at room temperature, membranes were probed with primary antibodies through overnight incubation at 4°C. Following three washes with TBST, membranes were incubated with species-matched HRP-conjugated secondary antibodies (1:5,000 dilution) for 1 h at ambient temperature. Protein bands were visualized using chemiluminescence (Thermo Fisher Scientific). Paired samples from each patient were always run on the same gel to ensure comparability. Antibodies used in this study included anti‐LDHA (19987-1-AP, Proteintech) and anti‐β‐actin (66009-1-Ig, Proteintech).

Analyses were performed in R v4.3.1 (https://www.r-project.org). A two-sided Wilcoxon test was adopted to compare glycolysis scores between the ICI response and non-response groups. Spearman correlation analysis assessed associations between Glyc.Sig and other biological signatures, including scores of HALLMARK pathways, immune-related genes, ITH, TMB, and immune infiltration. FDR was adjusted via Benjamini–Hochberg. Model training, validation, and testing were conducted based on the R package cancerclass and caret. The predictive performance of the models was evaluated by ROC and AUC (86). Survival differences were analyzed using Cox regression analysis.

FH-deficient RCC is an aggressive cancer syndrome driven by inactivation of the fumarate hydratase gene. To investigate the potential phenotype correlated with FH-deficient RCC, we utilized data from the Renji cohort in our prior investigation (14), which contains bulk RNA-seq data of primary cancer and adjacent normal tissues from three FH-deficient RCC patients ( Figure 2A ). After that, we compared the expression levels of different genes between cancer and adjacent normal tissues ( Figure 2B ), and some glycolysis-related genes, such as ENO1, PLOD2, and PKM, were upregulated. Considering loss-of-function mutations in FH trigger a metabolic crisis characterized by defective tricarboxylic acid (TCA) cycle flux, forcing cellular metabolic reprogramming toward aerobic glycolysis, we speculated that glycolysis might be more active in FH-deficient RCC cancer tissues. Therefore, we conducted the GSEA analysis of glycolysis pathways, and predictably, there was a significant positive enrichment of glycolysis pathways in cancer tissues of FH-deficient RCC ( Figure 2C ). Next, we investigated the relationship between glycolysis and ICI outcomes in FH-deficient RCC. Figure 2D showed the t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization of the FH-deficient RCC from previously published datasets. Patients with progressive disease (PD) exhibited significantly higher levels of glycolysis compared to those with partial response (PR) ( Figure 2E ). In particular, glycolysis was more enriched in cancer cells (epithelial cells in Figure 2F ) compared to other cell types. As glycolysis is a key metabolic pathway that is often upregulated in many cancers (87), we examined the findings above in another renal cancer type, clear cell renal cell carcinoma (ccRCC), and similar findings were observed, as shown in Figures 2G–I , further supporting the potential role of glycolysis in influencing the outcomes of ICI therapy.

To further investigate the role of glycolysis in ICI responses, the relationships between ICI outcomes and glycolysis intensity were explored in BCC and SKCM cancer datasets. In BCC, Figures 3A, B demonstrate that non-responders (NR) exhibited a higher glycolysis intensity compared to responders (R). A quantitative analysis shown in Figure 3C further confirms this finding, revealing that the NR subgroup had significantly higher glycolysis intensity than the R subgroup. For SKCM, an additional dataset was used (18), which excluded samples lacking data on malignant cells, to validate the findings. However, due to missing data for some responders, the comparison was made between treatment-naive (TN) patients and NR. It was hypothesized that treatment-naive patients might consist of both potential responders and non-responders. The analysis also demonstrated a significantly higher glycolysis intensity in the NR subgroup than the TN subgroup in the SKCM cohort, as shown in Figures 3D–F (p < 0.001).

Since the intensity of glycolysis was significantly correlated with ICI resistance, we suggested that a gene set containing genes reflecting the level of glycolysis, named as glycolysis signature (Glyc.Sig), might help in predicting the ICI response. Therefore, we utilized 41 pan-cancer scRNA-seq datasets to develop the Glyc.Sig. Spearman correlation was initially adopted between enrichment scores of tumor glycolysis and gene expression levels. The Gx gene set contained genes with a positive correlation with scores of glycolysis intensity (R > 0 and FDR < 1e−05). The Gy gene set comprised significantly overexpressed genes in malignant cells. We intersected Gx and Gy to formulate the Gn gene set within each dataset (n = 1, 2,…, 41) (18), which contained upregulated genes of tumor specificity with a positive correlation with glycolysis intensity. For each gene in the G1–G41 gene sets, we calculated the geometric mean of the Spearman R value, and genes with a final R mean greater than 0.3 (88) were selected to form Glyc.Sig. The detailed gene list is demonstrated in Supplementary Table S6 . The simplified process of the generation of Glyc.Sig could be intuitively observed in Figure 4A . Additionally, pathway analysis based on Reactome ( Figure 4B ) and KEGG ( Figure 4C ) was adopted to explore the biological functions of Glyc.Sig. The overrepresented pathways were primarily associated with glycolysis, hypoxia, the KEAP1-NFE2L2 pathway, G1/S DNA damage checkpoints, glutathione metabolism, and carbohydrate metabolism processes.

We first investigated the relations between Glyc.Sig and 75 immune-related genes (58), consisting of the HLA set, stimulatory set, and inhibitory set. A general negative correlation was found across almost all 75 genes in 30 different cancer types ( Figure 4D ). Figure 4E demonstrates the situation of immune cell infiltration. A higher level of glycolysis intensity was negatively correlated with various categories of immune cell infiltration, including cytotoxic cells (CD8 T cell, cytotoxic lymphocyte, NK cells), B lineage cells, and myeloid dendritic cells. The analysis above suggested a decreased antitumor immunity in tumors with a high intensity of glycolysis.

Next, we explored the links between enrichment of HALLMARK pathways and expression intensity of Glyc.Sig. Both of them were calculated through GSVA. All of the top 10 HALLMARK pathways, including DNA repair, MYC signaling, and glycolysis, were positively correlated with the expression intensity of Glyc.Sig ( Figure 5A ). All of these pathways might be associated with a weak immune response based on previous studies (89–91). We also explored the relations between enrichment of tertiary lymphoid structure (TLS)-related genes and the expression intensity of Glyc.Sig. As shown in Figure 5B , TLS scores as well as most of the TLS-related genes were negatively related to the expression intensity of Glyc.Sig. Notably, high TLS scores and related genes usually indicated an abundant immune cell infiltration. Additionally, we investigated the relations between the median score of Glyc.Sig and the median ITH as well as the median TMB. ITH was a feature that facilitated immunosuppression (92). As anticipated, the expression intensity of Glyc.Sig was positively correlated with ITH (R = 0.64, p = 0.00014, Figure 5C ). However, a positive link between TMB and Glyc.Sig was found (R = 0.8, p = 9.1e−8, Figure 5D ), which seemed to go against our current understanding of TMB. Generally, higher TMB means better immune response because of abundant antigenicity, while a high level of Glyc.Sig does the opposite. Based on the median GSVA score of Glyc.Sig and the median TMB as grouping criteria, patients were divided into four distinct subgroups: high Glyc.Sig/high TMB (HGHT), high Glyc.Sig/low TMB (HGLT), low Glyc.Sig/high TMB (LGHT), and low Glyc.Sig/low TMB (LGLT). First, we clarified that both high Glyc.Sig and low TMB had decreased cytotoxic lymphocyte infiltration ( Figure 5E , Supplementary Figure S1A ), which was consistent with our current understanding. We further compared immune cell infiltration among the LGLT, LGHT, HGLT, and HGHT subgroups. LGHT was identified as the group with the highest infiltration of cytotoxic lymphocytes and NK cells, while HGLT was found to have the lowest infiltration of most immune cells, including cytotoxic lymphocytes, monocytic lineage, CD8 T cells, NK cells, T cells, and B lineage ( Figure 5F , Supplementary Figure S1B ). As for the LGLT and HGHT subgroups, the situation seemed to be unclear, possibly for the reason that these two subgroups had both immune-stimulating factors and inhibiting factors.

To explore the prognostic value of Glyc.Sig, bulk RNA-seq data, and the corresponding clinical data from 10 ICI cohorts were curated and analyzed to construct the models for ICI response prediction based on Glyc.Sig. To better utilize these data, we divided them into three subsets: the training set for model training and parameter tuning based on five-fold cross-validation repeated 10 times (except cancerclass); a validation set for model comparison and selection according to AUC; and the independent testing set for final model diagnosis. The flowchart reflected the whole process ( Figure 6A ). Seven different models were incorporated and trained at the beginning; naive Bayes was finally chosen as the Glyc.Sig model with the highest AUC of 0.69 ( Figure 6B ). Next, we utilized the independent testing set to make further evaluation of the Glyc.Sig model by predicting the ICI response. An AUC of 0.66 was finally observed ( Figure 6C ).

To evaluate the predictive value of Glyc.Sig on overall survival (OS), patients were divided into a high-risk group (predicted non-responders, NR) and a low-risk group (predicted responders, R) according to the prediction of ICI response given by the Glyc.Sig model. Kaplan–Meier survival analysis demonstrated that in both validation and testing cohorts, a high-risk group had a significantly shorter OS than the low-risk group (p < 0.01, Figure 6D ). In the validation cohorts, the low-risk group exhibited a median OS of 29.9 months (HR = 1.83, 95% CI 1.17–2.86), nearly doubling the OS (14.42 months) observed in high-risk counterparts. In the testing cohort, low-risk patients had a 31.7-month median OS, while the OS in high-risk subjects was 11.23 months (HR = 2.28, 95% CI: 1.27–4.09). Subsequently, we examined the performance of the Glyc.Sig predictive model in each of the 10 ICI RNA-seq cohorts. The AUC across these 10 cohorts ranged from 0.53 to 0.91 ( Supplementary Figure S2 ), with the AUC of the SKCM cohort in 2017 by Riaz reaching the highest (AUC = 0.91, 95% CI: 0.82–0.99), followed by a UC cohort in 2018 by Mariathasan (AUC = 0.90, 95% CI: 0.87–0.94). The GBM cohort in 2019 by Zhao exhibited the lowest predictive accuracy (AUC = 0.53, 95% CI: 0.24–0.82), which could be attributed to its restricted sample size. As for survival analysis, the GC cohort in 2018 by Kim was excluded due to the lack of OS data. For the remaining nine cohorts, high-risk patients predicted by the Glyc.Sig model showed significant survival benefits. Notable observations were made in the SKCM cohort from 2015 (Van Allen), the UC cohort from 2017 (Synder), the UC cohort from 2018 (Mariathasan), the SKCM cohort from 2019 (Liu), and the SKCM cohort from 2019 (Gide) ( Supplementary Figure S3 ).

Comparison between Glyc.Sig and other well-constructed predictive gene signatures was made by calculating the AUC for each gene signature. In comparison with other gene signatures (81–85), Glyc.Sig performed best in the validation cohort ( Figure 6E ) with an AUC of 0.69. Additionally, Glyc.Sig manifested the best performance in most of the 10 individual ICI cohorts, while most of the other signatures performed ideally in only one or two cohorts ( Figure 6F , Supplementary Table S7 ). We also compared the average AUC of these gene signatures over the 10 individual ICI cohorts. Notably, Glyc.Sig dominated the list with an average AUC of 0.76, surpassing the second signature (INFG.sig) with an AUC of 0.62 ( Figure 6G ).

Immune response data corresponding to gene knockouts were collected from seven CRISPR cohorts. These cohorts were further categorized into 17 datasets based on the model cells and treatments adopted. A total of 22,505 genes were documented across these datasets. By ranking the genes according to their average Z-scores, we identified the top-ranked genes as immune-resistant, suggesting that their knockout could enhance antitumor immunity. Conversely, the bottom-ranked genes were classified as immune-sensitive, indicating that their knockout might suppress antitumor immunity. The ranking process is depicted in Figure 7A . Out of the 22,505 genes, the numbers of genes in the top 5%, 10%, and 20% rankings were 1,125, 2,250, and 4,501, respectively. To further investigate the immunotherapeutic value of Glyc.Sig, we subsequently calculated and compared the proportion of the top 5%, 10%, and 20% ranked genes in Glyc.sig and other previously reported immune-resistant gene signatures, including TcellExc.Sig, ImmuneCells.Sig, IMS.Sig, LRRC15.CAF.Sig, and CRMA.Sig (18, 83, 93–95). The IPRES signature was excluded due to its unique composition of pathways rather than genes. As anticipated, Glyc.Sig had the highest proportion of overlapped top-ranked genes compared to other gene signatures ( Figure 7B ). The immune-resistant genes (those in the top 20%) were significantly enriched in Glyc.Sig (p = 0.02, Fisher’s exact test). Twenty-seven genes in Glyc.Sig, including PSMF1, GCLC, CES1, MAF1, LDHA, KLHL24, FCGR2B, POLR2K, MCCC1, GPI, COMMD8, C19orf24, DNAJC12, YDJC, HSPB7, DERL1, NDUFA4L2, AGTRAP, RNF128, MRPL10, EIF5A2, B3GNT3, SSR1, DNAJC2, SDCCAG8, TPI1, and CKAP4, ranked in the top 20%. The relations between tumor immune regulation and these glycolysis-related genes were verified using a series of independent CRISPR datasets ( Figure 7C , Supplementary Figure S4 ), indicating their potential as therapeutic targets for combination with immune checkpoint blockade (ICB).

Given the upregulation of glycolysis and the predictive value of Glyc.Sig in FH-deficient RCC and pan-cancer data, we hypothesized that the hub gene of CRISPR cohorts, LDHA, might play a key role in this process. LDHA encodes the final step of glycolysis, where pyruvate is converted into lactate, generating ATP. Analysis of the FH-deficient RCC scRNA-seq dataset revealed that LDHA expression was notably higher in FH-deficient RCC epithelial cells ( Figure 7D ) compared to other cell types. Additionally, when comparing LDHA expression across metastasis, primary tumors, and adjacent normal tissues using the GSE157256 dataset (15), we found significantly lower LDHA levels in adjacent normal tissues ( Figure 7E ). To further confirm these results, immunoblotting analysis and IHC staining from patients at Renji Hospital showed significantly higher LDHA expression in tumor tissues compared to normal adjacent tissues ( Figures 7F, G ).