Integrated transcriptome analysis and machine learning to construct a homeostatic model of acetylation for bladder cancer and validate the key gene CES1

Introduction

Bladder cancer (BLCA), among the most prevalent cancers that impact the urinary system, has a rising incidence and mortality rate annually, posing a substantial challenge to public health (1). At present, intravesical chemotherapy and transurethral bladder lesion excision are the standard treatments for non-muscle-invasive BLCA (2). For muscle-invasive BLCA, radical cystectomy remains the mainstay of treatment (3). Although therapeutic strategies for BLCA have continued to evolve in recent years (4), particularly with advances in molecular targeted therapies and immune checkpoint inhibitors, these approaches have not yet effectively halted disease progression or significantly improved overall survival (OS) (5-7). Therefore, identifying reliable prognostic factors is crucial for predicting progression, choosing treatment, and ameliorating OS. Meanwhile, further investigations into the biological mechanisms of BLCA and new treatment targets remain urgently needed.

Post-translational modifications profoundly influence protein function and are essential in nearly all cellular biological processes (BPs). Protein acetylation (PA) is a widespread post-translational modification entailing the attachment of an acetyl moiety, derived from coenzyme A (CoA), onto specific lysine residues within target proteins, a process catalyzed by acetyltransferases (8). Following acetylation, changes in protein conformation and charge properties occur, thereby facilitating the regulation of diverse cellular functions. PA influences many cellular processes, including regulating gene transcription, modulating metabolic enzyme activity, controlling protein degradation and localization, and maintaining cellular homeostasis (9-11).

Highlighted PA is a critical factor in cancer, as it modulates tumor-suppressive and carcinogenic signaling pathways (12). Aberrant acetylation disrupts cellular homeostasis, thereby influencing tumor initiation, progression, and metastasis (13). For instance, it can modulate transcription factors and metabolic enzymes, ultimately affecting metabolism and the tumor microenvironment (TME). Therefore, targeting PA pathways is promising in cancer therapy (14). Intracellular CoA is closely linked to acetylation processes, and CoA alterations directly influence PA, thus affecting gene transcription, DNA replication, and other linked processes (15,16).

PA is important in the biological behavior of various malignancies. In renal cancer, SFMBT2 is stabilized through auto-acetylation, thereby suppressing tumor growth (17). In colorectal cancer, bacterial infection promotes tumorigenesis by regulating CDC42 acetylation (18). In breast cancer, acetyl-CoA carboxylase 1-dependent PA influences metastasis and recurrence (19). In cholangiocarcinoma, enhanced PA reduces cellular proliferation and migration (20). However, systematic investigations of PA and acetylation homeostasis in BLCA remain limited.

This study first established and validated a novel acetylation homeostatic model by integrating transcriptome analysis and machine learning for BLCA. The expression patterns of key genes in the signature were further explored at the single-cell level. Enrichment analysis (EA), immune infiltration analysis (IIA), and drug sensitivity analysis (DSA) were performed to reveal molecular differences between high‑and low‑risk groups. Furthermore, we conducted in vitro functional experiments to validate CES1 as the key gene in regulating PA homeostasis and malignant phenotypes of BLCA cells. Overall, our findings elucidate the vital functions of acetylation-related genes in BLCA progression and provide promising biomarkers for prognosis evaluation and targeted therapy of BLCA. We present this article in accordance with the MDAR and TRIPOD reporting checklists (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0435/rc).

Methods

Data collection and processing

RNA sequencing (RNA-seq) and clinical information came from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). The study included RNA-seq datasets TCGA-BLCA, containing 431 clinical samples, and GSE32894, based on the GPL6947 platform and comprising 308 clinical samples, and GSE13507, based on the GPL6102 platform, consisting of 165 clinical samples. In addition, two single-cell RNA-seq datasets were retrieved: GSE130001 generated on the GPL16791 platform (two BLCA samples) and GSE135337 established on the GPL24676 platform (eight BLCA samples). Transcripts per million (TPM) values were extracted from each dataset for subsequent analyses. Duplicated gene entries were discarded, and genes showing no expression (TPM =0) in all samples were excluded from subsequent analysis.

