Construction and validation of prognosis and treatment outcome models based on plasma membrane tension characteristics in bladder cancer

Data acquisition

After consulting relevant literature, we selected 41 MTRGs for analysis (Table S1) (Diz-Muñoz, Fletcher & Weiner, 2013; Simunovic et al., 2019; Issa & Noureddine, 2017; Ha & Chi, 2012). RNA sequencing data and clinical annotations for BLCA samples were obtained from The Cancer Genome Atlas (TCGA) (TCGA-BLCA, https://www.cancerimagingarchive.net/collection/tcga-blca/) and Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds). The TCGA database includes 19 normal samples and 412 tumor samples, whereas the GEO (GSE13507) database contains 165 samples with clinical data (Kim et al., 2010). Combining the TCGA-BLCA dataset with the GEO dataset, we eliminated batch effects using the Combat algorithm prior to subsequent analysis. Somatic mutation data were downloaded from TCGA, which contained 407 BLCA samples. Immunohistochemical (IHC) images of proteins taken from Human Protein Atlas (https://www.proteinatlas.org/) (Uhlén et al., 2015).

Differential expression analysis of MTRGs and clinical correlation of molecular subtypes

We examined the variance in expression among 41 MTRGs and explored the association between overall survival (OS) and MTRGs. Using the R “Consensus Cluster Plus” software package, we conducted consensus unsupervised cluster analysis based on the differentially expressed MTRGs. This method strengthens intra-group links while weakening inter-group links. As a result, BLCA patients were classified into different molecular subtypes. To explore the clinical significance of these molecular subtypes, we utilized Kaplan–Meier curves to examine their association with overall survival (OS) and clinicopathological characteristics. Clinicopathological features included age, sex, grade, and TNM stage.

Gene set enrichment analysis and immune profiling analysis

To further investigate the pathways and biological processes of MTRGs molecular subtypes, we conducted a series of enrichment analyses, including GSVA and GSEA analyses. Gene set required for enrichment analysis obtained from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/) (Liberzon et al., 2011). Additionally, we employed the R package “CIBERSORT” to calculate the differences in tumor immune cells (TICs) infiltration between different molecular subtypes (Newman et al., 2015).

Gene subtypes and prognostic model construction

Differentially expressed genes (DEGs) between molecular subtypes based on MTRGs were identified using the R package “limma” (Ritchie et al., 2015). A significance level of P < 0.05 and —logFC—>1.5 was considered meaningful. Following that, we conducted univariate Cox regression analysis to identify DEGs linked to prognosis. Based on consensus cluster analysis of the expression of prognostic DEGs, patients were stratified into different gene subtypes. We further investigated the clinical significance of these gene subtypes and differences in MTRGs expression. The TCGA-BLCA and GSE13507 datasets were divided into training and test sets at a ratio of 1:1 using the R package “caret”. MTRG risk scores were constructed using logistic minimum absolute contraction and selection operator (LASSO) and COX regression in the training set and validated in the test and merged sets. The MTRG risk score was calculated as the sum of the product of the risk coefficient (coef) and gene expression (Exp) for each prognostic model gene, where coef represents the risk coefficient and Exp represents the expression of the prognostic model gene. We visualized risk scores and prognostic outcomes using mulberry graph in relation to different molecular or gene subtypes. The training, test, and merged sets were all stratified into high-risk and low-risk groups based on median risk scores. Survival analyses, heat maps, and receiver operating characteristic (ROC) curves were generated for each group. Additionally, using the R package “rms” and significant clinicopathological features (p < 0.05), we developed a nomogram incorporating age, sex, TMN stage, and MTRGs risk score. This nomogram was used to estimate short-, medium-, and long-term overall survival in patients with BLCA. Calibration plot of the prognostic model was generated to verify the accuracy of the model.

Tumor microenvironment (TME), tumor mutation burden (TMB) and drug sensitivity analysis

