The glioma and other cancer cohorts used in this study were sourced from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/), Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn/), Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), and Kaplan-Meier Plotter (https://kmplot.com/) databases, with detailed information provided in Table 1 . According to previously published study (13, 14), a total of 354 lactylation-related genes were included and were presented in Supplementary Table S1 .

To construct a lactylation-related gene signature, we first utilized the TCGA glioma cohort to screen for genes associated with patient prognosis, conducting batch survival analysis using the R package “survival”. The resulting prognostic genes were intersected with 354 lactylation-related genes, followed by further refinement of the overlapping genes through LASSO regression analysis using the R package “glmnet”. Univariate and multivariate Cox proportional hazards regression analyses were employed to develop the LRGs model, with Cox regression performed using the R package “survival” and forest plot visualization achieved through the R package “ggplot2”. The chromosomal locations of the four LRGs core genes (KIF2C, CALD1, HSPE1, and IFI16) were visualized using the R package “circlize”, and Spearman’s correlation analysis was employed to examine the relationships among these genes. Glioma patients were stratified into LRGs-high and -low groups using the median LRGs score as a cutoff value. Kaplan-Meier survival analysis and visualization were performed using the R packages “survival” and “survminer” respectively. Scatter plot was generated with the R package “ggplot2”, while time-dependent ROC was produced and visualized using the R packages “timeROC” and “ggplot2”.

To validate the expression patterns of the LRGs, we utilized the TCGA and GSE16011 cohorts to compare the differences of the LRGs score and its core genes between normal brain and glioma tissues. Subsequently, multiple independent cohorts, such as CGGA301, CGGA325, CGGA693, GSE4412, and GSE43378, were employed to evaluate the prognostic value of the LRGs. The methodologies for generating Kaplan-Meier survival curves, scatter plots, and time-dependent ROC curves were consistent with those described previously.

The TCGA and CGGA693 glioma cohorts were employed to compare differences in gender, age, WHO grade, IDH status, and 1p/19q co-deletion between LRGs-high and -low groups of glioma patients. Based on the TCGA glioma cohort, we incorporated five variables, including age, WHO grade, IDH status, 1p/19q co-deletion, and LRGs score, to construct a predictive nomogram model. Cox regression analysis was performed using the R package “survival”, while the nomogram was developed and calibration curves were plotted using the R package “rms”. Time-dependent AUC and Decision curve analysis (DCA) were analyzed with the R packages “timeROC” and “stdca.R”, respectively, with visualization implemented using the R package “ggplot2”.

Glioma patients in the TCGA cohort were divided into LRGs-high and -low groups according to the median LRGs score as the cutoff value, and differential expression analysis of the original Counts matrix was performed using the R package “DESeq2”. Differentially expressed genes (DEGs) were subjected to enrichment analysis through the Metascape platform (https://metascape.org/). Concurrently, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) were conducted on the DEGs using the R package “clusterProfiler” for analysis and R package “ggplot2” for visualization.

The immune cell infiltration scores for glioma patients in the TCGA cohort were obtained from the CIBERSORTx online platform (https://cibersortx.stanford.edu/) based on the Cibersort core algorithm. The infiltration of myeloid-derived suppressor cells (MDSCs) and cancer-associated fibroblasts (CAFs) were derived from the Tumor Immune Dysfunction Exclusion (TIDE) online platform (http://tide.dfci.harvard.edu/). The Stromal, Immune, and ESTIMATE Scores were sourced from the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) database (https://bioinformatics.mdanderson.org/estimate/). We compared the differences in the aforementioned parameters between the LRGs-high and -low groups, and evaluated the correlations between the LRGs core genes and these parameters. The immune subtypes in glioma were sourced from a previously published study (15). Kaplan-Meier survival analysis was performed to assess the significance of immune subtypes in predicting patient survival. The single-cell resolution data from the GSE131928 cohort were sourced from the TISCH2 platform (http://tisch.comp-genomics.org/). We analyzed the differentially expressed genes across various malignant cell types, as well as the expression abundance of LRGs core genes in different cells. Additionally, the biological functions of highly expressed genes in malignant cells were investigated through online analysis using the Metascape platform.

The Cancer-immunity cycle was obtained from the Tracking Tumor Immunophenotype (TIP) database (http://biocc.hrbmu.edu.cn/TIP/), which consists of seven steps: (1) release of cancer cell antigens, (2) cancer antigen presentation, (3) priming and activation, (4) trafficking of immune cells to tumors, (5) infiltration of immune cells into tumors, (6) recognition of cancer cells by T cells, and (7) killing of cancer cells. Next, we examined the correlation between the LRGs core genes and the expression of immune checkpoints, such as CTLA4 and PDCD1 etc., using sequencing data from glioma patients in the TCGA cohort. Subsequently, glioma patients were stratified into four groups based on immune checkpoint expression levels combined with LRGs score, and the prognostic significance of these combinations was assessed through Kaplan-Meier survival analysis. Furthermore, we evaluated the predictive capability of the LRGs for immunotherapy efficacy using the TIDE algorithm and multiple cohorts, including glioma (PRJNA482620), melanoma (GSE91061), and atezolizumab pan-cancer (KM Plotter) datasets.

Integrated mutation data for glioma samples were obtained from the GDC portal (https://portal.gdc.cancer.gov/) and analyzed using the R package “maftools” to identify the top 15 most frequently mutated genes between LRGs-high and -low groups. Subsequently, we compared overall survival differences between patients with mutant- and wild-type variants of genes, such as TTN and EGFR etc., and analyzed expression differences of the LRGs core genes between mutant- and wild-type groups.

Tumor mutation burden (TMB) data were calculated for TCGA glioma samples using the R package “maftools”. Microsatellite instability (MSI) data were referenced from the study by Russell Bonneville et al. (16). Data on tumor purity, tumor ploidy, homologous recombination deficiency (HRD), and neoantigens were obtained from the research by Vesteinn Thorsson et al. (15). Tumor stemness data were derived from previous studies (17), we acquired six stemness indicators based on mRNA expression and methylation signatures, including RNAss, EREG-METHss, DMPss, ENHss, and EREG-EXPss.

For chemotherapy sensitivity analysis, we extracted expression matrix data of glioma patients from the TCGA cohort and utilized the core algorithms of the R packages “oncoPredict” and “pRRophetic” to calculate drug sensitivity scores by integrating drug and cell line expression profiles provided by these tools. For radiotherapy sensitivity analysis, we extracted data from the TCGA glioma patients who received radiotherapy, including those with progressive disease/stable disease (PD/SD) and partial response/complete response (PR/CR).

siRNAs were purchased from Sangon Biotech. The sequences of the siRNAs were as follows: siCtrl: 5’-UUCUCCGAACGUGUCACGUTT-3’; siHSPE1-1#:5’-GCAGGACAAGCGUUUAGAATT-3’; siHSPE1-2#:5’-GAGUGCUGCUGAAACUGUATT-3’. Glioma cells were seeded in 12-well plates and transfection was initiated when the cell density reached approximately 70%. Transfection was performed using Lipofectamine 3000 according to the manufacturer’s instructions. Transfection efficiency was evaluated by qRT-PCR. Total RNA was extracted using Trizol reagent, and reverse transcription was conducted using M-MLV Reverse Transcriptase. Quantitative real-time PCR was carried out in triplicate using the SYBR Green Master Mix. The following primers were used for qPCR: GAPDH-F: 5’-GTCTCCTCTGACTTCAACAGCG-3’, GAPDH-R: 5’-ACCACCCTGTTGCTGTAGCCAA-3’; HSPE1-F: 5’-GCTGAAACTGTAACCAAAGGAGG-3’, HSPE1-R: 5’- TCTCCAACTTTCACGCTAACTGG-3’.

Western blot analysis was performed as previously described (18). Briefly, glioma cells were lysed using RIPA lysis buffer supplemented with a 1× protease inhibitor cocktail. Equal amounts of protein were separated by 12% SDS-PAGE and transferred to a 0.22 μm PVDF membrane. The membranes were blocked with 5% skimmed milk and then incubated overnight at 4°C with primary antibodies: HSPE1 (Affinity, AF0183) and β-actin (Abmart, P30002). Following primary antibody incubation, the membranes were incubated with the secondary antibody at room temperature for 1 hour. Protein bands were visualized using the Super ECL detection reagent.

Cell proliferation was assessed using the Cell Counting Kit-8 (CCK-8), following the manufacturer’s protocol. Absorbance at 450 nm was measured at the indicated time points to determine cell viability. For colony formation assays, 1×10³ cells were seeded in 6-well plates and incubated for 14 days. After incubation, cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet to visualize the colonies. For the invasion assay, 5×10⁴ cells were seeded into Transwell inserts with an 8 μm pore size, precoated with Matrigel. After 24 hours of incubation, the invasive cells were stained with 0.1% crystal violet, observed under a microscope, and the number of invasive cells was counted using ImageJ software.

Apoptosis was analyzed using Annexin V-FITC/propidium iodide (PI) double staining followed by flow cytometry. Briefly, harvested cells were resuspended in binding buffer and incubated with Annexin V-FITC and PI for 10 min at room temperature in the dark. The apoptotic cell population was quantified by flow cytometer, and data were analyzed using FlowJo software.

All statistical analyses were processed on R Studio (V4.2.1) or GraphPad Prism 8 software, and P value < 0.05 indicated statistically significant differences. The quantitative results are presented as the mean ± standard deviation (SD). Wilcoxon rank sum test was used for unpaired samples, t-test was used for paired samples, and ANOVA was used for comparisons between multiple groups. Log Rank P test was used for Kaplan-Meier survival analysis. Spearman test was used for Correlation analysis.

The workflow of this study was depicted in Figure 1 . To establish a LRGs signature predictive of glioma prognosis, we employed the TCGA glioma cohort as the training set. Initial analysis identified 7,622 risk-associated genes and 5,264 protective genes ( Figure 2A ). Intersection of the 7,622 risk-associated genes with 354 known lactylation-related genes yielded 147 overlapping candidates ( Figure 2B ). Subsequent LASSO regression refined this set to 28 potential prognostic genes ( Figures 2C, D ). Differential expression analysis revealed significant upregulation of nearly half these genes in glioma versus normal brain tissues ( Figure 2E ). Through integrated univariate and multivariate Cox proportional hazards regression analyses, four lactylation-related genes—KIF2C, CALD1, HSPE1, and IFI16—emerged as core biomarkers, forming the basis of the LRGs signature ( Figures 2F, G ). The specific Cox regression coefficients were applied in the formula: LRGs score = (0.20 × KIF2C expression) + (0.48 × CALD1 expression) + (0.75 × HSPE1 expression) + (0.17 × IFI16 expression) – 8 (constant) ( Figures 2H ). These four genes localized to distinct chromosomal loci: KIF2C (1p34), CALD1 (7q33), HSPE1 (2q33), and IFI16 (1q23) ( Figure 2I ). Heatmap analysis demonstrated significant positive co-expression correlations among all four genes within glioma tissues ( Figure 2J ). Scatterplot illustrated that high-LRGs score patients experienced both elevated mortality rates and reduced survival durations ( Figure 2K ). Consistent with this, Kaplan-Meier survival curve revealed markedly inferior overall survival for patients stratified into the LRGs-high group compared to the -low group ( Figure 2L ). The LRGs signature exhibited robust predictive capacity, with time-dependent ROC analysis yielding area under the curve (AUC) values of 0.848 (95%CI: 0.809 - 0.887), 0.894 (95%CI: 0.860 - 0.928), and 0.827 (95%CI: 0.775 - 0.880) for 1-, 3-, and 5-year survival, respectively ( Figure 2M ).

To validate the prognostic robustness of the LRGs signature, we first quantitatively evaluated expression patterns in TCGA and GSE16011 cohorts. This analysis confirmed significantly upregulation of LRGs score and its core genes (KIF2C, CALD1, HSPE1, and IFI16) in glioma versus normal tissues ( Figures 3A, B ). We subsequently validated the LRGs prognostic value in five independent glioma cohorts, including CGGA301, CGGA325, CGGA693, GSE4412, and GSE43378. Patients stratified into LRGs-high group exhibited consistently inferior overall survival and elevated mortality incidence compared to -low group. Time-dependent ROC analyses demonstrated sustained predictive accuracy, with AUC values exceeding 0.700 for 3- and 5-year survival probabilities in all validation cohorts ( Figures 3C–G ). Further quantification via Harrell’s concordance index yielded C-index values of 0.809 (95%CI: 0.797 - 0.822; TCGA), 0.747 (95%CI: 0.731 - 0.762; CGGA325), 0.704 (95%CI: 0.686 - 0.722; CGGA301), 0.659 (95%CI: 0.644 - 0.673; CGGA693), 0.658 (95%CI: 0.620 - 0.695; GSE43378), and 0.618 (95%CI: 0.581 - 0.655; GSE4412), affirming the model’s discriminative capacity ( Figure 3H ). Comparative analysis against established prognostic signatures revealed that our LRGs model achieved comparable prognostic efficacy with greater parsimony, outperforming multi-gene constructs including those by Zhang Q et al. (16 genes) (19), Zhang N et al. (5 gene pairs) (20), and Zhang M et al. (10 genes) (21) in feature economy ( Figures 3I–K ).

To optimize the precision of our LRGs signature, we systematically evaluated clinicopathological covariates influencing glioma outcomes. Comparative analysis of TCGA and CGGA693 cohorts revealed significant differences in multiple clinical characteristics between LRGs-high and -low glioma patient groups—including age, WHO grade, IDH status, and 1p/19q co-deletion—but showed no association with gender distribution ( Figures 4A, B ). Next, univariate and multivariate Cox regression analyses confirmed LRGs score, age, WHO grade, IDH status, and 1p/19q codeletion as independent risk factors ( Table 2 ). These variables were subsequently integrated into a comprehensive nomogram, with calibration curves demonstrating excellent agreement between predicted and observed survival probabilities ( Figures 4C, D ). The composite model achieved a Harrell’s concordance index (C-index) of 0.860 (95% CI: 0.850 - 0.871). Time-dependent ROC analysis further validated its predictive efficacy, with AUC values exceeding 0.800 for 1- to 5-year survival predictions ( Figure 4E ). DCA demonstrated that the nomogram provided superior clinical net benefit compared to using the LRGs signature alone across most threshold probabilities, particularly at 3- and 5-year time points, outperforming both all-treat and all-none reference strategies ( Figures 4F–H ). Collectively, the integrated nomogram demonstrated enhanced prognostic capability and clinical applicability for glioma risk stratification.

Given the significant prognostic association of the LRGs signature, which prompting further investigation into its functional underpinnings in glioma biology. Differential expression analysis between LRGs-high and -low groups identified 1,169 significantly upregulated and 803 downregulated genes ( Figure 5A ). Notably, these upregulated genes demonstrated significant enrichment in immune signaling pathways, including inflammatory response, cytokine signaling, and immune response functions, while downregulated genes were primarily implicated in synaptic function and neurotransmitter transmission ( Figures 5B, C ). GO analysis revealed distinct compartmentalization characteristics, with cellular components primarily localized to synaptic membranes, collagen-containing extracellular matrices, and transmembrane transporter complexes ( Figure 5D ). Molecular functions predominantly involved passive transmembrane transporter activity, gated channel function, and cytokine binding ( Figure 5D ). Biological process annotation demonstrated positive regulation of cytokine-mediated signaling, T cell activation, and interferon-γ response pathways, while negatively regulating trans-synaptic signaling, membrane potential modulation, and neurotransmitter transport ( Figure 5E ). KEGG pathway analysis further confirmed involvement in oncogenic and immune signaling cascades, including PI3K-AKT, IL-17, and JAK-STAT pathways ( Figure 5F ). GSEA of hallmark gene sets verified significant associations with cytokine signaling, cell cycle checkpoint regulation, and interferon response pathways ( Figures 5G–I ).

Given the established connection between the LRGs and immune pathways, we systematically profiled tumor-infiltrating immune cells using the Cibersort algorithm. Patients in the LRGs-high group exhibited significant enrichment of immunosuppressive populations, including M0, M1, and M2 macrophages, neutrophils, regulatory T cells (Tregs), γδ T cells, and resting memory CD4⁺ T cells. In contrast, other patients showed preferential infiltration of monocytes, activated NK cells, and activated mast cells ( Figure 6A ). The expression of LRGs core genes demonstrated significant positive correlations with pro-tumoral macrophages (particularly M2-polarized subtypes), neutrophils, Tregs, and γδ T cells, while being inversely correlated with monocytes, memory B cells, and naïve CD4⁺ T cells ( Figure 6B ). Complementary analyses using the TIDE and ESTIMATE algorithms further revealed increased infiltration of MDSCs and CAFs, along with elevated stromal, immune, and ESTIMATE scores in glioma with high LRGs score ( Figures 6C, D ). Each of these ESTIMATE metrics showed positive correlations with the expression of all four LRGs core genes ( Figures 6E ). Immunophenotypic stratification analysis revealed that patients in the LRGs-low group predominantly exhibited the immunologically favorable C5 subtype, while those in the -high group were significantly enriched for the poor-prognosis C4 subtype ( Figures 6F, G ).

To elucidate the cellular expression patterns of LRGs core genes, we utilized the TISCH2 platform to analyze the GSE131928 cohort, with a focus on genes specifically expressed in malignant cell populations ( Figure 6H ). The results revealed 27 significantly upregulated and 55 downregulated genes in AC-like malignant cells; 246 upregulated and 147 downregulated genes in MES-like malignant cells; 176 upregulated and 269 downregulated genes in NPC-like malignant cells, and 61 upregulated and 172 downregulated genes in OPC-like malignant cells ( Figure 6I ). We further extracted these upregulated genes for enrichment analysis and found that they were primarily involved in regulating cell cycle and oncogenic pathways, such as the FOXM1 and VEGFA pathways ( Figure 6J ). The LRGs core genes were specifically expressed in various malignant cell populations, with HSPE1 showing particularly prominent expression. Additionally, IFI16 was also expressed in monocyte/macrophage lineages and exhausted T cell compartments ( Figure 6K ). Based on these findings, we extracted the top 10 upregulated genes expressed in each malignant cell subtype and observed significant positive correlations with the four LRGs core genes, except for B4GALNT1, KIF5A, and TSFM in NPC-like malignant cells and PIP4K2A in OPC-like malignant cells ( Figure 6L ).

To evaluate the role of the LRGs in tumor immunity, we characterized cancer-immunity cycle dysregulation between LRGs-high and -low groups. Despite elevated tumor antigen release, patients in the LRGs-high group exhibited significantly impaired antigen presentation and processing capacity ( Figures 7A–C ). Concurrently, these patients demonstrated heightened immune cell recruitment to tumor peripheries but diminished intratumoral infiltration ( Figures 7D, E ). Although T-cell recognition of tumor cells remained comparable between the two groups, patients in the LRGs-high group exhibited significantly reduced cytotoxic efficacy against cancer cells ( Figures 7F, G ). Subsequent analysis revealed significant positive correlations between the LRGs core genes and multiple immune checkpoints, with these checkpoints demonstrating marked upregulation in the LRGs-high group ( Figures 7H, I ). Kaplan-Meier analysis demonstrated poorer survival outcomes in patients with concurrently elevated expression of both the LRGs and immune checkpoints, including CTLA-4, PDCD1, CD274, HAVCR2, PDCD1LG2, TNFRSF4, or TNFRSF18 ( Figures 7J–P ). Multi-platform validation incorporating TIDE algorithmic assessment, glioma, melanoma, and atezolizumab pan-cancer immunotherapy cohorts consistently showed that patients in the LRGs-high group exhibited elevated TIDE score, shortened survival duration, and attenuated therapeutic responsiveness to immune checkpoint blockade ( Figures 7Q-T ). These findings collectively established the LRGs signature as a biomarker for response to immunotherapy in glioma.

Given the established link between genetic alterations and oncogenesis, we examined LRGs-associated genomic variations in glioma progression using comprehensive mutational profiling ( Figure 8A ). Mutations in critical driver genes, such as TTN, EGFR, PTEN, NF1, and FLG, demonstrated significant association with adverse overall survival outcomes ( Figure 8B ). The expression levels of LRGs core genes were significantly elevated in most oncogenic mutation groups, with the exception of HSPE1 suppression in FLG-mut group ( Figure 8C ). Conversely, mutations in IDH1, CIC, FUBP1, or NOTCH1 conferred favorable prognoses, correlating with downregulation of these LRGs-associated markers, though IFI16 expression remained unaltered in NOTCH1-mut group ( Figures 8D, E ). These findings position LRGs as molecular integrators of genetic variation influencing glioma pathogenesis.

Patients with elevated tumor stemness scores tend to have poorer prognosis, and these scores can predict tumor metastasis and recurrence. Comprehensive analysis revealed that the LRGs and its core genes showed positive correlations with DNAss, EREG-METHss, DMPss, ENHss, and EREG-EXPss, but negative correlation with RNAss, suggesting that the LRGs may promote tumor cell dedifferentiation and enhance stemness ( Figures 9A–C ). Tumor patients exhibiting high heterogeneity levels face increased risks of recurrence and mortality. Results of heterogeneity analysis revealed that the LRGs and their core genes were positively correlated with TMB, HRD, and LOH, while showing negative correlations with MATH and MSI, suggesting that tumors with high LRGs score exhibited greater genomic instability, higher mutational burden, and potential homologous recombination deficiency ( Figures 9D–F ). Chemosensitivity analysis indicated that patients in the LRGs-high group exhibited enhanced sensitivity to temozolomide, cisplatin, oxaliplatin, irinotecan, teniposide, 5-fluorouracil, and gemcitabine, but developed resistance to olaparib, afatinib, paclitaxel, gefitinib, and erlotinib ( Figure 9G ). Notably, the LRGs score was significantly elevated in the radiotherapy non-responder (PD/SD) group compared to the responder (PR/CR) group, further establishing LRGs as biomarkers of therapeutic resistance Figure 9H .

Given that HSPE1, a core LRGs gene, demonstrated the highest hazard ratio in multivariate Cox analysis, exhibited high expression abundance in malignant cells, and has an undefined functional role in gliomagenesis, we selected HSPE1 for further investigation into its ability to affect glioma cells. In U87 and U251 glioma cell models, we successfully knocked down HSPE1 expression, as verified by qPCR and Western blot analysis ( Figures 10A, B ). Functional characterization demonstrated that HSPE1 ablation significantly attenuated malignant phenotypes, suppressing cellular proliferation in CCK-8 assay ( Figure 10C ), reducing clonogenic capacity in colony formation assay ( Figures 10D, E ), and markedly attenuating invasive potential in transwell assay ( Figures 10F, G ). However, HSPE1 knockdown showed no significant effect on the apoptosis of glioma cells ( Supplementary Figure S1 ). Complementing these deficits observed in in vitro functional assays, clinical validation utilizing the Human Protein Atlas (HPA) database revealed significant overexpression of HSPE1 in glioma specimens compared to normal brain tissue, which was further confirmed by immunohistochemical (IHC) analysis of tissue microarrays ( Figures 10H, I ). In summary, these integrated mechanistic and clinical findings establish HSPE1 as a critical mediator of gliomagenesis within the lactylation-related oncogenic network.