The TCGA-CESC dataset (accessed December 20, 2024) contains 306 CC tumor tissues and survival data from 293 samples. (https://portal.gdc.cancer.gov/) (18). Samples with incomplete or missing clinicopathological information were excluded. Differential gene expression analysis was performed using the R package ‘limma’ (version 3.64.3) (19). All datasets were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), including three microarray datasets [GSE122697 (20), GSE146114 (21) and GSE63514 (22)] for WGCNA analysis, the scRNA-seq dataset GSE44001 (23) and the validation set GSE30760 (24). The GSE7803 (25) dataset was used for independent external validation of the prognostic model. At present, no other publicly available scRNA-sequencing datasets with sufficient sample size, annotation quality and relevance to cervical cancer are available; therefore, GSE44001 was selected as the only suitable dataset for single-cell analysis in the present study. Detailed information about the datasets is provided in Table SI. Batch effects were mitigated using the ‘sva’ package (version 3.56.0) (26). Furthermore, 325 LAGs were collected from previously published studies, as detailed in Table SII (27–29).

To further assess the association between the prediction model and the immune response of the patients, the present study integrated data from the IMvigor210 cohort using the ‘IMvigor210CoreBiologies’ R package (version 1.15.0.) (30). The present study also obtained cervical cancer genome-wide association studies (GWAS) data and LAG expression quantitative trait loci (eQTL) GWAS data from the FinnGen database (https://risteys.finngen.fi/) for MR analysis (31). The cervical cancer GWAS dataset, ‘FinnGen R12 CD2 INSITU CERVIX UTERI EXALLC’ contains data on 20,104,836 single nucleotide polymorphisms (SNPs).

scRNA-seq data from both normal and cervical cancer tissues were obtained from the GSE168652 dataset. Data preprocessing and quality control were performed using the ‘Seurat’ R package (version 5.3.1; http://satijalab.org/seurat/) (32). Cells with low gene counts or high mitochondrial gene percentages were filtered out to ensure data integrity. Batch effects across samples were corrected using the ‘Harmony’ package (version 1.2.4; http://github.com/immunogenomics/harmony) (33). Following normalization and scaling, unsupervised cell clustering was performed at a resolution of 0.8 and dimensionality reduction was visualized using t-distributed stochastic neighbor embedding. To assess the activity of lactylation-related gene sets in individual cells, the ‘AddModuleScore’ function in Seurat was applied. In addition, intercellular communication networks were analyzed using the ‘CellChat’ R package to explore signaling interactions among different cell populations (version 2.2.0; http://github.com/jinworks/CellChat) (34).

Differentially expressed genes (DEGs) were identified using three GEO datasets: GSE122697, GSE146114 and GSE63514. These datasets were chosen for their high-quality, well-annotated expression profiles and adequate sample sizes, ensuring data reliability and reproducibility. The grouping details for each dataset are provided in Table SI. To minimize batch effects, the ComBat function from the ‘sva’ R package was applied (version 3.22) (26). Subsequently, a disease-specific co-expression network was constructed using the ‘WGCNA’ R package (version 1.73) (35). The optimal soft-thresholding power (β) was selected to ensure scale-free topology and modules most correlated with cervical cancer were identified for downstream analyses.

The highly disease-associated genes identified by WGCNA were cross-analyzed with differentially expressed genes and lactation-related genes to identify potential LAG-associated prognostic candidate genes. Subsequently, these genes were subjected to feature selection and prognostic modeling using multiple machine learning algorithms, including least absolute shrinkage and selection operator (LASSO) (36), Stepwise Cox (Backward), CoxBoost (37), random survival forest (RSF) (38), generalized boosted regression model and survival support vector machine (18–21). Each algorithm was independently evaluated for the prognostic performance of the candidate genes in TCGA training cohort and GSE30760 validation cohort. The concordance index (C-index) was used to quantify model performance and the models were ranked according to their average C-index values. The StepCox (Backward) model achieved the highest predictive performance [area under the curve (AUC)=0.82; C-index=0.79], outperforming CoxBoost and survival SVM while maintaining simplicity and interpretability. To minimize overfitting and ensure robustness, the StepCox (Backward) model was lastly selected for downstream analysis, which is consistent with previous studies demonstrating the stability and strong generalization ability of this model in cancer prognosis modeling (39).

After obtaining the LAGs risk score of each patient, patients were divided into high-risk and low-risk groups based on the median of the LAGs risk score distribution within the cohort. Survival outcomes, including progression-free survival (PFS), disease-specific survival (DSS) and disease-free interval (DFI), were compared between the two groups using Kaplan-Meier analysis implemented in the ‘survminer’ R package (version 0.5.1; http://rpkgs.datanovia.com/survminer/index.html). The predictive accuracy of the model was further assessed by calculating the AUC for 1-, 3- and 5-year overall survival (OS) using the ‘timeROC’ R package (version 0.4) (40). Differences in risk scores across clinical subgroups, such as tumor (T), lymph node (N) and metastasis (M) stage and histological grade (41), were examined using TCGA-CESC clinical data. To determine whether the LAGs-based risk model serves as an independent prognostic factor, univariate and multivariate Cox regression analyses were performed. Lastly, a nomogram integrating the LAGs risk score and key clinical variables was constructed to predict individual patient survival probabilities and enhance the model's clinical applicability. Simultaneously, the GEPIA2.0 database (http://gepia2.cancer-pku.cn/) (42) was used to validate the survival analysis.

To explore the biological pathways underlying differences between the high- and low-risk groups, Gene Set Enrichment Analysis (GSEA) was performed (43). Additionally, enrichment scores were quantified at the individual sample level using single-sample gene set enrichment analysis (ssGSEA) implemented in the R package GSVA (version 2.2.1) on the TCGA-CESC (Cervical squamous cell carcinoma and endocervical adenocarcinoma) dataset downloaded from The Cancer Genome Atlas (TCGA) via the GDC Data Portal, to assess pathway activity (44). Functional annotations were further investigated using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to identify biological processes and signaling pathways associated with the LAGs risk model. Pathways with P<0.05 were considered statistically significant.

To examine immune infiltration patterns across LAGs risk subgroups, the ESTIMATE analysis was performed to determine the malignant tumor stroma and immune cell (ESTIMATE) scores of patients (version 4.7.0) (45). The CIBERSORT algorithm (https://cibersortx.stanford.edu/) was then used to estimate the relative abundance of immune cell populations, while ssGSEA was applied to assess immune cell signature enrichment and immune-related functions. Concurrently, ssGSEA was applied to assess the enrichment of immune cell signatures and immune-related functions. Statistical differences between subgroups were evaluated using the Wilcoxon rank-sum test. In addition, to assess intratumoral heterogeneity, the Mutant-Allele Tumor Heterogeneity (MATH) algorithm was utilized to calculate mutation-derived heterogeneity scores (version 2.4.6) (46). These MATH scores were subsequently incorporated into survival analyses to examine their association with patient prognosis and genomic instability.

To evaluate potential therapeutic implications of the LAGs-based risk model, drug response prediction was performed using the ‘pRRophetic’ R package (version 0.5.1) (47). The half-maximal inhibitory concentration (IC₅₀) values of commonly used chemotherapeutic and targeted agents were estimated for each patient. Differences in predicted drug sensitivity between high- and low-risk groups were compared using the Wilcoxon rank-sum test, with P<0.05 considered statistically significant. These analyses aimed to identify potential compounds with differential efficacy across risk subgroups, providing a foundation for personalized treatment strategies in CESC.

To further elucidate the causal relationship between LAGs prediction models and cervical cancer, two-sample MR analysis was employed to assess whether the expression levels of prognostic LAGs exert a causal effect on cervical cancer susceptibility. SNPs were selected as instrumental variables (IVs) if they demonstrated significant association with the exposure variable (P<5×10⁻⁸). To ensure variant independence and minimize linkage disequilibrium (LD) effects, SNPs were clustered using an LD threshold of R²<0.001 and a 10,000 kb clustering window. Variants with minor allele frequency <0.01 or F-statistic <10 were excluded to avoid weak instrument bias. The inverse variance weighting (IVW) method was employed as the primary analysis to estimate causal effects between LAGs and cervical cancer. Heterogeneity among IVs was assessed using the Cochran Q test, while potential horizontal pleiotropy and directional bias were assessed via MR-Egger regression. Sensitivity analyses, including stepwise exclusion tests, validated the robustness of causal estimates. All MR analyses were performed using the ‘TwoSampleMR’ and ‘MRPRESSO’ R packages (version 1.0) (48,49). Sensitivity analyses were conducted for the two-sample MR studies, including leave-one-out tests, pleiotropy analysis via the MR-Egger intercept test and heterogeneity analysis using Cochran's Q test (50). These analyses ensure the robustness and validity of the MR findings by identifying potential biases or heterogeneity in the genetic instruments used. Final causal estimates were visualized with forest plots and scatterplots to display the direction and magnitude of the association.

Human cervical squamous cell carcinoma cell lines SiHa (cat. no. HTB-35) and HeLa (cat. no. CRM-CCL-2) were obtained from the American Type Culture Collection and the human keratinocyte cell line HaCaT was obtained from The Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (SCSP-5091; Beijing, China). Although short tandem repeat analysis was not performed, the HaCaT cells were purchased from a certified national cell bank with traceable origins and their authenticity was confirmed by morphological observation and mycoplasma testing. All cell lines were cultured in DMEM (cat. no. KGL1201-500; Jiangsu KeyGen Biotech Co., Ltd.) containing 10% FBS (cat. no. C04001-500; Shanghai Xiaopeng Biotechnology Co., Ltd.) and 1% penicillin-streptomycin-amphotericin B (cat. no. P7630; Beijing Solare Technology Co., Ltd.) at 37°C in a humidified incubator with 5% CO₂.

SiHa and HeLa cells were seeded into 6-well plates and transfected with HMGN1-specific small interfering RNA (siRNA) or non-targeting control (NC) siRNA using siTran 2.0™ (cat. no. TT320002; Origene Technologies, Inc.) according to the manufacturer's instructions. The sequences and catalog numbers of the HMGN1-specific siRNAs and NC siRNA are provided in Table SIII (Tsigke). Cells were cultured to approximately 60–70% confluence prior to transfection. A final concentration of 50 nM siRNA was used for each transfection. For each well of a 6-well plate, the appropriate amount of siRNA was gently mixed with 2.4 µl siTran 2.0™ reagent and 100 µl transfection buffer, followed by incubation at room temperature for 15 min to allow complex formation. The complexes were then added dropwise to the cell culture medium. After incubating the cells at 37°C for 12 h, the medium was replaced with fresh growth medium. Cells were collected after an additional 12 or 24 h incubation. RNA and proteins were extracted separately for subsequent analysis.

Total RNA was extracted from SiHa and HeLa cells using the RNeasy™ Kit (cat no. R0027; Beyotime Biotechnology) according to the manufacturer's instructions. The extracted total RNA was then reverse transcribed into complementary DNA (cDNA) using SuperMix (R323-01; Vazyme Biotech Co., Ltd.). After cDNA synthesis, qPCR was performed on QuantStudio 5 instrumentation using SYBR^® Green PCR Master Mix (cat. no. Q311-02; Vazyme Biotech Co., Ltd.). The thermocycling conditions were as follows: An initial pre-denaturation at 95°C for 10 min, followed by 40 cycles of denaturation at 95°C for 10 sec, annealing at 58°C for 30 sec and extension at 72°C for 30 sec. GAPDH was used as the internal reference gene for normalization. Relative gene expression levels were calculated using the 2^−∆∆Cq method (51). The human-derived mRNA sequences of the 14 genes for the lactylation prediction model are shown in Table SIV. All mRNA sequences were provided by Beijing Qingke Biotechnology Co., Ltd.

Cells were collected and the medium was discarded. Total protein was extracted using RIPA lysis buffer (cat no. P0013B; Beyotime Biotechnology) supplemented with protease and phosphatase inhibitors (cat no. P1045; Beyotime Biotechnology). Samples were centrifuged at 12,000 × g for 15 min at 4°C and the supernatant was collected for subsequent analysis. Protein concentration was quantified using the BCA Protein Assay Kit (cat no. 20201ES86, Shanghai Yeasen Biotechnology Co., Ltd). Equal amounts (10 µg) of protein were separated by 10% SDS-PAGE and transferred onto PVDF membranes (cat no. IPVH00010; MerckKGaA). Membranes were incubated in Rapid Blocking Buffer (cat. no. P30500; Suzhou Xinsaimai Biotechnology Co., Ltd.) for 20 min at room temperature to minimize non-specific binding. After blocking, membranes were incubated overnight at 4°C with anti-HMGN1 (1:1,000; cat no. PK87656; Abmart, Inc.) and anti-β-actin (1:1,000; cat no. A17910; ABclonal, Inc.; internal control). After washing with 0.1% TBS-Tween, membranes were incubated with HRP-conjugated anti-rabbit secondary antibody (cat. no. ZB-2306; ZSGB-Bio, Inc.) for 1 h at room temperature. The membranes were detected using an ultra-sensitive ECL kit (cat no. P10100; Suzhou CellPro Biotechnology Co., Ltd.) and quantified using ImageJ software (version 1.0; National Institutes of Health).

Cell proliferation was assessed using the CCK-8 Assay Kit (cat. no. C0005; TargetMol Chemicals Inc.). Cells were seeded in 96-well plates and absorbance was measured at 24, 48 and 72 h after adhesion. To assess cell proliferation, 10 µl of CCK-8 reagent was added to each well and the cells were then incubated at 37°C and 5% CO₂ for 1 h. The absorbance at 450 nm was then measured using a microplate reader to determine cell viability.

Cell cycle distribution was analyzed using propidium iodide (PI) staining. Transfected SiHa and HeLa cells were seeded into 6-well plates and cultured under standard conditions. Cells were collected, washed twice with PBS (cat. no. C0221A; Beyotime Biotechnology) and fixed in pre-cooled 70% ethanol (cat. no. Y263010; Beyotime Biotechnology) at 4°C for 8 h. Fixed cells were centrifuged, washed and resuspended in staining buffer containing RNase A (cat. no. ST576; Beyotime Biotechnology) and 5 µg/ml PI (cat. no. G1021; Wuhan Servicebio Technology Co. Ltd.). After incubation for 20 min in the dark at room temperature, samples were analyzed using a FACSCalibur flow cytometer (BD Biosciences) and data were processed with CytExpert 2.0 software (Becton; Dickinson and Company).

All statistical analyses were conducted using R software (version 4.4.1; Posit Software, PBC) and GraphPad Prism (version 10.6.1; Dotmatics). All results represent the mean ± standard deviation of at least three independent experiments. Intergroup comparisons were performed using the two-tailed unpaired Student's t-test, while comparisons among multiple groups employed one-way ANOVA followed by Tukey's post hoc test. The relationship between two continuous variables was assessed using Wilcoxon rank-sum. To compare differences in Kaplan-Meier survival curves, the two-tailed log-rank test was applied. P<0.05 was considered to indicate a statistically significant difference, unless otherwise specified.

The GSE44001 dataset, comprising transcriptomic data from 300 patients with cervical cancer, included scRNA-seq data for a total of 22,702 cells. After rigorous quality control and filtering, cells were categorized into eight major clusters: Columnar and squamous epithelial cells, CD8⁺ T cells, CD4⁺ T cells, stromal cells, cancer cells, basal blood cells, and endothelial cells (Fig. 2A). The proportions of these cell populations are shown in a bar chart (Fig. 2B). To evaluate the activity of LAGs across different cell types, the ‘AddModuleScore’ function of the ‘Seurat’ package was applied. The resulting scores demonstrated that squamous and columnar epithelial cells exhibited higher LAG activity compared with other clusters (Fig. 2C and D), suggesting that LAGs play a key role in cervical cancer pathogenesis, as these epithelial cells are involved in tumor initiation and progression. Elevated LAG levels in these cells may be associated with their influence on cervical cancer progression. Based on the distribution of LAG activity, cells were stratified into high- and low-LAG groups, yielding 5,443 DEGs for subsequent analysis. To validate these findings, three GEO datasets (GSE122697, GSE146114 and GSE63514) were integrated after batch correction using the ComBat algorithm (Fig. S1A-D). From these datasets, disease-related DEGs were identified and the top 50 genes were visualized in a heatmap (Fig. 2E).

Next, WGCNA was performed to identify gene modules associated with cervical cancer. After removing outlier samples, hierarchical clustering and correlation analyses between normal and tumor groups were performed using the ssGSEA algorithm. With a merging threshold of 0.25, a total of nine co-expression modules were identified (Fig. 2F and G). Among them, the turquoise module (MEturquoise) revealed the strongest positive correlation with cervical cancer progression (r=0.57; P<0.001). The strong association between gene significance and module membership within this module (r=0.96; P<0.001; Fig. 2H) further confirmed its biological relevance. Genes within the turquoise module were thus selected for further analysis. Differential expression analysis between tumor and normal tissues revealed significant gene expression changes (adjusted P<0.05; log₂FC >0.5), which were visualized in a volcano plot (Fig. 2I). By intersecting turquoise module genes, DEGs and LAGs, a total of 43 overlapping genes were identified as LAGs (Fig. 2J). GO enrichment analysis revealed that these LAGs were predominantly involved in ‘monosaccharide metabolic process’, ‘hexose metabolic process’ and ‘pyridine-containing compound metabolic process’. They were also enriched in cellular components such as the ‘nuclear envelope’ and ‘cell-substrate junction’ and exhibited molecular functions associated with ‘tubulin binding’ (Fig. 2K; Table SV).

To elucidate the prognostic relevance of LAGs in CESC, multiple machine learning algorithms were employed to construct and evaluate prognostic models based on 43 differentially expressed LAGs identified in prior analyses. TCGA-CESC cohort was used as the training dataset, whereas GSE30760 served as the independent validation cohort. A total of 101 predictive models were generated via cross-validation and their predictive efficacy was systematically compared using the C-index across both training and validation datasets (Fig. 3A). After excluding models exhibiting overfitting tendencies, the StepCox (Backward) model demonstrated the most stable and robust performance and was therefore selected as the optimal prognostic model. This refined model incorporated 14 feature genes with significant prognostic relevance. For each patient, a risk score was calculated by summing the expression levels of these 14 genes weighted by their respective Cox regression coefficients (Fig. 3B). Based on the median risk score, patients were stratified into high-risk and low-risk groups. Kaplan-Meier survival analysis revealed that patients in the high-risk group exhibited significantly worse OS compared with those in the low-risk group (P<0.001; Fig. 3C). Consistent findings were observed for DSS and PFI, with low-risk patients demonstrating notably increased clinical outcomes (P<0.0001; Fig. 3D and E). To further assess the discriminatory capability of the model, time-dependent receiver operating characteristic (ROC) curve analyses were conducted. In TCGA training cohort, the model achieved AUC values of 0.80, 0.78 and 0.77 for 1-, 3- and 5-year survival, respectively (Fig. 3F). Similar results were obtained in the GSE30760 validation cohort, yielding AUC values of 0.76, 0.67 and 0.67, indicating moderate generalization ability and demonstrating the acceptable performance of the model in short-term survival prediction (Fig. 3G). Furthermore, external validation using the GSE7803 dataset further confirmed the predictive reliability of the 14-gene model, with AUC values >0.5 for all included genes except PFKM. However, due to sample limitations, the predictive capability of genes other than PAK2 and TCOF1 may be limited (Fig. 3H). These findings indicate that LAG-based predictive models may be able to evaluate the prognosis of patients with cervical cancer and provide important reference for clinical stratification management.

To further evaluate the clinical applicability and prognostic independence of the LAGs risk model in cervical cancer, the present study examined the relationship between the LAGs-derived risk score and diverse clinicopathological characteristics using stratified and Cox regression analyses. This comprehensive evaluation aimed to assess the ability of the model to predict patient outcomes across different clinical subgroups. Patients were stratified according to clinical parameters, including age (>60 vs. ≤60 years), T stage (T1; T2; T3 and T4), N stage (N0, N1 and NX), M stage (M0, M1 and MX) and tumor stage (I; II; III and IV), as illustrated in the clinical distribution plot (Fig. 4A and B). Notably, patients in the stage M1, III–IV and T3-4 subgroups exhibited significantly higher risk scores compared with those patients in the M0, I–II and T1-2 subgroups, respectively (Fig. 4C and D). These findings indicate a strong association between elevated LAGs risk scores and advanced tumor stage as well as unfavorable prognosis in CESC.

ROC curve analysis further confirmed the predictive utility of the LAGs-based model, with an AUC value of 0.76 in distinguishing the prognosis of patients with different M stages, indicating moderate predictive capability (Fig. 4E). In addition, stratified analyses by age and tumor stage demonstrated that high-risk patients consistently exhibited worse survival outcomes (Fig. 4F). Kaplan-Meier analyses of the 14 model genes using the GEPIA2 database that low expression levels of transketolase (TKT) and phosphoglycerate kinase 1 (PGK1) were significantly associated with improved prognosis (P<0.05), while no significant associations were observed for other genes (Fig. S2A). To further validate the prognostic independence of the model, both univariate and multivariate Cox regression analyses were performed for OS, DFI, PFI and DSS in the training cohort. The results demonstrated that LAGs were independent prognostic factors (all P<0.05; Figs. 4G and H; and S2B and C). Lastly, a nomogram integrating the risk score with key clinical parameters was constructed to quantitatively predict survival outcomes in patients with CESC (Fig. 4I). The nomogram analysis confirmed that the risk score significantly contributed to the highest prognostic weight (P<0.001), underscoring its clinical utility as a robust and independent predictor of cervical cancer prognosis.

To further elucidate the molecular mechanisms that associate LAGs with the prognosis of cervical cancer, functional annotation analyses were performed across the identified LAGs-derived molecular subtypes. In the low-risk group, Gene Set Variation Analysis (GSVA) revealed significant enrichment of biological pathways associated with ‘receptor ligand activity, ‘signaling receptor activator activity’ and ‘cell adhesion molecule binding’ (Fig. 5A). KEGG pathway enrichment further demonstrated that the most enriched pathways were predominantly involved in energy metabolism, including ‘central carbon metabolism in cancer’, ‘biosynthesis of amino acids’ and ‘glycolysis/gluconeogenesis’ (Fig. 5B). Comparative GSVA between risk subgroups (Table SVI) indicated distinct molecular signatures: The high-risk group exhibited stronger activation of oncogenic and proliferative pathways, such as KRAS signaling downregulation, whereas the low-risk group displayed higher activity in immune- and stress-related pathways, including ‘hypoxia’, ‘unfolded protein response’ and ‘angiogenesis’ (Fig. 5C). Correlation analysis further confirmed that LAGs risk scores were significantly associated with these key signaling pathways (Fig. 5H). To assess the prognostic relevance of these pathways, Kaplan-Meier survival analyses were performed. Pathways negatively correlated with LAGs expression, such as cholesterol homeostasis, unfolded protein response, hypoxia and angiogenesis, were associated with favorable survival outcomes (Fig. 5D-G). By contrast, KRAS signaling downregulation, which was positively correlated with LAGs, was associated with poor prognosis (Fig. S2D). Collectively, these findings suggest that differential pathway activation between LAGs-defined subgroups reflects distinct metabolic and signaling landscapes, thereby elucidating potential molecular mechanisms by which LAGs contribute to cervical cancer progression.

To investigate the role of LAG in the tumor immune microenvironment at the single-cell transcriptome level, the expression profiles of LAG prognostic genes were analyzed across distinct cell populations (Fig. 6A). Based on the LAGs prognostic model, tumor cells were categorized into high- and low-risk groups, followed by differential interaction analysis to examine intercellular communication patterns within the TME. Distinct cell-cell communication networks were observed between the two risk groups (Fig. 6B and C). Notably, intercellular signaling occurred among nearly all cell types (Fig. S2E), with particularly strong interactions detected between columnar and squamous epithelial cells. Functional annotation of ligand-receptor signaling revealed that the collagen, laminin and Notch pathways were predominantly enriched in highly cancerous cells, squamous and columnar epithelial cells, CD4⁺ T cells and endothelial cells (Fig. 6D-F). These intercellular interactions are likely to regulate cell adhesion, migration and differentiation, thereby contributing to tumor progression and microenvironmental remodeling.

Since intratumoral heterogeneity is a major driver of tumorigenesis, progression, metastasis, therapeutic resistance and immune evasion, the present study next evaluated heterogeneity using the MATH score. Tumors with high LAGs risk exhibited markedly elevated MATH scores, indicating greater genomic instability (Fig. 6G). Furthermore, integrative analysis combining LAGs and MATH metrics revealed that patients in the high-risk and high-MATH score subgroup experienced significantly worse survival outcomes compared with those in the low-risk and low-MATH score subgroup (Fig. 6H; P<0.001). Collectively, these findings suggest that LAG activity is closely associated with alterations in cellular communication, tumor heterogeneity and the immune microenvironment, potentially influencing cervical cancer progression and patient prognosis.

Immune infiltration serves a key role in tumor initiation, progression and therapeutic response, profoundly influencing the clinical outcomes of patients with cancer. To characterize the immune landscape of cervical cancer, the present study calculated the immune, stromal and ESTIMATE scores for each sample (Fig. 7A). Although the high-risk group exhibited slightly elevated immune and stromal scores compared with the low-risk group, these differences did not reach statistical significance. Subsequently, ssGSEA was performed to evaluate immune pathway activities. The results demonstrated significant enrichment of the ‘antigen processing and presentation’ pathway and the ‘NOD-like receptor signaling’ pathway in the high-risk group (Fig. 7B). To further elucidate the immune cell composition, the abundance of tumor-infiltrating immune cells was quantified across samples. Activated B cells, eosinophils and immature B cells were predominantly enriched in the low-risk group, displaying significant intergroup differences (Fig. 7C). These findings were independently validated using the CIBERSORT algorithm, which yielded consistent results (Fig. 7D). Spearman correlation analysis identified nine immune cell subtypes significantly associated with LAGs risk scores (Fig. 7E). Furthermore, correlation analyses between LAGs prognostic genes and specific immune cell populations revealed that AT-rich interaction domain 3B, phosphoglycerate kinase 1 (PGK1) and protein phosphatase, Mg²⁺/Mn²⁺ dependent 1G (PPM1G) were positively correlated with M0 macrophages and activated mast cells and indicated widespread associations with monocyte-related gene expression (Figs. 7F and S2F). Lastly, the relationship between TME infiltration and OS was examined and six immune cell subtypes were significantly associated with patient prognosis (log-rank test; P<0.05; Fig. S2G-L), including M0 and M1 macrophages, resting mast cells, CD8, γΔ and regulatory T cells. Collectively, these findings indicate that the LAGs risk signature is closely associated with alterations in the immune landscape of cervical cancer, potentially influencing both tumor immunity and clinical outcomes.

To further assess the clinical utility of the LAGs-based risk model in guiding treatment strategies, the present study explored its association with immunotherapy response and drug sensitivity. Using the IMvigor210 cohort, the present study evaluated the predictive value of the LAGs risk score for immunotherapy efficacy. The distribution of clinical responses-including complete remission (CR), partial remission (PR), stable disease (SD) and progressive disease (PD)-was compared between risk subgroups using χ² tests. The high-risk group exhibited a markedly higher proportion of PD/SD cases, whereas CR/PR responses were significantly more frequent in the low-risk group (Fig. 8A). Consistently, patients achieving CR/PR demonstrated significantly lower risk scores compared with those with SD/PD, suggesting that lower LAGs scores are associated with more favorable immunotherapeutic outcomes (Fig. 8B and C). To evaluate chemotherapeutic drug sensitivity, the IC₅₀ values of various agents were compared between high- and low-risk subgroups. The risk score exhibited a significant positive correlation with the IC₅₀ values of paclitaxel (R=0.16; P<0.05) and rapamycin (R=0.26; P<0.05) (Fig. 8D and E) and a negative correlation with midostaurin (Fig. 8F; R=−0.32; P<0.05). Specifically, paclitaxel and rapamycin displayed lower IC₅₀ values in the low-risk group, indicating higher predicted sensitivity, while midostaurin demonstrated higher IC₅₀ values in the high-risk group, suggesting potential efficacy in these patients (Fig. 8G-I).

Furthermore, to validate the expression of LAGs in vitro, the present study measured mRNA expression levels of the 14 LAGs included in the prognostic model using RT-qPCR in HaCaT epithelial cells and cervical squamous carcinoma cell lines (SiHa and HeLa). The results revealed that the majority of genes, including G6PD, HMGN1, ZYX, PAK2, TCOF1 and TKT, were significantly upregulated in cervical cancer cells compared with normal epithelial cells (Fig. 8J). Therefore, these findings demonstrate that the LAGs-based prognostic model not only stratifies patients by survival risk but also provides notable insights into immunotherapy responsiveness and chemotherapeutic sensitivity, thereby offering a potential framework for personalized treatment optimization in cervical cancer in the future.

To further investigate the relationship between the LAG prognostic model and cervical cancer and to determine whether this relationship is causal or merely correlative, the present study performed a two-sample MR analysis. In this analysis, prognostic genes were used as exposure factors and cervical cancer was the outcome. The present study was unable to obtain identifiers for several genes, including glucose-6-phosphate dehydrogenase, proteasome subunit α type 7, treacle ribosome biogenesis factor 1, PGK1, PPM1G and PFKM, as certain genes lacked sufficient SNPs, leading to their exclusion from the analysis. Using the IVW method (Table I), the present study identified HMGN1 as a protective factor, with an odds ratio of 0.81 (95% CI, 0.70–0.94; P=0.0051; Fig. 9A and D). However, no significant causal relationship was observed between zyxin (ZYX) and TKT (P>0.05; Fig. 9B and C, E and F).

To assess the robustness of the present analysis, funnel plots were constructed using the β coefficients and SEs for each IV associated with the three genes studied (Fig. 9G-I). The distributions of these plots revealed certain asymmetry, possibly due to the limited number of IVs, which may introduce heterogeneity and potential bias. To account for this, the present study performed heterogeneity tests, all of which yielded P>0.05, indicating that heterogeneity did not significantly affect the present study results (Table SVII). Furthermore, the present study conducted leave-one-out analyses to assess the stability of the present study findings. The results remained consistent even after excluding individual SNPs, further supporting the robustness of the associations (Fig. 9J-L). The present study also performed a test for horizontal pleiotropy on HMGN1, with P>0.05, indicating no significant evidence of horizontal pleiotropy influencing the present study results (Table SVIII). The present study results provide strong evidence for a causal relationship between HMGN1 and cervical cancer, while no significant MR associations were observed for the other analyzed genes. These findings highlight the potential role of HMGN1 in the pathogenesis of cervical cancer.

To further investigate the role of HMGN1 in cervical cancer, the present study designed HMGN1-specific siRNA sequences using siRNA technology and transfected them into SiHa and HeLa cells. Subsequently, 24 h after transfection, the present study assessed the knockdown efficiency of HMGN1 using RT-qPCR, which confirmed that HMGN1 was successfully silenced in both cell lines (Fig. 10A). Subsequent western blotting was performed to assess the protein expression of HMGN1, further validating the effectiveness of the knockdown (Fig. 10B). Functional assays revealed that inhibition of HMGN1 significantly suppressed the proliferation of SiHa and HeLa cells (Fig. 10C), reduced the proportion of S-phase cells and prolonged the G₂ phase of the cell cycle (Fig. 10D). Therefore, these findings suggest that HMGN1 may be a risk factor for CESC.