The scRNA-seq data for bladder cancer, sourced from Gene Expression Omnibus (https://ww.ncbinlm.nih.gov/geo/) (GSE135337), was utilized in this study. Data pertaining to bulk RNA-seq was acquired from the Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/), which included clinical details (age, gender, ethnicity) and somatic mutation information for bladder cancer patients. By leveraging publicly accessible data, the quirement for ethical approval was bypassed.

The gene expression data was processed in R (v4.2.0) with the Seurat package (v4.3.0). Low-quality cells were filtered out based on stringent criteria, including nFeature range (300–5000), nCount range (500–50,000), and limits on mitochondrial ≤ (1) and erythrocyte gene expression ≤ (5 contributions to the total gene count. A total of 40,167 cells were obtained after quality control.

Seurat’s NormalizeData prepared the gene expression data. The FindVariableFeatures function pinpointed the 2,000 most variable genes (6–10). ScaleData functions (v3.1.4) prepared the gene expression data (11), guiding subsequent principal component analysis. The harmony R package (v.0.1.1) (12, 13) mitigated batch effects between sample.CellCycleScoring is used to calculate cell cycle phases (14). Clusters were formed based on the top 30 principal components (PCs). For UMAP visualization, the same 30 PCs were selectively used to depict gene expression patterns (15–18). A dimensionality reduction of 30 and a resolution of 1.2 were used.

Cell clusters were initially identified with Seurat’s FindClusters and FindNeighbors functions (18–21). Clusters were annotated based on marker genes expression averages.

Slingshot (v2.6.0) inferred cell lineages and pseudotimes, leveraging clustering-based MSTs to define lineage architecture. Synchronized master curves and branching curve fittings were applied, with getCurves used to derive smooth trajectory curves.

The AUCell method (22) is utilized for the identification of cells exhibiting active gene single-cell RNA-seq data underwent analysis to assess gene expression profiles. It takes gene sets as input and provides an assessment of the ‘activity’ of each gene set within every individual cell. In this particular study, it was employed to evaluate the stemness characteristics associated with different tumor cell subpopulations.

Differentially expressed genes (DEGs) were identified using FindAllMarkers, applying a Wilcoxon test with thresholds of min.pct and min.diff.pct at 0.25, and a logfc threshold of 0.25 (23). ClusterProfiler (v4.6.0) facilitated DEGs enrichment and analysis using GO (24), KEGG, and GSEA (25–27). Significant GO terms were identified with an adjusted P-value < 0.05.

CellChat (v1.6.1) enables quantitative inference and analysis of scRNA-seq data to decipher cellular interactions. CellChat’s netVisual diffInteraction function assesses alterations in cell-cell communication strength, while IdentifyCommunicationPatterns determines the count of distinct communication patterns. Scatter plots, heatmaps, and various visualization techniques are used to visually analyze the signals entering and exiting each cell. The CellChatDB database (http://www.cellchat.org/) is then used to identify signaling pathways and receptor pairs associated with specific types of TCs in cancer. A P-value threshold of 0.05 was set to identify statistically significant cell-cell interactions across various cell types.

Using pySCENIC (v0.10.0) in Python 3.7, we inferred single-cell regulatory networks and clustered tumor cell subpopulations. GRNBoost identified TF target genes, which were then refined by DNA-motif analysis. AUCell scored regulon activities, revealing the top 5 transcription factors (TFs) with most significant changes in expression, based on human gene motif rankings from https://resources.aertslab.org/cistarget/.

The purpose of this research was to evaluate the prognostic potential of specific prognostic genes associated with diverse bladder cancer subpopulations in predicting patient survival outcomes. Through a rigorous process involving univariate and multivariate Cox proportional hazards analysis (28, 29), combined with Least Absolute Shrinkage and Selection Operator (LASSO) regression, we pinpointed the most influential prognostic genes as key predictors for developing a robust prognostic model. Subsequently, we formulated a risk scoring system, where the risk score is calculated as the sum of the products between individual gene expression levels and their respective coefficients.

Employing an optimized cutoff threshold derived from the “surv_cutpoint” function, patients were stratified into distinct low risk and high STRS groups, allowing for a comparative analysis of prognostic differences across patient subpopulations. To visualize and statistically validate the predictive capabilities of our risk score model, we leveraged the “Survival” package (v3.3.1) in R for survival analysis and utilized the “ggsurvplot” function to generate survival curves (30–32). Moreover, we ensured the reliability of our model by assessing its accuracy and calibration through the generation of Receiver Operating Characteristic (ROC) curves (33–37) using the “timeROC” package (v0.4.0), providing a comprehensive evaluation of the model’s predictive performance.

We utilized the CIBERSORT R package (version 0.1.0) to compute immune-related scores for immune cells, providing a comprehensive evaluation of the patients’ immune milieu. Furthermore, we scrutinized the levels of immune cell infiltration and differential gene expression associated with immune checkpoints, and conducted investigations into the correlation between risk scores, immune cells, and model genes (38). Concurrently, we employed the Tumorlmmune Dysfunction program to appraise the response to tumor immunotherapy.

To distinguish between cancerous and non-cancerous cellular populations, we initially estimated the baseline copy number variation (CNV) across various genomic regions by analyzing disruptions in chromosome gene expression patterns. This was accomplished utilizing the inferCNV R package (accessible at https://github.com/broadinstitute/inferCNV/wiki), a tool specifically designed for CNV inference. Using endothelial cells (ECs) as a benchmark, we leveraged the inferCNV algorithm to characterize the CNV landscape within distinct cellular subpopulations. Subsequently, those EPCs subpopulations that exhibited marked alterations in their CNV profiles, indicative of significant genomic instability, were identified and classified as malignant cells, thereby facilitating the differentiation between cancerous and non-cancerous cell types.

The UM-UC-1 tumor-derived cell line was propagated in MEM medium under standardized conditions (37°C, 5% CO2 atmosphere, and 95% humidity) supplemented with 10% fetal bovine serum (FBS) and 1% antibiotics. Similarly, the VM-CUB1 cell line was maintained in DMEM medium under identical conditions, also enriched with 10% FBS and 1% antibiotics to ensure optimal growth and health.

RNA constructs sourced from GenePharma (Suzhou, China) facilitated the downregulation of KDELR2 expression. The cells were seeded onto a 6-well plate at a moderate density of 50%, subsequently subjected to transfection procedures involving KDELR2-specific knockdown constructs (si-KDELR2-1 and si-KDELR2-2), as well as a negative control construct (si-NC) for comparison. Lipofectamine 3000RNAiMAX (Invitrogen, USA) was used for transfection under manufacturer directions. Every si-RNA (RIbbio, China) was transfected into cells. siRNA sequences: si-1: GUAGUCCAGACCAUCCUAU; si-2: UCGUGCUUUGUAUCUUGUC. qRT-PCR Primers: F: TGGATCTGGCGCTTCTACTT; R: GCTGGCAAACTGAGCTTCTT.

After the transfected cells attained a 70% confluency level, they were lysed in RIPA buffer. The resulting lysates were clarified through centrifugation at 12,000 rpm for 15 minutes, preparing them for SDS-PAGE separation. The separated proteins were then transferred onto PVDF membranes, which were subsequently blocked with 5% BSA for 1.5 hours at ambient temperature. Following an overnight incubation with an Anti-KDELR2 antibody at 4°C, the membranes were further incubated with a secondary antibody for an hour. Ultimately, the presence of KDELR2 protein bands was detected using an ECL Western Blot substrate for visualization.

The RNA extraction process involved the utilization of Trizol reagent, Trizol reagent was used to lyse cells and release RNA, with chloroform and isopropanol employed to precipitate the RNA while suppressing RNase activity. Throughout the extraction process, it was essential to confirm that the workbench, tools, and water utilized were RNase-free to avoid RNA degradation. Followed by a reverse transcription step facilitated by the PrimeScript™ Kit. Subsequently, the quantitative Real-time Polymerase Chain Reaction was performed using SYBR Green premix as the fluorescent dye for amplification detection.

To assess the viability of UM-UC-1 and VM-CUB1 cells post-transfection, the Cell Counting Kit-8(CCK-8) assay was employed. Cells were seeded in 96-well plates at a density of 5×10³ cells per well and allowed to incubate for 24 hours. Subsequently, 10μL of CCK-8 labeling reagent (A311-01, Vazyme) was added to each well, and the plates were incubated in the dark at 37°C for two hours. From day one to day four, cell viability was quantitatively determined by measuring the absorbance at 450nm using an enzymatic marker (A33978, Thermo). The average optical density values were computed and plotted on a line graph to visually represent the cellular viability trends over time.

Prior to the experiment, the cells underwent a 24-hour serum starvation period in medium devoid of serum. Afterward, the cell suspension was mixed with Matrigel (BD Biosciences, USA) and seeded into the upper chamber of Costar plates, while the lower chamber was filled with serum-rich medium to create a chemoattractant gradient. The cells were then incubated for 48 hours in a cell culture incubator to allow for migration and invasion. Following incubation, the cells were fixed with 4% paraformaldehyde and stained with crystal violet to visually assess their invasive potential.

The stably transfected cells were plated in 6-well dishes placed in the incubator and cultured at 37°C with 5% CO₂ until the cells reached confluence. Using a sterile 200μL pipette tip, uniform scratches were generated across the cell monolayer in each well. Subsequently, the wells were gently washed with PBS to remove any dislodged cells or debris. The scratched areas were then subjected to incubation in serum-free medium to monitor cell migration. Images of the scratch wounds were documented at 0 hours and again after 48 hours of incubation, with the widths of the scratches measured utilizing Image-J software for quantitative analysis.

The UM-UC-1 and VM-CUB1 cell lines, post-transfection, were seeded at a concentration of 5×10³ cells per well in 6-well plates. Following a 24-hour incubation period at ambient temperature, the EdU working solution was introduced into the culture medium and allowed to incubate for 2 hours. Subsequently, the cells underwent a double wash with PBS and were fixed using a 4% paraformaldehyde solution for 15 minutes to stabilize them. After fixation, the cells were permeabilized and quenched with a mixture of 2 mg/ml glycine and 0.5% Triton X-100 for 15 minutes. Finally, the cells were stained with a combination of 1X Apollo solution (1 ml) and Hoechst staining reaction solution (1 ml), followed by a 30-minute incubation period. Fluorescence microscopy was then employed to assess cell proliferation by capturing images of the stained cells.

R and Python software packages are utilized for analyzing data from databases, while GraphPad Prism, specifically version 8.0.1, serves as the tool of choice for experimental data analysis. Throughout the analyses, two-tailed p-values are employed, and statistical significance is determined based on values falling below the threshold of 0.05.

Based on our research objectives and requirements, we conducted a comprehensive review of the scRNA-seq data related to bladder cancer in public database. To explore the cellular transcriptome characteristics of the bladder cancer TME, we performed a systematic and comprehensive analysis of the collected scRNA-seq data. The dataset was obtained from seven primary tumor tissue samples and one paracancerous sample from seven bladder cancer patients that included diverse demographic backgrounds and clinical stages of bladder cancer. This diversity enhanced the relevance of our findings to a broader patient population. After removing batch effects and primary quality control, we obtained 40,167 high-quality cells ( Figure 1A ). Dimension reduction cluster analysis was performed on the screened cells, and we found that they could be divided into 5 cell types, including B plasma cells, ECs, epithelial cells (EPCs), myeloid cells, and fibroblasts ( Figure 1B ). The bubble plot depicted the levels of expression for the most significant genes in distinct cellular subpopulations, categorized by cell types and tissue samples ( Figure 1C ). The results found that ECs highly expressed PLVAP, GNG11, SPARC, RGCC, IGFBP7; fibroblasts highly expressed CFD, LUM, DCN, GSN, MT2A; EPCs showed high expression of SPINK1, KRT19, LY6D, FXYD3 and S100A2; B plasma cells showed high expression of IGHG1, IGLC3, IGHA1, IGKC, IGLC2; myeloid cells showed high expression of HLA-DRA, CCL3, HLA-DPB1, C1QB and SPP1; bladder cancer tissue cells showed high expression of SPINK1, LY6D, KRT19, FXYD3, S100A2; paracancerous cells highly expressed genes including CFD, LUM, DCN, GSN, MT2A,RGCC and so on. In conclusion, we could find that EPCs marker genes had high consistency with the highly expressed genes in cancer tissues. At the same time, we could observe that the cell phases are different in these 5 cell types, and most of the cells were in G1 phase ( Figure 1D ). Moreover, most ECs and EPCs were mainly derived from bladder cancer tissue, while fibroblasts were mainly derived from paracancerous cells ( Figure 1E ). The above study were consistent with Ro/e analysis, that EPCs tended to be derived from bladder cancer tissues ( Figure 1H ). Through cell cycle analysis of cell subpopulations, each cell subpopulation tended to be distributed in different cycle phases ( Figures 1C, D, F ). In addition, the proportion of EPCs in bladder cancer tissue cells and paracancerous cells was much greater than that of other cells, especially in bladder cancer tissue cells, as high as 97.7% ( Figures 1E, G ).

Volcano plots results showed that the different genes among ECs, fibroblasts, EPCs, B plasma cells, and myeloid cells ( Figure 1I ). To investigate the biological processes of different cell types, by using GO-BP enrichment analysis ( Figure 1J ), we identified that EPCs were mainly enriched in ATP synthesis coupled electron transport, mitochondrial ATP synthesis coupled electron transport, cytoplasmic translation, proton motive force-driven mitochondrial ATP synthesis and mitochondrial respiratory chain complex assembly. Further GSEA ( Figure 1K ) of EPCs was performed according to GO-BP terms. The results showed that pathways such as mitochondrial electron transport, NADH to ubiquinone; oxidative phosphorylation and ATP synthesis coupled electron transport showed a positive enrichment trend, by comparison, pathways such as adaptive immune response and extracellular structure organization showed a trend of negative enrichment in this genome.

The findings of various studies have demonstrated that the presence of TCs can trigger significant alterations at a molecular, cellular, and physical level within the surrounding tissue, leading to the formation of a specialized environment known as the TME. During the initial stages of tumor development, there is an intricate interplay between cancer cells and various components of this microenvironment, facilitating tumor survival, local infiltration, and metastatic dissemination (39). Therefore, TCs are the real culprit of tumor development. We first used the infer CNV algorithm to successfully separate 19,621 bladder cancer cells from EPCs based on CNV ( Supplementary Figure S1 ). Next, we identified six subpopulations of bladder cancer cells based on the level of marker genes expression. The circle diagram showed six cell subpopulations. These were C0 FABP4+ TCs, C1 S100A8+ TCs, C2 TFF2+ TCs, C3 CRH+ TCs, C4 BIRC5+ TCs and C5 IL32+ TCs. In addition, we used UMAP plots to show the nCount-RNA, nFeature-RNA, G2M.Score and S.Score of each subpopulation ( Figure 2A ). Subsequently, the bubble plots showed the marker genes of different tissue cells and subpopulations ( Figure 2B ). According to the comparison of named gene expression of each subpopulation of tumor cell in bar charts, C0 was characterized by high expression of FABP4. The high expression of S100A8 characterizes C1, while C2 was characterized by the elevated expression of TFF2. Similarly, C3 was distinguished by its heightened expression of CRH, and C4 exhibits a notable increase in BIRC5 expression. Lastly, the high expression of IL32 served as a defining feature for C5. UMAP plots showed the named gene expression of each subpopulation, which was consistent with the results of bar charts ( Figure 2C ).

For further analysis of the subpopulations of TCs, we conducted deeper visualization study. First, we performed sample sources and cell cycle phases analysis of TCs, and combined the UMAP plots to observe the main distribution of different TCs ( Figure 2D ). The results showed that most of the C1 subpopulation of TCs in G2M stage, compared to the other subpopulations, the highest proportion, as high as 50.90%,at the same time, the tissue classification results show that the C1 subpopulation is different from other subpopulations, its mainly come from bladder cancer tissue cells ( Figures 2E, F ), and the CNV score was relatively high ( Figure 2G ). In addition, we conducted CNV score analysis of bladder cancer tissue cells and paracancerous cells, the cancer tissue cell expression level of CNV score was higher than the paracancerous cells. There were significant differences between them. UMAP plot showed the CNV score distribution characteristics. Additional examination of the stemness characteristics of specific subpopulations of cells in bladder cancer demonstrated that the C1 subpopulation had the highest Area Under the Curve (AUC) score for stemness. This meant that C1 had the characteristics of low differentiation degree and strong differentiation potential. The stemness score of bladder cancer cells exhibited a slightly higher level compared to that of paracancerous cells; however, the difference did not reach statistical significance ( Figure 2H ).

Subsequently, we conducted an analysis of gene expression related to cell stemness in the subpopulations of TCs. And we found that the stemness gene expression level was more significant in the C1 subpopulation compared with the other cell subpopulations ( Figure 3A ). The bar charts confirmed the above conclusion that the C1 subpopulation is higher in CTNNB1, MYC, HIF1A and BMI1 expression levels than other subpopulations cells ( Figure 3B ). Research had indicated that the atypical activation of CTNNB1 was linked to the development of various types of tumors, including but not limited to colorectal cancer, ovarian cancer, prostate cancer, hepatoblastoma, and hepatocellular carcinoma (40). HIF1A, hypoxia-inducing factor, promotes angiogenesis and is important for the vascular system in the embryo and for cancer tumors (41, 42). As a proto-oncogene, BMI1 was an important component of the polycomb gene family. High levels of BMI1 expression have been found to be significantly associated with the onset, progression, and prognosis of diverse malignancies. BMI1 could participate in tumorigenesis by inhibiting multiple gene loci, synergizing with other proto-oncogenes, and enhancing telomerase activity (43). The synergistic effect of these genes precisely demonstrateed the pro-tumor effect of C1 subpopulations. The UMAP plots showed the expression distribution of CTNNB1, MYC, HIF1A, BMI1 genes ( Figure 3C ).

Cell stemness was closely related to cell differentiation and development. To explore the differentiation trajectory of tumor cell subpopulations, we performed a pseudotime analysis and presented it using UMAP plots. First, we showed the lineage trajectory of bladder cancer tissue cells and paracancerous cells ( Figure 3D ), and the cell cycle lineage trajectory ( Figure 3E ). The initial cell differentiation trajectory began in the G1 phase and gradually transitioned to the S and G2M phases according to the cell cycle sequence. We then analyzed six different TCs subpopulations, and the UMAP plots showed two main cell lineage tracks ( Figure 3F ), including lineage 1:C0→C2→C5→C4→C3; lineage 2:C0→C2→C4→C1. The difference between the two trajectories mainly existed in the later stage. Combined with Figure 3D , we could find that the ended of lineage 1 ended mainly C3 subpopulation, which contained cells from both tumor tissues and cells from paracancerous cells. The end of lineage 2 ended mainly C1 subpopulation, which was a tumor cell subpopulation derived entirely from a bladder cancer tissue sample. Combined with the evolution of TCs and based on the tissue origin, CNV score and cell stemness of C1 subpopulation, the conclusion that C1 subpopulation has a high degree of malignancy and is closely related to the progression of bladder cancer was reaffirmed.

Furthermore, we conducted an analysis on the temporal expression patterns of marker genes belonging to 6 distinct subpopulations. The findings indicated that. different from other cell subpopulations. The high expression of marker genes for the C1 subpopulation was predominantly observed during the later stage ( Figure 3G ). To confirm the above findings, we also conducted a temporal analysis of the expression of cell stemness genes. The primary distinction between lineage1 and lineage2 was observed in the later stage, with a notable increase in stemness gene expression detected in lineage2 compared to lineage1 ( Figure 3H ). This was consistent with previous conclusions. Subsequently, we analyzed the subpopulations using GO-BP enrichment analysis to verify the related biological processes of the two lineages ( Figure 3I ), and the dynamic timing showed the expression changes of the different genes of TCs along the two trajectories within the pseudotiming.

To explore the biological functions of TCs in bladder cancer, we first analyzed the different expression genes of six tumor cell subpopulations. The results showed that the main upregulated genes in C1 S100A8+ TCs were RPL32, RPS8, BEST3, IGF2 and NUPR1 ( Figure 4A ). Among them, RPL32 played a crucial role in modulating cellular signaling pathways. It had found that RPL32 could participate in multiple signaling pathways during tumorigenesis and significantly contribute to the initiation and progression of neoplastic growth (44, 45). RPS8 was mainly involved in protein folding and stability, and in tumor studies, RPS8 has the potential to serve as a biomarker specific to tumors due to its tendency for elevated expression levels in tumor tissues and cells compared to normal tissues (46). BEST3 was a coding gene that encoded proteins with good histocompatibility and low expression levels. BEST3 had been observed to be present in multiple types of cancerous tissues, and its expression level appears to be linked with the tumor’s capacity for invasion. Furthermore, the correlation between the expression level of BEST3 and the efficacy of tumor therapy has prompted its investigation as a potential therapeutic target or biomarker (47). The role of IGF2 was pivotal in regulating cellular proliferation, development, motility, differentiation, and viability. Furthermore, IGF2 played a role in a variety of cancer development (48). NUPR1 had been linked to the onset and progression of cancer, and it was commonly observed to be upregulated in different cancer forms. Moreover, NUPR1 could act as a transcriptional regulator affecting the expression of genes involved in cell cycle regulation, apoptosis, and stress response. These genes had a significant correlation with the onset and progression of neoplasms (49, 50).

Furthermore, we conducted enrichment analysis using GO-BP and KEGG to identify enriched biological processes and pathways associated with the DEGs in the subpopulations of TCs. And the heatmap showed the results of the genes enrichment in the TCs subpopulations ( Figure 4B ). C0 FABP4+ TCs was mainly associated with cytoplasmic translation, translational initiation, intrinsic apoptotic signaling pathway, ATP synthesis coupled electron transport and mitochondrial ATP synthesis coupled electron transport biological processes related in GO-BP, and was associated with the ribosome, coronavirus disease-COVID-19,chemical carcinogenesis-reactive oxygen species, oxidative phosphorylation and parkinson disease pathways in KEGG.C1 S100A8+ TCs was mainly associated with cytoplasmic translation, protein folding, chaperone-mediated protein folding, regulation of apoptotic signaling pathway and the biological process of ribonucleoprotein complex biogenesis in GO-BP, and was associated with the fluid shear stress and atherosclerosis, ribosome, protein processing in endoplasmic reticulum, prion disease and parkinson disease pathways in KEGG. The association of C2 TFF2+ TCs primarily pertained to the biological processes related to hypoxia response, decreased oxygen levels response, oxygen levels response, cellular response to decreased oxygen levels, and cellular response to hypoxia in GO-BP, and was associated with the estrogen signaling pathway, HIF-1signaling pathway, legionellosis, lipid and atherosclerosis and the glycolysis pathways correlation in KEGG. The association of C3 CRH+ TCs primarily pertained to the process of ATP production through electron transport, mitochondrial ATP synthesis coupled with electron transport, aerobic electron transport chain, and oxidative phosphorylation in GO-BP, and was associated with the parkinson disease, ribosome, huntington disease, oxidative phosphorylation and the amyotrophic lateral sclerosis pathways correlation in KEGG.C4 BIRC5+ TCs was mainly associated with chromosome segregation, ATP synthesis coupled electron transport, mitochondrial ATP synthesis coupled electron transport, nuclear chromosome segregation, respiratory electron transport chain in GO-BP, and was associated with the parkinson disease, huntington disease, amyotrophic lateral sclerosis, prion disease and alzheimer disease pathways correlation in KEGG. C5 IL32+ TCs was mainly associated with leukocyte mediated cytotoxicity, leukocyte mediated immunity, cell killing, lymphocyte mediated immunity and regulation of T cell activation biological processes in GO-BP, and was associated with the natural killer cell mediated cytotoxicity, graft-versus-host disease, antigen processing and presentation, primary immunodeficiency and the Epstein-Barr virus infection pathways in KEGG.

In addition, we performed a further step GSEA of the six tumor cell subpopulations, and we could observe that C1 S100A8+ TCs were mainly enriched in intracellular zinc ion homeostasis ( Figure 4C ). GSEA results ( Figure 4D ) also showed that the C1 subpopulation was positively enriched in the regulation of extent of cell growth, regulation of axon extension, positive regulation of mononuclear cell migration, the intracellular zinc ion homeostasis, while it was negatively regulated in positive regulation of peptide hormone secretion, antigen processing and presentation of exogenous antigen, glycolytic process through glucose-6-phosphate and glycolytic process through fructose-6-phosphate.

In summary, C1 subpopulation was mainly associated with biological pathways such as protein folding and ribosomes, which promoted and accelerated tumor cell migration and invasion. Cells were very sensitive to physiological conditions, and when the physiological conditions changed, the cells would change accordingly. When the loaded proteins exceed the folding capacity of the endoplasmic reticulum stress (ERS), ERS will be caused. The core of ERS lied in the occurrence of protein misfolding, which significantly influenced the growth and viability of cancerous cells and may facilitate the progression of bladder carcinoma (51, 52). Ribosomes, as a kind of nutrient, can escort the smooth and rapid translation of viruses and provide conditions for tumor proliferation, thus affecting the progression and development of bladder cancer (53).

To gain a thorough and organized comprehension of intricate cellular reactions, our aim was to conduct an examination on intercellular connections and networks involved in ligand-receptor communication, with the intention of visualizing interactions between cells. By means of CellChat analysis, we successfully established a network for cellular communication encompassing the majority of cells, including ECs, fibroblasts, B plasma cells, myeloid cells and six tumor cell subpopulations ( Figure 5A ). Then, we determined the quantity of interactions through the ‘line’ connection linking two cell types, with thicker lines indicating a higher number of interaction paths. Additionally, we represented the strength of interactions using line weight, where thicker lines indicated stronger interactions. Among different cell subpopulations. Next, we combined the CellChat maps to explore how different cell subpopulations interact with each other through these pathways. First, we categorized cell communication patterns into three main types. In addition, the displayed heatmaps illustrated the expression of cell interaction proteins in the three modes. The two heatmaps on the left showed the outgoing signal patterns of different cell subpopulations and the active interacting proteins in the signal patterns, and the two heatmaps on the right showed the incoming signal patterns and interacting proteins of cell subpopulations ( Figure 5B ). We could observe that in the outgoing signal pattern diagrams, C1 S100A8+ TCs were mainly in pattern 3. Correspondingly, ncWNT, CDH, MPZ, MK, and OCLN showed high expression levels in pattern 3. Additionally, the incoming signal patterns of C1 S100A8+ TCs exhibited pattern 1 including but not limited to the pathways represented by such as EPHA, CDH, MPZ, IGF, and OCLN. The Sankey diagrams and heatmaps ( Figure 5C, D ) showed the communication pattern and target signal of each cell subpopulations, which was consistent with the results in Figures 5B .

Subsequently, we analyzed the number of cell interaction strength between the C1 S100A8+ TCs and other cell subsets in depth ( Figures 6A, B ). We could find that both in the incoming signal or the outgoing signal, the number and intensity of the interaction between myeloid cells and C1 S100A8+ TCs were more significant than other cells, so we could infer that C1 S100A8+ TCs cells and myeloid cells had a strong interaction. Subsequently, we performed signaling network analysis and showed that C1 S100A8+ TCs and myeloid cells mainly communicate with MIF signals pathway ( Figure 6C ). Then we analyzed the role of the two in the signaling pathway ( Figure 6D ), which showed that C1 S100A8+ TCs cells mainly played the role of signal sender, while myeloid cells mainly played the role of signal sender and influencer.

By comparison, we could infer that the MIF ligand of C1 S100A8+ TCs cells acted on the CD74-CXCR4 and CD74-CD44 receptor of myeloid cells ( Figures 6E, F ). The cell interaction circle diagrams ( Figures 6G, H ) showed the interaction relationship between C1 S100A8+ TCs cells and myeloid cells in MIF-(CD74+CXCR4) and MIF-(CD74+CD44) cell signaling pathway, verifying the above statement. In order to make the results more visual, we used the interaction hierarchy diagrams ( Figures 6I, J ) to show the relationship between C1 S100A8+ TCs cells and myeloid cells (the thickness of the line represents the interaction strength, the thicker of the line, the relationship is more significant). The findings indicated that the intercellular communication between C1 S100A8+ TCs and myeloid cells primarily took place through paracrine, leading to signal crosstalk. This could disrupt normal signal transduction, leading to alterations in the tumor immune microenvironment, thereby facilitating the proliferation, metastasis, and invasion of TCs (54).

TFs act on genes to regulate gene transcription by binding specific nucleotide sequences upstream of the gene, thereby affecting cell biological functions. Firstly, we conducted cluster analysis of bladder cancer tissue cells according to genes expression ( Figure 7A ). In addition, the UMAP plots displayed the expression of TFs in six tumor cell clusters ( Figure 7B ). Subsequently, based on the heatmap of TFs correlation ( Figure 7C ), we divided the TFs with similar functional and expression into three main modules, namely M1, M2 and M3. Next, we visualized the analysis of the divided modules, and observed the expression levels of each TCs subpopulation in different modules through the UMAP plots. The results showed that C1 S100A8+ TCs cells had significant expression in the M1 module ( Figure 7D ). To further verify the above statement, we performed transcription factor regulatory activity score analysis ( Figure 7E ), and we can observe that in the M1 module, C5 IL32+ TCs cells have the highest transcription factor regulatory activity score, followed by C1 S100A8+ TCs cell. The bar charts also verify the above statement that C1 S100A8+ TCs cells have high expression levels in the M1 module, second only to C5 IL32+ TCs cells ( Figure 7F ). Subsequently, we conducted an analysis on the top 5 TFs in the subpopulations of TCs ( Figure 7G ). The specificity score identified HES1, AHR, TBL1XR1, IRF4, and SPI1 as the top 5 TFs in C1 S100A8+ TCs. Finally, we used the UMAP plots to show the TFs distribution characteristics of C1 S100A8+ TCs cells ( Figure 7H ).

To investigate the prognostic factors in patients, we employed univariate Cox regression analysis to discover 21 genes that exhibited an association with prognosis ( Figure 8A ). We could observe HES1 HR<1, while HR values of other genes were>1. Therefore, HES1 was a protective factor that favors patient prognosis, while others were risk factors and unfavorable to patient prognosis. To mitigate the issue of multicollinearity among genes, we employed LASSO regression analysis to identify 10 genes that were associated with prognosis ( Figure 8B ). Multivariate Cox regression analysis was then performed on the above genes and used to calculate the genes risk coefficient for these genes ( Figures 8C, D ). The curve graph and scatter plot revealed variability in risk scores and survival outcomes between low STRS (S100A8+ tumor risk score) group and high STRS group, while high STRS group was associated with worse outcome ( Figure 8E ). The heatmap showcased the distinctive expression patterns of the prognostic genes in both cohorts, setting them apart from the rest of the gene pool ( Figure 8F ), HES1 was more significantly expressed in low STRS group and HES1 favors the prognosis of the patient, consistent with previous conclusions. The ROC curves and AUC values ( Figure 8G ) for 1-year, 3-years, and 5-years periods exceeded the threshold of 0.6, indicating that the prediction model was specific and valuable. The Kaplan-Meier survival curve further confirmed the low survival profile of high STRS group ( Figure 8H ). Combined with Figure 8D , the KDELR2 risk score was the highest. Hence, we conducted an analysis using the Kaplan-Meier survival curve ( Figure 8I ) and observed a notable decline in survival rates among patients exhibiting elevated KDELR2 gene expression.

To provide clarity on the distinction between the two scoring cohorts, an analysis was conducted on the genes that exhibited differential expression ( Figure 9A ). The volcano plot showed the up-regulation and downregulation trends of DEGs ( Figure 9B ). Subsequently, to understand the biological processes for the above genes, we performed various enrichment analyses. The first was GO enrichment analysis, which reveals the main biological process, cell composition and molecular function ( Figure 9C ). In GO-BP, genes were mainly enriched in biological processes like epidermis development, keratinization, keratinocyte differentiation, intermediate filament organization, skin development. In GO-CC, its predominant enrichment were observed in the cornified envelope, intermediate filament, and cytoskeletal components of intermediate filaments as well as keratin filament. In GO-MF, mainly enriched in structural constituent of skin epidermis, serine-type endopeptidase activity, serine hydrolase activity, endopeptidase activity and receptor ligand activity. In addition, we performed KEGG pathway enrichment analysis of differential genes ( Figure 9D ), and showed that the related pathways were staphylococcus aureus infection, retinol metabolism, ECM-receptor interaction, steroid hormone biosynthesis and cytokine-cytokine receptor interaction. Finally, we performed a GSEA of the enriched pathways ( Figure 9E ), the results showed that keratinization, keratinocyte differentiation, skin development and the epidermis development related pathways showed a positive enrichment trend, mitochondrial respiratory chain complex assembly, uronic acid metabolic process, the flavonoid metabolic process and ATP synthesis coupled electron transport related pathways showed a negative enrichment trend.

We performed immune infiltration analysis on high STRS group and low STRS group. The box plot revealed higher expression levels of macrophages M0 and M1, and neutrophils in the high STRS group ( Figure 10A ). Those cell types belong to myeloid cells. It was noteworthy that macrophages exhibited a higher abundance of immune infiltration and exerted a more pronounced effect on the high STRS group. Subsequently, we further investigated the correlation between risk score and immune infiltration ( Figure 10B ). In the correlation lollipop chart, macrophages also showed a close positive correlation with risk score. In addition, the immune-score, stromal-score, and ESTIMATE-score of the high STRS group were significantly higher than those of the low STRS group, while the tumor-purity was significantly lower than that of the low STRS group ( Figure 10C ). Scatter plots demonstrated positive associations between KDELR2 expression levels with immune-score, stromal-score, and ESTIMATE-score as well as macrophages M0 ( Figure 10D ). Elevated KDELR2 expression was linked to increased infiltration of macrophages and a more prominent immunosuppressive microenvironment. Therefore, target mining of macrophages might be a potential strategy for bladder cancer immunotherapy. The heatmap provided visual comparisons that supported our previous findings ( Figure 10E ).

Furthermore, the bubble plot ( Figure 10F ) showed the correlation between immune checkpoints and prognostic genes, OS, Risk. Among them, CD274, CD276, CD80, CD86, HAVCR2, LAIR1, NRP1, PDCD1LG2, TNFRSF8, TNFRSF9, TNFRSF4, and VTCN1 were significantly positively correlated with KDELR2. This suggested that patients with high KDELR2 expression were potential candidates for these immune checkpoint inhibitors (ICIs). In order to optimize ICIs utilization and improve therapeutic efficacy further analysis was conducted on gene expression levels of different immune checkpoints between high STRS group and low STRS group using a box plot ( Figure 10G ). The results indicated significantly higher expression levels for most genes in the high STRS group compared to those in the low STRS group. These findings suggested that application of these ICIs to patients classified as high STRS group might result in improved treatment outcomes.

In an effort to delve deeper into the impact of KDELR2 on the patient outcome in cases of bladder cancer, we conducted in vitro experiments. Firstly, we divided tumor cell lines of different sources into three groups: si-NC, si-KDELR2-1 and si-KDELR2-2. The findings indicated notable variances in the expression levels of protein and mRNA for UM-UC-1 and VM-CUB1. Specifically, the si-KDELR2 group exhibited lower expression levels compared to the si-NC group ( Figure 11A ). In addition, we investigated the effect of KDELR2 on tumor cell activity. First, tumor cell activity was tested by CCK-8, and fluorescence staining showed that the tumor cell activity in si-NC was higher than that in si-KDELR2-1 and si-KDELR2-2. Subsequently, we used line plots to quantify indicators. The OD values of si-KDELR2-1 and si-KDELR2-2 were found to be significantly reduced compared to the si-NC group, aligning with the outcomes obtained from CCK-8 staining ( Figure 11B ). Next, we studied the proliferation ability of TCs. The colony formation assay ( Figure 11C ) similarly corroborated the aforementioned statements. Through EDU experiments, the microbiota density was lower after KDELR2 was knocked out. Subsequent data analysis showed that si-KDELR2 inhibited cell proliferation and cloning ( Figures 11D, E ). In order to investigate the invasion and migration ability of the cells, we conducted wound healing assay and transwell assay ( Figures 11F, G ). After 48 hours of assay, si-KDELR2 TCs’ distance was greater than that of si-NC, while the density was lower than that of si-NC, which proved that si-KDELR2 TCs had low invasion and migration ability. The bar plots were consistent with the experimental conclusion that UM-UC-1 and VM-CUB1 cells inhibited the invasion, migration and wound healing of TCs after the knockout of si-KDELR2 ( Figure 11H ). In summary, KDELR2 had obvious tumor promoting effect, which played a key role in the development of TCs. In conclusion, our research provided a more comprehensive understanding of the KDELR2 gene. Through a series of studies, we discovered that KDELR2 promoted tumor cell development. It could be inferred that KDELR2 played a key role in driving tumor cell proliferation, migration, and other related processes. Accordingly, we could achieve the purpose of improving prognostic survival by inhibiting KDELR2, making it possible focus for therapeutic intervention in cases of bladder cancer.