We assessed the expression of MTRGs in high- and low-risk groups and calculated differences in the TME across risk groups. Correlation analysis further revealed the relationship between TICs and the eight scoring genes. Mutation frequencies in different risk groups were assessed using the R package “maftools”. To investigate differences in drug susceptibility between risk groups, we calculated the half-inhibitory concentration (IC50) values for common therapeutic drugs using the R package “pRRophetic”. The method uses drug sensitivity data of Cancer Cell Line Encyclopedia (CCLE) in combination with gene expression data of BLCA samples to estimate 50% inhibitory concentration based on prediction model.

Immunohistochemistry (IHC) and single cell analysis

Two genes (HTRA1 and DCBLD2, p < 0.05) were identified as being positively associated with prognosis by analyzing the association between eight prognostic score genes and overall survival (OS). Immunohistochemical (IHC) images from the Human Protein Atlas database (HPA, https://www.proteinatlas.org/) were used to compare HTRA1 and DCBLD2 expression levels between cancer tissues and adjacent tissues. The expression and distribution of HTRA1 and DCBLD2 in various cell clusters were assessed in the BLCA- GSE130001 single-cell dataset using the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org/).

Verification of DCBLD2 expression in BCLA tissue samples

On one hand, immunohistochemistry showed that HTRA1 was highly expressed in BLCA tissues, although not significantly. On the other hand, considering relevant literature in recent years, we ultimately selected DCBLD2 for further validation. Ten pairs of tumor and adjacent tissues were obtained from BLCA patients who recently underwent radical surgery at the First Affiliated Hospital of Nanchang University. All participants provided written informed consent. Total RNA extraction was performed using TRIzol reagent, followed by cDNA synthesis using the Takara PrimeScript RT kit. Real-time quantitative PCR was conducted using SYBR Green (Roche, Basel, Switzerland) reagent. Relative gene expression levels were determined using the 2$\stackrel{ˆ}{}$-ΔΔCt method, with β-actin as the endogenous reference. The primer sequences used were as follows:

DCBLD2_F: CCTGCAAAAGCAGTGGACCATG.

DCBLD2_R: CTCCTACCAGTGGCTGAGCATA.

Cellular cultivation and transfection

Human BLCA cell lines BIU and T24 were cultured in RPMI-1640 and DMEM media, respectively. All cells were incubated in media supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin in a humidified incubator at 37 °C and 5% CO2. These cell lines were obtained from Procell Life Science & Technology Co., Ltd (Wuhan, China) and authenticated by the cell bank of the Chinese Academy of Sciences Type Culture Collection Center (Shanghai, China). siRNA targeting DCBLD2 and controls were synthesized by GenePharma, Ltd (Shanghai, China). Cells were seeded into 6-well plates at 50%–60% density per the manufacturer’s instructions, and siRNA was transiently transfected using Lipofectamine 2000 (Invitrogen). The siRNA sequences used were as follows:

DCBLD2 siRNA-1 (5′-GUUCCUCCUGCUCUUACUUTT-3′)

DCBLD2 siRNA-2 (5′-GGCCUCAUACUCUGUUAUATT-3′)

Wound healing and transwell assay

T24 and BIU cells were seeded into six-well plates at 50% density and allowed to grow to near confluence after transfection with siRNA. A cross was scratched onto the dish using a 200 µl pipette tip, and the debris was washed away with PBS. After changing the medium to serum-free medium, images were captured at 0 h under a microscope. The six-well plate was then placed in the incubator for an additional 24 h, and images were taken again after 24 h. The wound healing rate was quantified using ImageJ software.

For the Transwell assay, 40,000 transfected T24 and BIU cells were mixed with 200-µl of serum-free medium and evenly inoculated into the upper chamber. The lower chamber contained 800 µl of medium containing 20% FBS. After incubation for 24 h, cells were fixed with 4% paraformaldehyde and stained with 1% crystal violet. Cells on the lower surface of the chamber were photographed under a light microscope and counted to measure the extent of cell migration.

Plate cloning and EDU experiments

T24 and BIU cells transfected with siRNAs targeting DCBLD2 and controls were incubated for 24 h, seeded at 1,000 per well into 6-well plates, and cultured in medium containing 10% FBS. After 14 days, they were removed and fixed with 4% paraformaldehyde for 30 min and stained with 1% crystal violet for 30 min. Transfected cells were inoculated into 96-well plates at 5,000 cells per well and incubated for 24 h, then the medium was replaced with EDU reagent and incubated for another 2 h. Subsequent staining was then performed according to the instructions provided in the Edu assay kit (C10310-1; RiboBio), and images were taken under a fluorescence microscope after all staining.

Statistical analysis

All bioinformatics results were analyzed using R 4.3.2 software, and experimental data were visualized using GraphPad Prism 9.5.1. Student t-test was used to compare means between groups, while one-way/two-way ANOVA was used to compare means between groups. For multiple hypothesis testing, the Benjamini–Hochberg (BH) method was used to adjust the p-value to control for false discovery rate (FDR). Adjusted p-values <0.05 were considered statistically significant.

Ethics approval and consent to participate

All organizations included in this study obtained ethical approval from the Ethics Committee of the First Affiliated Hospital of Nanchang University (Approval ID: (2024) CDYFYYLK (08-040)), and patient participation was contingent upon informed consent.

Differentially expressed MTRGs and molecular subtypes

To determine changes in the expression of MTRGs in BLCA patients, we analyzed the differences in the expression of 41 MTRGs in the TCGA-BLCA dataset. Among them, 24 MTRGs showed different expression levels between tumor and normal tissues (p < 0.05, Fig. 1A). The TCGA-BLCA dataset was merged with the GSE13507 dataset, and the Combat algorithm was utilized to eliminate batch effects. Subsequent analyses were conducted on the merged dataset. Survival analysis revealed that 23 MTRGs were significantly associated with overall survival (OS) in BLCA patients (Figs. S1 and S2). Following this, one-way COX regression analysis indicated that 9 MTRGs were associated with survival in BLCA patients (Table S2). Crossover analysis of survival and univariate Cox regression showed that all 9 genes were significantly associated with prognosis in BLCA patients. Next, a plasma membrane tension network was performed to fully demonstrate the association between MTRGs and their prognostic value in BLCA patients. The thickness of the link between two genes represents the strength of the correlation. The results show that there are close and universal interactions between MTRGs.

To further verify the association of MTRGs with the disease, we employed a consensus clustering algorithm and found that K = 2 was the optimal choice (high within-group correlation, low between-group correlation, Figs. 1C–1D). Consequently, BLCA patients were classified into two molecular subtypes, with patients in subtype B showing a better prognosis than those in subtype A (p < 0.05, Fig. 1E). Principal component analysis (PCA) unveiled substantial disparities between subtypes A and B (Fig. S3A). Additionally, we analyzed the association of MTRGs with clinicopathological characteristics according to different subtypes, including gender, age, TNM stage, grade, and dataset origin (Fig. 2A, clustering details: Rows and columns were hierarchically clustered using the Euclidean distance and complete linkage method.).

Association between TME and molecular subtypes

To explore the characteristics of different subtypes of the TME, we conducted enrichment analyses using KEGG GSVA and GSEA methods. KEGG GSVA analysis revealed that subtype B was significantly enriched in peroxidase and glycerolphospholipid metabolic pathways (Fig. 2B). Conversely, subtype A exhibited significant enrichment in multiple pathways, including natural killer cell-mediated cytotoxicity, JAK-STAT signaling, T cell receptor signaling, MAPK signaling, actin cytoskeletal regulation, and focal adhesion (Fig. 2B). Additionally, employing the R package “CIBERSORT”, we assessed differences in the expression of immune cells between different subtypes. Our analysis indicated notable infiltration differences between subtype A and B, with Activated B cells, Activated CD4 T cells, macrophages, and Regulatory T cells showing higher expression levels in subtype A compared to subtype B (Fig. 2C).

Identification of MTRGs associated gene subtypes

Using the R “limma” software package, we identified 952 DEGs between subtypes A and B (Table S3) and performed GEVA enrichment analysis of the DEGs. GO GEVA enrichment analysis showed that the pathways were mainly concentrated in extracellular matrix components, leukocyte chemotaxis, endoplasmic reticulum lumen, cell adhesion and collagen binding. The KEGG GEVA enrichment assay revealed more cytokine receptors, chemokines, and cell adhesion-related pathways. In conclusion, we hypothesize that DEGs play an important role in the regulation of cell structure and adhesion. Further univariate COX regression analysis of DEGs was used to determine which DEGs were associated with prognosis. Finally, we screened out 447 DEGs (Table S4). Subsequently, we performed consensus cluster analysis of 447 prognostic relevant DEGs (Fig. 3A, Fig. S3B). Based on the results of the analysis, we reclassified the patients into two gene subtypes (gene subtype A, gene subtype B). Similarly, OS differed between gene subtypes, with patients with gene subtypes B having a better prognosis than those with gene subtypes A (p < 0.05, Fig. 3B). We reassessed the association of the three (gene subtype, clinicopathological features, and molecular subtype) as shown in the heat map (Fig. 3C). We also assessed the differences in MTRGs expression across gene subtypes and showed that a total of 29 MTRGs showed differences across gene subtypes (p < 0.05, Fig. 3D).

Identification of MTRG scoring system and nomogram

The combined TCGA-BLCA and GSE13507 dataset (n = 531) was split into a training dataset (n = 266) and a test dataset (n = 265) at a 1:1 ratio using the R package “caret”. Subsequently, LASSO and multivariate COX regression (Figs. 3E–3F) were performed on the 447 prognosis-related DEGs in the training set. Consequently, eight risk score genes (HTRA1, GOLT1A, DCBLD2, UGT1A1, FOSL1, DSC2, IGFBP3, and TAC3) were obtained. Based on the risk scores of these eight genes, we constructed the MTRG risk scoring system as follows: (1)${\displaystyle \text{MTRG Risk Score}=\sum _{i=1}^{n}Ex{p}_i{\beta }_i.}$

In the given equation, “Exp” and “β” denote the adjusted expression value and regression coefficient of gene i, respectively. Utilizing patient clinical data, we constructed a mulberry graph to depict the correlation between molecular subtype , gene subtype , risk group, and final outcome (Fig. 3G).

Similarly, we calculated risk scores for each BLCA sample for different molecular subtypes and gene subtypes. The risk scores of molecular subtype A and gene subtype A were significantly higher than those of molecular subtype B and gene subtype B (Fig. S3C), respectively, consistent with the previous Kaplan–Meier (KM) curve (Figs. 1E and 3B). Examining the expression of MTRGs between the high-score and low-score groups, we noted that 11 MTRGs were upregulated in the low-score group, while 18 MTRGs were upregulated in the high-score group (Fig. 4A).

To validate the accuracy of the scoring system, we employed K–M curves, ROC curves, heat maps, risk score distributions, and survival scatter plots. Validation was performed separately in the merged set, training set, and test set. The Kaplan–Meier curves showed that patients with lower scores showed better OS (Figs. 4B–4D). ROC curves indicated that the MTRG risk score predicted 1-, 3-, and 5-year survival with high sensitivity and specificity (Figs. 4E–4G). Heatmaps depict differences in scoring genes among different scoring groups, with HTRA1 and DCBLD2 significantly expressed in the high-risk group (Figs. 4H–4J). Scatter plots of risk score distribution and survival showed that higher risk scores corresponded to shorter survival times (Fig. S3D).

Consistent results across the merged set, training set, and test set validated the scoring system’s high accuracy. By combining relevant clinical characteristics and the risk scoring system, we constructed a nomogram (Fig. 4K) to assess OS at 1, 3, and 5 years in BLCA patients. Through calibration curve analysis (Fig. 4L), we confirmed the predictive ability of our constructed nomogram credible.

Correlation analysis of MTRG risk score with TME and TMB

To further explore the relationship between the MTRG risk score system and TME, We identified eight risk score genes associated with TIC abundance. As shown in Fig. 5A, eight genes are strongly associated with most TICs, especially T cells regulatory (Tregs). Immunocyte correlation analysis showed that macrophages M0, activated mast cells and neutrophils were positively correlated with risk score, while monocytes, T cells CD8 and regulatory T cells (Treg) were negatively correlated with risk score (p < 0.05, Fig. S3E). The high score group showed higher levels of StromalScore, ImmuneScore, and ESTIMATEScore than the low score group (p 0.001, Fig. 5B). Hence, we speculate that the MTRG risk score could be linked to the TME of BLCA. Some studies have shown that tumor stem cells (CSC) are the fundamental factors of tumorigenesis, drug resistance, recurrence and metastasis, and also an important reason for tumor treatment failure (Patel, Oh & Galsky, 2020; Dobruch et al., 2016). Therefore, we further explored the correlation between the MTRG scoring system and the CSC index. The results showed a negative correlation between score and CSC index (R = − 0.13, P < 0.001, Fig. 5C). Mutation frequencies in different risk groups were analyzed by the R “maftools” program. The analysis revealed that the mutation frequencies of TP53, KMT2D, MUC16, SYNE1, STAG2, ELF3, RB1, RYR2, KMT2C, MACF1, OBSCN, and CSMD3 were elevated in the high-score group (Figs. 5D–5E). However, further analysis of TMB reflecting the number of tumor mutations found no clinically significant difference in TMB between the high and low scoring groups (R = 0.039, p = 0.46, Figs. S4B–S4C).

Identification of drug susceptibility and scoring gene expression levels

Using the R “pRRophetic” software package, we computed IC50 values for several common drugs utilized in BLCA treatment. A lower IC50 value indicates that the sample is more sensitive to a particular drug. We combined this information with the results of the high-risk and low-risk groups to assess the potential therapeutic effect of different drugs in each risk group. The findings indicated that certain drugs such as cisplatin, paclitaxel, doxorubicin, and docetaxel exhibited lower IC50 values in the high-scoring group (Figs. 5F–5I), suggesting lower drug resistance in the high-risk group. To delve deeper into the biological function of risk score genes, we individually analyzed the association of the eight score genes with OS. The analysis revealed a negative correlation between DCBLD2 and HTRA1 expression and prognosis (p < 0.05, Figs. 5J–5K, Fig. S4A). The remaining six genes showed no clinically significant impact on OS and were therefore excluded from subsequent analyses. BLCA samples were queried in the Human Protein Profile (HPA) database to access immunohistochemical images of DCBLD2 and HTRA1 (Fig. 6A). The results indicated that DCBLD2 exhibited moderate to high expression in BLCA compared to normal bladder tissue, whereas HTRA1 showed lower expression levels in BLCA relative to DCBLD2. To further validate the distinct expression of DCBLD2 and HTRA1 at the single-cell level, we utilized a single-cell dataset (BLCA-GSE130001) from the TISCH database to explore gene expression differences among different subtypes. The analysis revealed that DCBLD2 was predominantly expressed in fibroblasts and myofibroblasts, while HTRA1 was primarily upregulated in endothelial cells and fibroblasts (Figs. 6B–6D). Figure 6E shows the composition of various cell types, and the correlation between DCBLD2 and HTRA1 is shown in Fig. S3D.

Cell identification in vitro

PCR analysis of DCBLD2 using 10 tissue pairs confirmed its significant expression in BLCA tissues (Fig. 7A). Also, using PCR, we verified the knockdown efficiency of Si1-DCBLD2 and Si2-DCBLD2 in T24 and BIU cell lines (Fig. 7B). Through wound healing experiments and Transwell assays, we found that decreased DCBLD2 levels was associated with reduced migration and invasion of T24 and BIU cells (Figs. 7C–7D). However, colony formation and EDU experiments showed contrasting results. The findings indicated that decreased DCBLD2 expression did not affect the proliferation of T24 and BIU cells (Figs. 7E–7F). Therefore, based on the results of in vitro cell experiments, we can infer that DCBLD2 is pivotal in facilitating BLCA cell migration and invasion but does not notably impact proliferation.