PA-linked genes were obtained from the Human Gene Database (https://www.genecards.org/). Genes with an acetylation relevance score <2 were removed, and 1,657 PA-linked genes were encompassed. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Acetylation-linked DEG identification

Differential expression analysis (DEA) of the TCGA-BLCA dataset was performed via limma. DEGs across bladder tumor and adjacent normal tissues were identified per the criteria P<0.05 and |logFC| ≥1. Acetylation-linked DEGs were obtained by intersecting the acetylation-related gene set with the identified DEGs.

EA

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway EAs were carried out utilizing clusterProfiler.

Protein-protein interaction (PPI) network construction

The PPI network was constructed leveraging STRING (https://string-db.org/). Acetylation-linked DEGs were uploaded to STRING for the interaction network with a confidence score threshold of 0.7. Cytoscape was employed for PPI network visualization.

Acetylation homeostatic model construction and validation

Univariate Cox regression analysis (CRA) was enabled by survminer and survival to preliminarily identify prognosis-related genes, with P<0.05 signifying statistical significance. Tumor samples in the TCGA-BLCA dataset were randomized into a training set (70%) and a validation set (30%). Machine learning-based least absolute shrinkage and selection operator (LASSO) CRA in the former set was enabled via glmnet. The random seed was set to 2025. A penalized function was fitted to construct the LASSO-Cox model, and a cross-validation LASSO (CV-LASSO) plot was generated to determine the optimal minimum λ value and the corresponding model genes. Multivariate CRA was subsequently carried out. Genes with P<0.05 were noted as significant prognostic risk genes.

A risk-score formula was constructed: Riskscore=∑i=15Expi×Coefi. Based on the median score, tumor samples were split into high- and low-risk cohorts. KM survival curves were obtained via survminer for OS comparison across groups. Time-dependent receiver operating characteristic (ROC) curves and the area under the curve (AUC) were derived at 1, 3, and 5 years via timeROC to rate the model’s accuracy and predictive performance. We further validated the acetylation homeostatic model performance in an internal training set (70%) and an internal validation set (30%) from TCGA-BLCA, as well as in an independent external validation dataset GSE32894.

Assessment of independent prognostic predictors

Univariate CRA was performed to assess the effects of risk scores and clinical characteristics, including age, gender, and tumor-node-metastasis (TNM) stage, on prognosis. Multivariate CRA was then conducted to identify independent prognostic factors of the acetylation homeostatic model.

Nomogram construction and analysis

Nomograms were constructed using the rms package. One nomogram integrated the risk score and clinical characteristics (age and TNM stage) with survival status to predict the 1-, 3-, and 5-year OS. The predictive value of this model was validated using calibration curves. The other nomogram integrated the expression levels of the prognostic model genes with survival status for the prediction of 1-, 3-, and 5-year OS, and its predictive performance was also verified by calibration curves.

scRNA-seq data processing

scRNA-seq data were processed via Seurat and SingleR. A gene expression matrix was generated utilizing Read10X, and a Seurat object was created with CreateSeuratObject. To ensure high-quality cells, only genes expressed in three or more individual cells were kept. Cells were excluded if they had fewer than 200 or over 4,000 detected genes, or if mitochondrial gene expression exceeded 25%. The top 2,000 highly variable genes were identified employing the vst approach. Principal component analysis (PCA) was carried out for dimensionality reduction, followed by cell clustering using FindNeighbors and FindClusters (resolution =0.5). t-SNE was used for dimensionality reduction and cell cluster visualization. Cell type reference data came from the human primary cell atlas in Celldex. SingleR was utilized for annotating cell types and examining key prognostic model gene expression across cell types.

EA of DEGs across risk groups

To explore differences in signaling pathways across risk cohorts, DEA was performed using limma. Genes with P<0.05 and |logFC| ≥1 were identified as DEGs across cohorts. GO and KEGG EAs were enabled via clusterProfiler. Leveraging GO and KEGG gene sets from GSEA (https://www.gsea-msigdb.org/), gene set enrichment analysis (GSEA) was performed by clusterProfiler and enrichplot. Enrichment scores were quantified, allowing identification of signaling pathways enriched in two risk cohorts.

IIA

IIA was enabled by CIBERSORT for comparing infiltration levels of 22 immune cell sorts across risk cohorts. ESTIMATE was used to derive ESTIMATE, immune, as well as stromal scores for tumor microenvironment (TME) evaluation.

DSA

DSA was completed via pRRophetic. The half-maximal inhibitory concentration (IC₅₀) of commonly used drugs was calculated for different risk groups, and differences in IC₅₀ were used to compare chemotherapeutic sensitivity across groups. P<0.01 signified statistical significance.

Immunohistochemistry (IHC) and prediction of modification sites

IHC images of bladder tissues came from the Human Protein Atlas (https://www.proteinatlas.org/). In addition, the predicted PA modification sites of CES1 were retrieved from PhosphoSitePlus (https://www.phosphosite.org/).

Cell culture and transfection

CES1, the key gene with the lowest P value in multivariate CRA, was selected as the key gene for experimental verification. Human BLCA cell lines, T24 and UMUC3, were utilized. Acquired from Procell Life Science & Technology Co., Ltd. (Wuhan, China), both were confirmed to be mycoplasma-negative and were identity-verified by STR genotyping. T24 was cultured in complete DMEM, while UMUC3 was cultured in complete MEM at 37 ℃ with 5% CO₂. After trypsinization, log-phase cells were seeded in plates with 6 wells. Transfection was performed with Lipofectamine 3000 when cells attained approximately 80% confluence. Both were transfected with si-CES1#1 or si-CES1#2 to knock down CES1 expression, while cells transfected with si-control served as negative controls. CES1 expression levels were measured 48 h after transfection. The siRNA sequences were: si-CES1#1 (ID:1066 siRNA-1693): 5'-CCAAGAAGGCAGUGGAGAA-3', si-CES1#2 (ID:1066 siRNA-1516): 5'-GCAAGAUGGUGAUGAAAUU-3', si-control: 5'-UUCUCCGAACGUGUCACGU-3'.

Quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was extracted utilizing TRNzol Universal Extraction Reagent (DP424, Tiangen, Beijing, China). qRT-PCR was performed employing 2× SYBR Green qPCR Master Mix (High ROX) (G3322, Wuhan, China). With β-Actin as the internal normalization control, relative gene expression was derived via the 2^−ΔΔCT approach. The primer sequences were: CES1_F: ATCCACTCTCCGAAGGGCAACT and CES1_R: GACAGTGTCGTCTGTTCCTCCT.

Western blot (WB) analysis

T24 and UMUC3 cells were collected 48 hours after transfection and lysed on ice with RIPA buffer containing phosphatase and protease inhibitors to obtain total cellular protein. The proteins were subsequently resolved via SDS-PAGE and electro-transferred onto PVDF membranes. Following a 1-hour blocking at room temperature utilizing 5% non-fat milk, the membranes were probed overnight at 4 ℃ with CES1 monoclonal antibody(Manufacturer: Proteintech, Catalog No.: 67079-1-Ig, Clone No.: 1C11C1) diluted at 1:5,000 and β-actin (Manufacturer: Proteintech, Catalog No.: 66009-1-Ig, Clone No.: 2D4H5) diluted at 1:10,000. After washing thrice with TBS-T, the membranes were incubated for one hour at room temperature with an HRP-conjugated goat anti-mouse IgG (H+L) secondary antibody. Protein bands were visualized via ImageQuant LAS 500, and an enhanced chemiluminescence (ECL) detection kit, and quantitative analysis was performed via ImageJ.

Wound healing assay

Cell migration was examined via a wound healing assay. T24 and UMUC3 cells were plated in 6-well plates at 30% confluence with 2 mL of culture medium and allowed to adhere overnight. Twenty-four hours post siRNA transfection, cells were scratched with a 10-µL pipette tip to create a linear wound, followed by phosphate-buffered saline (PBS) washes to clear debris and the addition of fresh medium before returning the culture to the incubator. The same visual field was imaged immediately after scratching (0 h) and again after 24 hours. Cell migration was analyzed utilizing ImageJ.

Cell Counting Kit-8 (CCK-8) assay

Cell proliferation was assessed utilizing the CCK-8 assay. T24 and UMUC3 were transfected and plated in 96-well plates with 4,000 cells per well. After 24, 48, and 72 hours, 10 µL of CCK-8 solution was added. The optical density (OD) at 450 nm was recorded using a microplate reader.

EdU proliferation assay

Cell proliferation was further examined via an EdU (5-ethynyl-2'-deoxyuridine) incorporation assay. Cells were plated in 96‑well plates (3×10³ cells/well) and adhered for 24 h. An EdU solution (BeyoClick™ EdU Cell Proliferation Kit with AF555) was added and incubated for 2 h. Thereafter, cells were fixed with 4% paraformaldehyde and permeabilized using PBS containing 0.5% Triton X‑100 for 10 min. After incubation with the appropriate staining solution for approximately 30 min, nuclei were counterstained with 4',6-diamidino-2-phenylindole (DAPI). Proliferating cells displayed red or green fluorescence, while nuclei appeared blue under DAPI staining. Fluorescence images were acquired via a fluorescence microscope, with the percentage of proliferating cells determined via ImageJ.

For quantitative analysis of EdU-positive cells, at least five randomly selected, non-overlapping microscopic fields were captured from each biological replicate using fluorescence microscopy. Quantitative statistical analyses based on field counting are presented in the results.

Colony formation assay

A colony formation assay was conducted for cellular proliferation assessment. Cells were plated in 6-well plates before siRNA transfection. Following transfection, the cells were cultured for 10 days. The resulting colonies were washed with PBS, fixed using 4% paraformaldehyde, and stained with 0.1% crystal violet. Images of colonies were obtained using an HP Scanjet G4050 scanner for subsequent analysis.

CoA content measurement

Intracellular CoA levels were measured via a CoA Assay Kit (Beijing Solarbio Science & Technology Co., Ltd.).

Statistical analysis

Statistical analyses and data visualizations were enabled by R 4.5.0. All in vitro experiments were performed using three independent biological replicates, with three technical replicates included for each biological replicate. Data are presented as the mean ± standard deviation (SD). Statistical significance between two groups was analyzed using an unpaired two-tailed Student’s t-test, and comparisons among multiple groups were performed by one-way analysis of variance (ANOVA). P<0.05 signified statistical significance (*P<0.05; **P<0.01; ***P<0.001; ****P<0.0001).

Results

Acetylation-related DEG analysis

First, RNA-seq data from 412 BLCA tumors and 19 matched adjacent normal tissue samples from TCGA were subjected to DEA, resulting in the identification of 1,627 DEGs (Figure 1A). Subsequently, these 1,627 DEGs were intersected with 1,657 acetylation-related genes retrieved from the Human Gene Database, yielding 163 candidate acetylation-related DEGs (Figure 1B). Their importance in the acetylation regulatory network of tumors is noteworthy. A heatmap presented the expression levels of the 163 candidate acetylation-related DEGs, illustrating their expression heterogeneity (Figure 1C). PPI network analysis indicated that MYC, MCM4, CDK1, CCNA2, H2BC21, TOP2A, CHEK1, CDK2, CENPA, JUN, and CDC20 may function as hub genes (Figure 1D).

Figure 1 Analysis of acetylation-related DEGs. (A) Volcano plot of DEA in BLCA. (B) Venn diagram showing 163 genes at the intersection of acetylation-related genes and DEGs. (C) Heatmap showing the expression levels of intersecting genes in normal tissues and BLCA tissues. (D) PPI network of the intersecting genes. (E,F) GO and KEGG EAs of the intersecting genes. BLCA, bladder cancer; BP, biological process; CC, cellular component; DEA, differential expression analysis; DEGs, differentially expressed genes; EA, enrichment analysis; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; PPI, protein-protein interaction.

To learn more about the biological characteristics of these 163 candidate genes, GO and KEGG EAs were carried out. The top enriched BPs were nucleosome assembly and negative regulation of the cell cycle; the most significantly enriched cellular component was the minichromosome maintenance protein complex; and the top enriched molecular functions included N-acetylgalactosaminyltransferase activity and polypeptide N-acetylgalactosaminyltransferase activity (Figure 1E). The most significantly enriched were the p53 signaling pathway and cellular senescence (Figure 1F). These pathways may regulate tumor cell proliferation, invasion, and migration, and may be critical in tumor initiation, progression, and prognosis.

Acetylation homeostatic model construction

Univariate CRA was first performed to preliminarily screen acetylation-related genes linked to prognosis. 39 acetylation-related genes were markedly related to OS (P<0.05) (Figure 2A). Subsequently, Machine learning-based LASSO-CRA was applied to reduce the number of feature genes and eliminate redundancy among these 39 candidates. After preprocessing the gene expression data, a LASSO-Cox regression model was established. To maximize predictive performance while controlling model complexity, tenfold cross-validation was carried out for the optimal minimum λ value (0.0296) and its corresponding model coefficients (Table S1). Genes with non-zero coefficients were kept, and 15 acetylation-linked genes were found (Figure 2B,2C). They were then subjected to multivariate CRA, which identified CES1, ATAD3A, H1F2, FASN, and GALNT14 as the final genes of the acetylation homeostatic model. (P<0.05) (Figure 2D). A nomogram incorporating the five prognostic genes was developed for 1-, 3-, and 5-year OS prediction (Figure 2E). Calibration curve analysis demonstrated good agreement across forecast and real OS at 1, 3, and 5 years, indicating reliable predictive performance for survival prognosis (Figure 2F).

Figure 2 Construction of the acetylation homeostatic model based on transcriptome and machine learning. (A) Univariate CRA initially identified 39 acetylation-related genes linked to prognosis. (B,C) LASSO-CRA further screened 15 candidate genes. (D) Multivariate CRA identified five prognostic genes with significant associations. (E) Nomogram constructed based on the prognostic genes to predict 1-, 3-, and 5-year OS. (F) Calibration curves for the nomogram. *, P<0.05; **, P<0.01 ; ***, P<0.001. AIC, Akaike information criterion; CI, confidence interval; CRA, Cox regression analysis; LASSO, least absolute shrinkage and selection operator; OS, overall survival.

Assessment of independent prognostic predictors of the acetylation homeostatic model

Univariate CRA showed that the risk score [P<0.001, hazard ratio (HR) =1.04, 95% confidence interval (CI): 1.03–1.05], clinical stage (P<0.001), pT stage (P=0.001, HR =1.544, 95% CI: 1.19–2.00), and age (P<0.001, HR =1.037, 95% CI: 1.02–1.06) were significantly associated with prognosis (Figure 3A). We further performed multivariate CRA, which revealed that the risk score (P<0.001, HR =1.035, 95% CI: 1.023–1.048), pN stage (P=0.041, HR =1.191, 95% CI: 1.007–1.410), and age (P=0.03, HR =1.023, 95% CI: 1.003–1.044) were independently associated with prognosis (Figure 3B). Among these factors, the risk score exhibited the highest significance as an independent prognostic indicator. A nomogram was constructed based on the risk score and clinical characteristics (age, TNM stage) to predict 1-, 3-, and 5-year OS (Figure 3C). Calibration curve analysis demonstrated good consistency between the predicted and actual OS at the 1-, 3-, and 5-year time points (Figure 3D). Collectively, these findings indicate that the prognostic model based on the risk score exhibits reliable predictive performance and can be used as a valid tool for clinical prognosis evaluation.

Figure 3 Independent prognostic value of the acetylation homeostatic model. (A) Forest plot of univariable Cox regression analysis. (B) Forest plot of multivariate Cox regression analysis. (C) Nomogram of predicting 1-, 3-, and 5-year overall survival rates. (D) Calibration curve. **, P<0.01 ; ***, P<0.001. C.L., confidence limit; CI, confidence interval; HR, hazard ratio; M, metastasis; N, node; T, tumor.

Acetylation homeostatic model validation and evaluation

Samples from TCGA-BLCA, GSE32894, GSE13507, as well as the internal training and validation subsets of TCGA-BLCA, were split into high- and low-risk cohorts. KM analysis revealed markedly poorer OS in the high-risk cohort than in the low-risk cohort. Therefore, the prognostic model effectively predicts patient outcomes (Figure 4A-4C, Figure S1A,S1B). Subsequently, the distributions of key gene expression levels, survival time, status, and risk scores were analyzed. In the TCGA-BLCA cohort, increasing risk scores were linked to shorter survival times and a higher number of deaths, and patients with higher expression levels of key genes tended to have higher risk scores and mortality rates (Figure 4D). Similar trends were observed in the GSE32894 cohort, the GSE13507 cohort, and the internal training and validation subsets of TCGA-BLCA (Figure 4E,4F, Figure S1C,S1D). Overall, these findings indicate that the prognostic model demonstrates good accuracy and stability across different datasets. The model’s predictive accuracy was examined via time-dependent ROC curves and ROC curves based on gene risk scores. In the TCGA-BLCA cohort, 1-, 3-, and 5-year AUCs were 0.681, 0.672, and 0.685, respectively, with a gene risk score AUC of 0.656 (Figure 4G,4H). In the GSE32894 cohort, the corresponding AUCs were 0.796, 0.774, and 0.775, alongside a gene risk score AUC of 0.765 (Figure 4I,4J). For the GSE13507 cohort, the 1-, 3-, and 5-year AUCs were 0.721, 0.687, and 0.663, respectively. The AUC of the risk score was 0.629 (Figure 4K,4L). For the internal training and validation subsets of TCGA-BLCA, the corresponding AUCs were 0.696, 0.686, 0.704, and 0.655, 0.641, 0.659, respectively, while the gene risk score AUC values were 0.652 and 0.666, respectively (Figure S1E-S1H). Therefore, the acetylation homeostatic model has good predictive accuracy and can effectively discriminate survival outcomes among BLCA patients.

Figure 4 Validation of the acetylation homeostatic model in multiple cohorts. (A) Kaplan-Meier OS curves for high- and low-risk groups in the TCGA-BLCA cohort. (B) Kaplan-Meier OS curves for high- and low-risk groups in the GSE32894 cohort. (C) Kaplan-Meier OS curves comparing the high- and low-risk groups in the GSE13507 cohort. (D) Risk score distribution, survival status scatter plot, and expression levels of prognostic genes in the TCGA-BLCA cohort. (E) Risk score distribution, survival status scatter plot, and expression levels of prognostic genes in the GSE32894 cohort. (F) Distribution of risk scores, survival status, and heatmap showing the expression of prognostic genes in the GSE13507 cohort. (G) ROC curves for 1-, 3-, and 5-year OS in the TCGA-BLCA cohort. (H) ROC curve of the risk score in the TCGA-BLCA cohort. (I) ROC curves for 1-, 3-, and 5-year OS in the GSE32894 cohort. (J) ROC curve of the risk score in the GSE32894 cohort. (K) Time-dependent ROC curves for predicting 1-, 3-, and 5-year OS in the GSE13507 cohort. (L) ROC curve evaluating the predictive performance of the risk score in the GSE13507 cohort. AUC, area under the curve; BLCA, bladder cancer; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

EA of DEGs between high- and low-risk groups

DEGs across risk groups were identified and subjected to GO and KEGG EAs, as well as GSEA. DEGs were enriched in many distinct signaling pathways (Figure 5A,5B). The high-risk cohort was mainly enriched in pathways involving extracellular matrix organization (e.g., collagen assembly), tissue development (including angiogenesis, skin development, and cartilage formation), chemokine-mediated regulation of immune cell migration, and immune defense responses. In contrast, the low-risk cohort was predominantly enriched in those related to energy production (mitochondrial metabolism), protein synthesis (ribosomal function), and the metabolism of lipids like steroids and retinol (Figure 5C,5D). These distinct pathways may have important implications for tumor prognosis and could provide insights for future treatment strategy formulation.

Figure 5 EA between high- and low-risk groups of the acetylation homeostatic model. (A,B) GO and KEGG EAs of DEGs across risk cohorts. (C,D) GSEA showing pathways related to acetylation homeostasis. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; EA, enrichment analysis; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; NES, normalized enrichment score.

Relation of the acetylation homeostatic model to immune cell infiltration (ICI)

ICI is critical in the progression of BLCA. Our study first compared the distribution and infiltration levels of immune cells across risk cohorts. The high-risk cohort displayed higher M0/M2 and memory B cell infiltration. The low-risk cohort exhibited increased infiltration of dendritic, CD4⁺ memory T, and regulatory T cells (Figure 6A,6B). These findings highlight the heterogeneity of ICI within the TME and reveal immune cell subpopulations correlating with tumor progression and prognosis. The relation of acetylation homeostatic model genes to ICI levels was further examined (Figure 6C), providing insights into how specific genes may influence ICI and function. In addition, correlations among immune cell types were analyzed to reveal potential interactions within the immune microenvironment (Figure 6D). Finally, ESTIMATE, immune, and stromal scores were computed. All were markedly higher in the high-risk cohort than in the low-risk cohort (Figure 6E), suggesting enhanced infiltration of certain immune cell subsets, which may be linked to poorer prognosis.

Figure 6 Immune infiltration characteristics of the acetylation homeostatic model. (A,B) Distribution and infiltration differences of immune cell subpopulations between high- and low-risk groups. (C) Correlation analysis between prognostic genes and ICI levels. (D) Correlation analysis among immune cell infiltrates. (E) Boxplots of ESTIMATE scores, immune scores, and stromal scores between high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. ICI, immune cell infiltration.

DSA

To further enhance the acetylation homeostatic model’s clinical utility for BLCA treatment, the IC₅₀ of commonly used BLCA-related drugs was calculated to compare drug sensitivity across risk cohorts (Figure 7). The high-risk cohort was more sensitive to MK2206, omipalisib, ZSTK474, lapatinib, PHA-665752, crizotinib, Z-LLNle-CHO, chloromethyl ketone, PAC-1, and 5-fluorouracil. The low-risk cohort showed greater sensitivity to vandetanib, all-trans retinoic acid, and pozaponib. Therefore, risk scores derived from the prognostic model may be used to guide personalized drug selection and optimize treatment.

Figure 7 Drug sensitivity analysis based on the acetylation homeostatic model. IC₅₀ differences between high- and low-risk groups for common chemotherapeutic and targeted agents. **, P<0.01; ***, P<0.001. IC₅₀, half-maximal inhibitory concentration.

Single-cell analysis of the acetylation homeostatic model

GSE130001 scRNA-seq data were analyzed to elucidate the roles of different cell sorts in specific BPs. A quality control step was initially applied to filter out cells containing fewer than 200 or over 4,000 detected genes, along with cells where mitochondrial genes represented over 25% of total gene expression (Figure 8A). Subsequently, clustering analysis was conducted to identify distinct cell populations, resulting in 10 clusters (Figure 8B). Cell type annotation classified these clusters into epithelial, endothelial, fibroblasts, and tissue stem cells (Figure 8C). Epithelial cells comprised the predominant population, and their characterization established a basis for biological investigations. Thereafter, the single-cell expression profiles of key genes from the prognostic model were analyzed, enabling the characterization of cell-type-specific gene expression features across different cell populations (Figure 8D).

Figure 8 Single-cell RNA-seq analysis of the acetylation homeostatic model genes. (A) Quality control of scRNA-seq data from the GSE130001 cohort. (B,C) Cell clustering and cell-type annotation of the GSE130001 dataset. (D) Cell-type-specific expression patterns of prognostic genes in the GSE130001 cohort. (E) Quality control of scRNA-seq data from the GSE135337 cohort. (F,G) Cell clustering and cell-type annotation of the GSE135337 dataset. (H) Cell-type-specific expression patterns of prognostic genes in the GSE135337 cohort.

The GSE135337 scRNA-seq dataset was analyzed to characterize the functional heterogeneity of distinct cell populations within the BLCA TME. Quality control was first performed by excluding cells with fewer than 200 or more than 4,000 detected genes, as well as cells in which mitochondrial gene expression exceeded 25% of total transcripts (Figure 8E). After quality filtering, unsupervised clustering identified 21 distinct cell clusters (Figure 8F). Cell annotation classified these clusters into epithelial cells, fibroblasts, dendritic cells, and B lymphocytes (Figure 8G). Epithelial cells and fibroblasts represented the predominant cell populations, providing the basis for subsequent biological analyses. The single-cell expression patterns of prognostic signature genes were further examined across different cell types to characterize their cell-type-specific expression profiles (Figure 8H).

CES1 was validated as the key gene regulating acetylation homeostasis in BLCA

KM survival analysis showed that patients with high CES1 expression had significantly poorer survival than those with low CES1 expression, indicating that elevated CES1 levels may be linked to worse clinical outcomes (Figure 9A). Next, CES1 expression between normal and tumor tissues was compared. CES1 was typically upregulated in tumor tissues in comparison to normal ones (Figure 9B,9C), suggesting CES1’s influence on tumor initiation and progression.

Figure 9 CES1 was validated as the key gene regulating acetylation homeostasis in BLCA. (A) KM survival curves showing the relation of CES1 expression to OS. (B) Boxplots showing differential CES1 expression between normal and tumor tissues. (C) IHC images showing CES1 expression in normal and tumor tissues (×200). Representative immunohistochemistry image was obtained from the Human Protein Atlas (https://www.proteinatlas.org). (D,E) qRT-PCR and WB analyses demonstrating downregulation of CES1 mRNA and protein levels after CES1 knockdown in tumor cells. (F) CCK-8 assay measuring OD values to assess cell proliferation. (G) Wound healing assay comparing wound closure rates to evaluate cell migration. (H) EdU assay comparing the proportion of EdU-positive cells to assess cell proliferation (×200). (I) Colony formation assay comparing colony numbers to evaluate cell proliferation and migration. Cells were fixed and stained with 0.5% crystal_violet staining (×10). (J) Measurement of acetyl-CoA levels in control and experimental groups to assess the effect of CES1 on the protein acetylation process. **, P<0.01; ***, P<0.001; ****, P<0.0001. BLCA, bladder cancer; CCK-8, Cell Counting Kit-8; CoA, coenzyme A; IHC, immunohistochemistry; OD, optical density; OS, overall survival; qRT-PCR, quantitative real-time polymerase chain reaction; WB, Western blot.

To further unveil the functional role of CES1 in BLCA, CES1 expression was knocked down in two BLCA cell lines (T24 and UMUC3) using siRNA. Knockdown efficiency was confirmed by qRT-PCR and WB analyses, which demonstrated significant downregulation of CES1 mRNA and protein levels in CES1-silenced cells (Figure 9D,9E; Figure S2). CES1 knockdown markedly reduced OD values and cell viability in both cell lines (Figure 9F). Wound healing assays demonstrated that CES1 silencing markedly slowed wound closure and impaired cell migratory capacity (Figure 9G). Upon CES1 knockdown, EdU incorporation assays showed a marked drop in EdU-positive cells (Figure 9H). Colony formation assays further showed that CES1 silencing reduced colony formation rates and significantly decreased the number of colonies (Figure 9I). Therefore, CES1 promotes BLCA cell growth, invasion, and migration.

In addition, intracellular CoA was measured in control and CES1-knockdown cells. CoA levels were markedly reduced in both cell lines following CES1 knockdown (Figure 9J), suggesting that CES1 may influence PA processes by regulating intracellular CoA availability. Furthermore, acetylation sites were predicted among the predicted post-translational modification sites of CES1 (Figure S3).

High CES1 expression predicted poor prognosis and was upregulated in BLCA tissues. CES1 knockdown significantly suppressed cell proliferation, migration, and invasion while reducing intracellular CoA levels. However, the current experimental evidence is insufficient to determine whether CES1 directly regulates PA homeostasis or indirectly influences acetylation status through global metabolic reprogramming. Consequently, the causal relationship among CES1, intracellular CoA metabolism, and PA remains to be established. Therefore, CES1 is a candidate critical gene associated with the acetylation homeostatic signature rather than a confirmed direct regulatory factor.

Discussion

BLCA represents a highly prevalent and lethal malignancy of the urinary tract worldwide, with its incidence and mortality continuously rising, thus imposing a substantial burden on global public health. Despite continuous advances in surgical techniques, chemotherapy, targeted therapy, and immune checkpoint inhibition (21), clinical outcomes for patients with advanced or metastatic BLCA remain unsatisfactory, highlighting the urgent need to identify robust prognostic biomarkers and effective therapeutic targets. PA, a pivotal post-translational modification, governs diverse cellular BPs and is critically implicated in tumor initiation and progression by modulating oncogenic signaling pathways and metabolic homeostasis (22,23). Dysregulation of PA has been linked to malignant transformation, metabolic reprogramming, and TME remodeling in multiple cancer types (24,25). However, systematic investigations into the clinical significance, prognostic value, and biological roles of acetylation-related genes in BLCA are still limited.

In this study, we constructed and validated a novel acetylation homeostatic model for BLCA using integrated transcriptome analysis and machine learning. This model effectively stratified patients into high- and low-risk groups and served as an independent prognostic factor. High-risk patients exhibited distinct molecular pathways, immune microenvironment, and drug sensitivity patterns, which were closely linked to acetylation homeostasis.

EAs of DEGs across risk groups indicated that these DEGs were involved in multiple distinct signaling pathways, including extracellular matrix organization, lipid metabolism (LM), and mitochondrial metabolism. The biological pathways identified in this study were inferred solely from bioinformatic enrichment analyses. No functional experiments were performed to validate their biological roles in BLCA, so the causal links between these pathways and tumor prognosis remain unclear and will be explored in our follow-up research. Immunotherapy is a novel treatment modality and has improved survival outcomes for patients with BLCA worldwide (26). Accordingly, immune infiltration analyses were carried out in both risk cohorts to further investigate how ICI within the TME influences tumor development and progression. The results revealed marked heterogeneity in ICI, with the high-risk cohort displaying higher infiltration of certain immune cell subpopulations. Therefore, the extent and composition of ICI are possibly related to tumor progression and prognosis. Antibody-drug conjugates have brought transformative advances to the treatment of BLCA (27) and are promising in treating metastatic BLCA (28). DSA showed that both risk cohorts were sensitive to many targeted agents, indicating that risk-stratification-guided therapeutic decision-making may further prolong survival in patients across different risk categories.

Finally, CES1 was validated as the key gene in this acetylation homeostatic model. CES1 primarily modulates cellular LM (29). LM is essential to virtually all BPs, and LM dysregulation is among the most prominent metabolic alterations in cancer (30). Members of the CES family have been extensively studied in oncology; for example, CES1 is related to poor prognosis in head and neck squamous cell carcinoma (31), CES2 is closely linked to lung cancer (32), and CES3 is a biomarker and possible treatment target in colon adenocarcinoma (33). In this study, CES1 modulated intracellular CoA levels and further regulated PA homeostasis, thereby promoting BLCA cell proliferation, invasion, and migration. These findings provide a mechanistic link between CES1, CoA metabolism, and acetylation homeostasis in BLCA.

Single-cell transcriptomic analyses of the GSE130001 and GSE135337 datasets demonstrated that CES1 was predominantly expressed in fibroblasts but showed relatively low expression in malignant epithelial cells. This cell-type-specific expression pattern suggests that CES1 may participate in regulating stromal composition and stromal-tumor communication within the BLCA TME, providing mechanistic support for our prognostic model. Future studies will isolate primary bladder fibroblasts and tumor epithelial cells to independently manipulate CES1 expression, thereby elucidating how fibroblast-derived CES1 influences extracellular matrix remodeling, ICI, and the malignant progression of adjacent tumor cells.

Nevertheless, there are limitations. First, the molecular mechanisms by which CES1 regulates PA remain incompletely understood. Although intracellular CoA concentrations were measured following CES1 knockdown, key regulators of PA, including histone acetyltransferases (HATs) and histone deacetylases (HDACs), were not evaluated, nor were lysine acetylation levels of downstream target proteins examined. Therefore, it remains unclear whether CES1 directly regulates acetylation homeostasis or indirectly affects PA through global metabolic reprogramming. The exact causal regulatory cascade shall be clarified through multi-omics and in-depth molecular experiments. Future studies shall unveil the molecular mechanisms of CES1-mediated PA leveraging multi-omics approaches, like proteomics and metabolomics, and isotope tracing. Second, the focus is merely on CES1. Although CES1 has been relatively underreported among acetylation-related DEGs and its functional role remains insufficiently characterized, whether it can comprehensively represent acetylation-related mechanisms in BLCA warrants further investigation. Third, the translational relevance of our findings is primarily constrained by their derivation from limited cell lines without in vivo validation. Fifth, only CES1 was subjected to in vitro functional validation. The remaining four prognostic signature genes, together with the hub genes identified from the PPI network (CDK1, CCNA2, CHEK1, and TOP2A), have not yet been experimentally validated. This limitation reduces the overall robustness of the acetylation homeostatic signature. Separate gene knockdown and overexpression experiments will be conducted for these molecules in subsequent research to complete systematic verification of the full prognostic signature. Furthermore, model robustness may have been affected by the limited sample size, potentially influencing the model’s stability and generalizability. Reliability and precision shall be increased by expanding sample diversity and incorporating in vivo validation. Finally, the absence of validation using clinical samples, particularly multicenter cohort or clinical tissue-based studies, limits the clinical applicability. Further validation of this model with clinical specimens and in multicenter studies is necessary to facilitate its clinical applicability and translation. Nonetheless, our acetylation homeostatic model and key gene CES1 provide new insight into BLCA prognostic evaluation and targeted therapy.

Conclusions

In summary, we successfully developed and validated a novel acetylation homeostatic model for BLCA using integrated transcriptome analysis and machine learning, which exhibits reliable performance for survival prediction and risk stratification. This model is independently associated with patient prognosis and may help guide individualized therapeutic strategies and drug selection. Furthermore, cellular experiments confirmed that CES1 is closely associated with CoA metabolism and the malignant progression of BLCA. Its potential association with PA homeostasis needs further mechanistic exploration. Collectively, our findings highlight the importance of acetylation homeostasis in BLCA and identify CES1 as a promising biomarker and therapeutic target.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the MDAR and TRIPOD reporting checklists. Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0435/rc

Data Sharing Statement: Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0435/dss

Peer Review File: Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0435/prf

Funding: This study was supported by the National Natural Science Foundation of China (grant No. 82503143).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0435/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. No animal experiments or human tissue specimens were involved. Therefore, additional ethics committee approval was not required.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA Cancer J Clin 2024;74:12-49. [Crossref] [PubMed]
2. Flaig TW, Spiess PE, Abern M, et al. NCCN Guidelines® Insights: Bladder Cancer, Version 2.2022. J Natl Compr Canc Netw 2022;20:866-78. [Crossref] [PubMed]
3. Stein JP, Skinner DG. Radical cystectomy for invasive bladder cancer: long-term results of a standard procedure. World J Urol 2006;24:296-304. [Crossref] [PubMed]
4. Shen J, Liu X, Li C, et al. Immune-related genes can accurately predict survival in bladder cancer: a retrospective study via two independent immunotherapy cohorts. Transl Androl Urol 2025;14:1661-78. [Crossref] [PubMed]
5. Sheng X, Wang L, He Z, et al. Efficacy and Safety of Disitamab Vedotin in Patients With Human Epidermal Growth Factor Receptor 2-Positive Locally Advanced or Metastatic Urothelial Carcinoma: A Combined Analysis of Two Phase II Clinical Trials. J Clin Oncol 2024;42:1391-402. [Crossref] [PubMed]
6. Li R, Linscott J, Catto JWF, et al. FGFR Inhibition in Urothelial Carcinoma. Eur Urol 2025;87:110-22. [Crossref] [PubMed]
7. Fléchon A, Morales-Barrera R, Powles T, et al. Association of Tumor Mutational Burden and PD-L1 with the Efficacy of Pembrolizumab with or without Chemotherapy versus Chemotherapy in Advanced Urothelial Carcinoma. Clin Cancer Res 2024;30:5353-64. [Crossref] [PubMed]
8. Pan S, Chen R. Pathological implication of protein post-translational modifications in cancer. Mol Aspects Med 2022;86:101097. [Crossref] [PubMed]
9. Baumann K. Post-translational modifications: Crotonylation versus acetylation. Nat Rev Mol Cell Biol 2015;16:265. [Crossref] [PubMed]
10. Gong F, Chiu LY, Miller KM. Acetylation Reader Proteins: Linking Acetylation Signaling to Genome Maintenance and Cancer. PLoS Genet 2016;12:e1006272. [Crossref] [PubMed]
11. You L, Nie J, Sun WJ, et al. Lysine acetylation: enzymes, bromodomains and links to different diseases. Essays Biochem 2012;52:1-12. [Crossref] [PubMed]
12. Meng W, Li N, Chen B, et al. Reliable prognostic definition and immunotherapy response prediction in bladder cancer, based on a novel aging-associated 5-gene signature model. Transl Androl Urol 2024;13:193-208. [Crossref] [PubMed]
13. Geffen Y, Anand S, Akiyama Y, et al. Pan-cancer analysis of post-translational modifications reveals shared patterns of protein regulation. Cell 2023;186:3945-3967.e26. [Crossref] [PubMed]
14. Li Y, Seto E. HDACs and HDAC Inhibitors in Cancer Development and Therapy. Cold Spring Harb Perspect Med 2016;6:a026831. [Crossref] [PubMed]
15. Houston R, Sekine S, Calderon MJ, et al. Acetylation-mediated remodeling of the nucleolus regulates cellular acetyl-CoA responses. PLoS Biol 2020;18:e3000981. [Crossref] [PubMed]
16. Menzies KJ, Zhang H, Katsyuba E, et al. Protein acetylation in metabolism - metabolites and cofactors. Nat Rev Endocrinol 2016;12:43-60. [Crossref] [PubMed]
17. Xie Q, Hu B, Li H. Acetylation- and ubiquitination-regulated SFMBT2 acts as a tumor suppressor in clear cell renal cell carcinoma. Biol Direct 2024;19:37. [Crossref] [PubMed]
18. Wang DN, Ni JJ, Li JH, et al. Bacterial infection promotes tumorigenesis of colorectal cancer via regulating CDC42 acetylation. PLoS Pathog 2023;19:e1011189. [Crossref] [PubMed]
19. Rios Garcia M, Steinbauer B, Srivastava K, et al. Acetyl-CoA Carboxylase 1-Dependent Protein Acetylation Controls Breast Cancer Metastasis and Recurrence. Cell Metab 2017;26:842-855.e5. [Crossref] [PubMed]
20. Saisomboon S, Kariya R, Mahalapbutr P, et al. Augmented Global Protein Acetylation Diminishes Cell Growth and Migration of Cholangiocarcinoma Cells. Int J Mol Sci 2024;25:10170. [Crossref] [PubMed]
21. Liu J, Tang X, Cui Y, et al. Development and validation of a gemcitabine sensitivity-related long noncoding RNA signature for predicting the prognosis and subtypes of bladder cancer. Transl Androl Urol 2026;15:60. [Crossref] [PubMed]
22. Pietrocola F, Galluzzi L, Bravo-San Pedro JM, et al. Acetyl coenzyme A: a central metabolite and second messenger. Cell Metab 2015;21:805-21. [Crossref] [PubMed]
23. Guertin DA, Wellen KE. Acetyl-CoA metabolism in cancer. Nat Rev Cancer 2023;23:156-72. [Crossref] [PubMed]
24. Lisek M, Bochenska N, Tomczak J, et al. Epigenetic Regulation in Ischemic Neuroprotection: The Dual Role of HDACs and HATs in Neuroinflammation and Recovery. Antioxidants (Basel) 2025;14:1015. [Crossref] [PubMed]
25. Schug ZT, Peck B, Jones DT, et al. Acetyl-CoA synthetase 2 promotes acetate utilization and maintains cancer cell growth under metabolic stress. Cancer Cell 2015;27:57-71. [Crossref] [PubMed]
26. Bellmunt J, de Wit R, Fradet Y, et al. Putative Biomarkers of Clinical Benefit With Pembrolizumab in Advanced Urothelial Cancer: Results from the KEYNOTE-045 and KEYNOTE-052 Landmark Trials. Clin Cancer Res 2022;28:2050-60. [Crossref] [PubMed]
27. Stecca C, Abdeljalil O, Sridhar SS. Metastatic Urothelial Cancer: a rapidly changing treatment landscape. Ther Adv Med Oncol 2021;13:17588359211047352. [Crossref] [PubMed]
28. Lock D, Soni YS, Hwang L, et al. Safety of Combining Radiation Therapy and Antibody Drug Conjugates in Advanced Urothelial and Other Cancers. Adv Radiat Oncol 2026;11:101926. [Crossref] [PubMed]
29. Selvaraj R, Lian J, Watts R, et al. CES1 Increases Hepatic Triacylglycerol Synthesis Through Activation of PPARγ, LXR and SREBP1c. Cells 2025;14:1548. [Crossref] [PubMed]
30. Bian X, Liu R, Meng Y, et al. Lipid metabolism and cancer. J Exp Med 2021;218:e20201606. [Crossref] [PubMed]
31. Jiang C, Liu C, Yao XI, et al. CES1 is associated with cisplatin resistance and poor prognosis of head and neck squamous cell carcinoma. Oncol Res 2024;32:1935-48. [Crossref] [PubMed]
32. Xie Q, Cao Z, You W, et al. Ganodermanontriol Suppresses the Progression of Lung Adenocarcinoma by Activating CES2 to Enhance the Metabolism of Mycophenolate Mofetil. J Microbiol Biotechnol 2024;34:249-61. [Crossref] [PubMed]
33. He L, Zhao C, Xu J, et al. A potential novel biomarker: comprehensive analysis of prognostic value and immune implication of CES3 in colonic adenocarcinoma. J Cancer Res Clin Oncol 2023;149:13239-55. [Crossref] [PubMed]