A tumor microenvironment preoperative nomogram for prediction of lymph node metastasis in bladder cancer

Background: Growing evidence suggests that tumor metastasis necessitates multi-step microenvironmental regulation. Lymph node metastasis (LNM) influences both pre- and post-operative bladder cancer (BLCA) treatment strategies. Given that current LNM diagnosis methods are still insufficient, we intend to investigate the microenvironmental changes in BLCA with and without LNM and develop a prediction model to confirm LNM status.

Method: "Estimation of Stromal and Immune cells in Malignant Tumors using Expression data" (ESTIMATE) algorithm was used to characterize the tumor microenvironment pattern of TCGA-BLCA cohort, and dimension reduction, feature selection, and StrLNM signature construction were accomplished using least absolute shrinkage and selection operator (LASSO) regression. StrLNM signature was combined with the genomic mutation to establish an LNM nomogram by using multivariable logistic regression. The performance of the nomogram was evaluated in terms of calibration, discrimination, and clinical utility. The testing set from the TCGA-BLCA cohort was used for internal validation. Moreover, three independent cohorts were used for external validation, and BLCA patients from our cohort were also used for further validation.

Results: The StrLNM signature, consisting of 22 selected features, could accurately predict LNM status in the TCGA-BLCA cohort and several independent cohorts. The nomogram performed well in discriminating LNM status, with the area under curve (AUC) of 75.1% and 65.4% in training and testing datasets from the TCGA-BLCA cohort. Furthermore, the StrLNM nomogram demonstrated good calibration with p >0.05 in the Hosmer-Lemeshow goodness of fit test. Decision curve analysis (DCA) revealed that the StrLNM nomogram had a high potential for clinical utility. Additionally, 14 of 22 stably expressed genes were identified by survival analysis and confirmed by qPCR in BLCA patient samples in our cohort.

Conclusion: In summary, we developed a nomogram that included an StrLNM signature and facilitated the preoperative prediction of LNM status in BLCA patients.

1 Introduction

Bladder cancer (BLCA) is one of the most common malignant tumor of the urinary tract (1). The incidence and mortality of BLCA are increasing due to factors such as increasing age, environmental pollution, smoking, and other factors (2). Lymph node metastasis (LNM) is a major means of metastasis of BLCA and several theories have been proposed to explain the mechanism of LNM in BLCA, including lymphangiogenesis, epithelial-mesenchymal transition (EMT), cell invasion, and tumor microenvironment (3). Previous studies have shown that LNM status is an independent unfavorable prognostic factor in BLCA patients (4–6), and the five-year overall survival (OS) is considerably lower in patients with LNM than in patients without LNM (28% vs. 54%) (7, 8).

Given that there is such a remarkable difference in prognosis between LN-positive and LN-negative patients, In clinical practice, preoperative assessment of LN status is important. Current diagnostic procedures based on traditional radiology, including computerized tomography (CT) and magnetic resonance imaging (MRI).Though the specificity of diagnosis LNM by these method is (91.9% for CT and 89%–98% for MRI), but the sensitivity is relatively low (36.9% for CT and 76–83% for MRI) to preoperatively predict the lymphatic metastasis of BLCA patients (9, 10). Various groups have attempted to develop a nomogram based on CT as well as other radiomics to predict the status of LNM in BLCA patients. For example, Wu et al. established a nomogram based on radiomics features extracted from arterial-phase CT images. However, this nomogram lacked external validation and was constructed without the use of genetic markers, thus its clinical significance needed to be confirmed further (11). As for positron emission tomography (PET), although PET/CT technology plays an increasingly important role in the preoperative evaluation of clinical patients, both F18-PET/CT (12) and 11C-choline PET/CT (13) were not doing well in the sensitivity of lymph node assessment. What’s more, considering the gap between economic development and medical resources in different countries and regions, it is unrealistic to routinely use PET/CT as a preoperative auxiliary examination for BLCA patients. Thus, it is necessary to establish a new noninvasive and effective preoperative assessment to predict LNM in patients.

Furthermore, the interaction between BLCA cells and other cells in the TME contributes to lymph node metastasis. Recognizing that the TME is characterized in cancer development has changed our understanding of cancer development from a cancer cell-centric perspective to tumor growth and metastasis supporting advanced tumor ecosystems (14). Rather than working alone, tumor cells build particular TME through intimate interactions with the extracellular matrix (ECM) as well as stromal cells (15). Lymphatic endothelial cells could express CCL21 and attract CCR7-expressing tumor cells to lymphatic vessels and present antigens (16). Numerous immune and non-immune cells in the TME infrastructure, along with the factors they secrete, contribute to chronic inflammation, immunosuppression, and the gastrointestinal milieu (17).

With the rapid development of next-generation sequencing technology and the comprehensive study of human transcriptome and gene changes, it has become possible to discover the differences of BLCA in LNM from the genomic level. In recent decades, we have seen a shift from a single analysis of multiple biomarkers to the combinatorial analysis of a panel of biomarkers to construct markers, which is regarded as an useful and effectivetool for clinical management methods (18). Our team previously built a nomogram combining an EMT-LNM signature and somatic genetic mutations for prediction of LNM in BLCA which was pretty useful for clinical practice (19). And in the current study, we focused on TME patterns and tried to establish a nomogram which was centered on an StrLNM signature for preoperative prediction of LNM status in BLCA patients.

2 Material and methods

2.1 Data collection and processing

The TCGA-BLCA dataset was obtained from the TCGA Genomic Data Commons (GDC) (https://portal.gdc.cancer.gov/) and was used as the training and internal validation cohort (Table S3, S4). Since multiple ENSEMBL IDs mapped to a single gene symbol, the highest expressed ENSEMBL ID was used. TCGA-BLCA fragments per kilobase million (FPKM) values were downloaded from the TCGA GDC and transformed into transcripts per kilobase million (TPM) (20). Details of the clinicopathological features of each dataset were obtained in our previous study (21). We included a total of 256 samples in our study including 128 lymph node metastases (LN+) and 228 non-lymph node metastases (LN-)as the total TCGA-BLCA cohort. Three distinct BLCA cohorts, GSE13507 cohort (22), GSE31684 cohort (23), and GSE106534 cohort (24), were acquired for external validation. “Affyy” R packages was used in R for the log2 transformation, background correction, annotation, and quantile normalization to process the raw data (25).

2.2 Identification of differentially expressed genes (DEGs)

The differential expression genes of LN+ and LN- samples was evaluated, as well as high Stromal score samples and low Stromal score samples, in the TCGA-BLCA cohort by R package “DESeq2”. The significant cutoff for defining DEGs was set as |log2Fold Change (FC)| >1.0 and false discovery rate (FDR) < 0.05. Among all the DEGs both upregulated or downregulated genes in LN+/LN- groups and high/low Stromal groups were defined as the final DEGs.

2.3 Generation of StrLNM signature

Primary predictive features was built by the Least absolute shrinkage and selection operator (LASSO) logistic regression analysis, and screened final DEGs was used to build an StrLNM signature by the R package “glmnet” (26). The best value of the penalty parameter λ was determined through 10-fold cross-validation error for dimension reduction to reduce noise or redundant genes. The risk value for each sample’s StrLNM signature was calculated from a linear combination of selected features, weighted with an appropriate coefficient. The equation of risk value=${\sum }_{i=1}^n$ (coefi × Expri), where Expri is the relative expression of the genes in the signature for patients i and coefi was the relevant coefficients of the genes. And we further used U test to explore the correlation of StrLNM and LNM in testing dataset.

2.4 Development of an individualized prediction model

Univariable and multivariable logistic regression analysis were further used to development the model with candidate features including the StrLNM signature and C3orf70. StrLNM nomograms were generated using the “rms”, “nomogrammex”, and “regplot” R packages as a quantitative tool for clinicians to predict individual LNM probabilities.

2.5 Validation of StrLNM signature and nomogram

After the StrLNM signature was established, internal validation was performed on the TCGA-BLCA cohort test dataset according to the risk assessment formula specified in the training dataset. The nomogram was then tested for predictive accuracy and stability using receiver operating characteristic (ROC) curves and calibration curves by using R packages “pROC” and “rms”. Moreover, three independent GEO cohorts (GSE13507, GSE31684, and GSE106534) were used for the external validation. A supervised hierarchical clustering method was used to explore clusters with k = 2 based on 1-Pearson’s correlation distance.

Tissue specimens

As previously described, a total of 55 BLCA specimens were randomly recruited and paired adjacent normal tissue samples were obtained at the Beijing Friendship Hospital, Capital Medical University (Beijing, China) between January 2021 and March 2022. The clinical characteristic of our specimens was listed in Table 1. And there was no attrition during the study.

TABLE 1

Table 1 Characteristic of patients.

The inclusion criteria for our specimens were patients with a clinical and histological diagnosis of BLCA.And the exclusion criteria were patients who had received any chemotherapy, neoadjuvant chemotherapy, radiotherapy and immunotherapy. Patients with severe urinary tract infection, severe renal failure, serious underlying diseases of other systems, or who had participated in other studies were also excluded. This study was approved by the Ethics Committee of the Affiliated Friendship Hospital of Capital Medical University (NO.2021-P2-159). The patients provided their written informed consent to participate in this study.

2.6 RNA extraction, reverse transcription, and quantitative real-time PCR (qRT-PCR)

During the whole process of RNA extraction and PCR, operators follow the principle of blind method. Total RNA of clinical samples was extracted by the RNeasy plus mini kits (74136, Qiagen, Germany) Subsequently, NanoDrop instrument (NP80, Implen, Germany) was used to exam the quality of the extracted RNA. Then, the RNA was used as a template for cDNA synthesis using the ReverTra Ace qPCR RT Kit (FSQ-201, Toyobo, Japan). Finally, forward and reverse primers of our key genes and iQTM SYBR^® Green Supermix (1708880, Bio-Rad, US) were mixed, and performed the qRT-PCR. The primer sequences were listed in Table S1. The relative expression level of the targeted genes was normalized to GAPDH as described by Schmittgen TD et al. 2-△△CT=[(CT gene of interest - CT internal control) sample A - (CT gene of interest - CT internal control) sample B] (27).

2.7 Statistical analyses

The statistical significance of data between two groups was tested by Student’s t-test, Mann-Whitney U test, Fisher’s exact and χ2 tests according to the data type. We used R package “survival” and “survminer” to generate Kaplan-Meier curves and Cox regression for survival analysis (28). The significant difference between survival curves belonging to different defined groups were determined with the log-rank (Mantel-Cox) test. And we used R packages “rms,” “nomogramEx,” and “regplot” to build nomogram and calibration curves (21). Decision curve analysis (DCA) was performed to determine whether our established nomogram was of clinical usefulness according to Iasonos et al.’s suggestion (29). The package “pROC” in R was used to plot and visualize ROC curves. All statistical analyses were performed with R software 3.5.3. Statistical cutoff point was set at p < 0.05.

3 Results

3.1 Characterization of patients in BLCA

We first assessed Stromal score, Immune score, and Estimated score in 356 patients with or without LNM in the TCGA-BLCA cohort (Table S2). Stromal score were significantly higher in patients with LNM (Figure 1A, p = 0.0022); However, Immune score and Estimated score did not differ significantly between these patients (Figures 1B, C). And we further examined overall survival in patients with or without LNM, and it is clear that LNM is an unfavorable prognostic factor (Figure 1D). Also, patients with a higher Stromal score had a potential for poorer overall survival (OS) (Figure 1E). Therefore, we divided all patients into four groups according to the Stromal score and the LNM status. Interestingly, regardless of the level of Stromal score, Kaplan-Meier survival curves showed significantly improved survival in LN- patients, and patients without LNM and low Stromal score level had a better prognosis than other groups (log-rank test, p < 0.05, Figure 1F). And the Recurrence free survival (RFS) curve was similar to OS and was shown in Supplementary Figure S1.

FIGURE 1

Figure 1 Association between LN metastasis status and TME related score of BLCA patients in the TCGA-BLCA cohort for overall survival (OS). (A) Correlation between LN metastasis and Stromal score. (B) Correlation between LN metastasis and immune score. (C) Correlation between LN metastasis and ESTIMATE score. (D) Kaplan-Meier survival curves of LN metastasis (E) Kaplan-Meier survival curves of Stromal score (F) Kaplan-Meier survival curves show the difference in prognosis advantage among four groups stratified by LN metastasis and Stromal score.

3.2 DEGs between LN+/- and high/low Stromal score patients

When comparing LN+ and LN-, we screened 121 down-regulated DEGs (DEG_LN_DN) and 230 up-regulated DEGs (DEG_LN_UP) by the package “DESeq2” in R (Figure 2A). In addition, we identified 3552 down-regulated differentially expressed genes (DEG_stromal_DN) and 1096 up-regulated differentially expressed genes (DEG_stromal_UP) between patients with high and low Stromal score (Figure 2B). Among all these DEGs there were 48 DEGs both involved in DEG_stromal_DN and DEG_LN_DN, and 105 DEGs participated in both DEG_stromal_UP and DEG_LN_UP (Figure 2C). Therefore, we defined these 153 genes as the most important genes among LN+/LN- patients. The heatmap showed the expression levels of 153 DEGs in patients was characterized by LN+/LN- and Stromal score (Figure 2D).

FIGURE 2

Figure 2 DEGs between LN+/- and high/low Stromal score patients. (A) DEGs between LN+ and LN-. (B) DEGs between high/low Stromal score patients. (C) DEGs both involved in Stromal score and LN. (D) Heatmap of DEGs both involved in Stromal score and LN.

3.3 Establishment of StrLNM-Signature

We further performed LASSO regression analysis on 153 DEGs and constructed a signature (StrLNM signature) to distinguish LNM status in the TCGA-BLCA cohort training dataset (Figures 3A, B). The coefficient of each feature in the StrLNM signature was shown in Table S5. Furthermore, the StrLNM signature was significantly higher in LN+ tumors than in LN tumors in both training(p <0.001) and test datasets(p = 0.04) of the TCGA-BLCA cohort, and the entire dataset also matched well with the same results (Figures 3C–E).

FIGURE 3

Figure 3 Feature selection using LASSO logistic regression model. (A) Tuning parameter (λ) selection 10-fold cross-validation error curve. The misclassification error was plotted vs. log (λ). The dotted vertical lines were drawn at the optimal values by the minimum criteria and the 1-SE criteria. (B) LASSO coefficient profiles of the 19 candidate EMT-related genes. A coefficient profile plot was produced against the log (λ) sequence. A vertical line was drawn at the value selected by 10-fold cross-validation, where the optimal λ resulted in 17 nonzero coefficients. (C–E) Difference in the StrLNM signature risk score between LN+ and LN- tumors in the training dataset (C), testing dataset (D) and entire set (E) of the TCGA-BLCA cohort.

3.4 Supervised hierarchical clustering according to StrLNM signature

In addition, we performed hierarchical clustering to see whether StrLNM signature could discriminate LNM in the other three external cohorts. We found cluster C1 was mainly concentrated in LN+ tumors, whereas the cluster C2 gathered in LN- patients in the GSE106534 cohort (Figure 4A). When we further explored the StrLNM signature score in LN+ and LN- patients, it was obvious that LN+ patients had a higher StrLNM signature score (Figure 4B). The area under curve (AUC) of the StrLNM signature score in the GSE106534 cohort was 0.88 (Figure 4C). We then further verified the StrLNM cluster in the other three datasets (GSE13507, GSE31684 and GSE48075) and confirmed that cluster C1 which almost all consisted of LN+ had a significantly worse OS and cancer-specific survival (CSS) than cluster C2 patients in GSE13507(Figures 4D–F), in GSE31684 (Figures 4G–I), and in GSE48075 (Figures 4J–L). These data suggested that stromal-rich BLCA patients were prone to have LNM and worse outcomes.

FIGURE 4

Figure 4 Supervised hierarchical clustering for EMT-LN signature. (A) Dendrogram showed that two clusters created by supervised hierarchical clustering could significantly distinguish LN metastasis status in the GSE106534 cohort. (B) different StrLNM signature score in LN+ and LN- patients. (C) ROC curves of StrLNM signature score in the GSE106534 cohort. (D) Dendrogram showed that two clusters created by supervised hierarchical clustering were strikingly associated with OS, RFS, and CSS in the GSE13507 cohort. (E, F) Kaplan-Meier survival curves showed the difference in OS (log-rank test, p =0.006, E) and CSS (log-rank test, p < 0.001, F) between LN+ and LN- tumors in the GSE13507 cohort. (G) Dendrogram showed that two clusters created by supervised hierarchical clustering were strikingly associated with disease metastasis survival (DMS), lymph node recurrence survival (LRS), and urinary tract recurrence survival (URS) in the GSE31684 cohort. (H, I) Kaplan-Meier survival curves showed the difference in OS (log-rank test, p =0.088, H) and RFS (log-rank test, p =0.056, I) between LN+ and LN- tumors in the GSE31684 cohort. (J) Dendrogram showed that two clusters created by supervised hierarchical clustering were strikingly associated with P53, RAS, RB1 and FGFR3 mutation as well as OS and CSS in the GSE48075 cohort. (K, L) Kaplan-Meier survival curves showed the difference in OS (log-rank test, p < 0.001, K) and CSS (log-rank test, p < 0.001, L) between LN+ and LN- tumors in the GSE48075 cohort.

3.5 Individualized prognostic nomogram: Development and validation

Logistic regression analysis was used to determine the predictive characteristics of preoperative lymph node metastasis. Considering that the pathological staging and detailed classification of TNM are more obtained by postoperative pathological examination than preoperative evaluation, in this study we included the StrLNM signature and mutation of C3orf70, which were proved to be obviously highly mutated in LN+ tumors when compared with LN- tumors (Table S6) (19), as the candidate features to build the preoperative nomogram.

Logistic regression found the StrLNM signature and C3orf70 mutation were statistically significant (p < 0.05) (Table S6). Therefore, we synthesize these features to form the StrLNM nomogram (Figure 5). The scores for each prognostic parameter for each patient were summed to obtain a complete score according to the nomogram. The higher the complete score, the more chance of LNM. The nomogram calibration curve then shows that the model performs similarly to the ideal model in the training set (p = 0.737) (Figure 6A). ROC curves (Figure 6B) showed that the StrLNM nomogram could effectively predict the LNM, with high AUCs in both training dataset (AUC: 0.751 [0.696−0.803]) and testing dataset (AUC: 0.654 [0.462−0.813)). Then, DCA curve showed that the net benefit of the StrLNM nomogram was higher than that of the “complete treatment” or “no treatment” strategy, suggesting that the StrLNM nomogram had a higher potential clinical application value (Figure 6C).

FIGURE 5

Figure 5 Development of preoperative StrLNM nomogram. By incorporating the StrLNM signature and genomic mutation of C3orf70, the StrLNM nomogram was built in the training dataset of the TCGA-BLCA cohort.

FIGURE 6

Figure 6 Clinical utility of the StrLNM nomogram. (A) Calibration curve of the StrLNM nomogram in the training dataset of the TCGA-BLCA cohort, which depicts the calibration of the fitted model between the predicted risk of LN metastasis and actual LN metastasis rate. The x-axis represents the predicted LN metastasis risk, and the y-axis represents the actual LN metastasis rate. The pink solid line represents the performance of the StrLNM nomogram, of which a closer fit to the diagonal dotted blue line represents an ideal prediction. The difference in the two models was measured with the Hosmer-Lemeshow test. (B) ROC curves showed the prediction accuracy of the StrLNM signature in prediction of the LN metastasis in training and testing datasets of the TCGA-BLCA cohort. (C) Decision curve analysis (DCA) for the StrLNM nomogram. The y-axis measures the net benefit. The pink line represents the StrLNM nomogram, the blue line represents the assumption that all patients have LN metastases, and the black line on the bottom represents the assumption that no patients have LN metastases.

3.6 Expression level of StrLNM candidate genes in BLCA patient samples

A total of 55 BLCA specimens and paired adjacent normal tissue samples, including 16 LN+ samples and 39 LN- samples were obtained. We validated the expression of 14 stably expressed candidate genes (over 50% samples detected in the TCGA-BLCA cohort) in BLCA patients’ samples by qPCR and examined the association between these genes and survival in the TCGA-BLCA cohort (Figure 7). AZGP1, C11orf86, CLDN9, CLIC3, EPHB6, ERVV.2, HS3ST2, HSD17B2, KCNK13, and STEAP4 showed a significant correlation with survival while the p-value of SHH was also close to 0.05 (Table S7). We then further explored the expression of these genes between LN+ and LN- BLCA patients, which showed that the expression level of AZGP1, C11orf86, HSD17B2, and SHH in LN+ patients were significantly lower than in LN- patients, meanwhile, CLIC3 and KCNK13 were both highly expressed in LN+ patients.

FIGURE 7

Figure 7 Overall survival of 14 stably-expressed genes. Kaplan-Meier survival curves showed the different correlation between OS and expression level of candidate genes.

4 Discussion

In the present study, we focused on stromal cells of TME and tried to establish a nomogram combining an StrLNM signature for preoperatively predict LNM and confirmed that AZGP1, C11orf86, CLIC3, HSD17B2, KCNK13, and SHH were potentially predictive biomarkers in BLCA. These results indicate that stromal infiltration involvement in LNM and the potential clinical utility of our StrLNM signature.

Cancer research has recently shifted from focusing solely on tumor cells to comprehensively understanding the surroundings of core tumor cells, known as the TME (30–32). Tumorigenesis may be a complex and dynamic process with three stages: initiation, progression, and metastasis, and the physiological state of the TME was found to be closely related to the various stages of tumorigenesis. TME immune cell infiltration in situ has been identified as a piece of important and valuable information for predicting the prognosis and immunotherapy of various cancers in clinical studies of ICIs (33–36). Multiple studies have demonstrated the important role of TME in determining prognosis and LNM in BLCA (24, 33, 37). These studies demonstrated the enormous potential of TME in BLCA prognosis and preoperative evaluation, which was consistent with our findings that higher stromal infiltration was found in LN+ BLCA and associated with worse outcomes.

Multiple studies have found that lymph node dissection can improve the prognosis of BLCA patients at various stages of tumor growth, implying the importance of a more accurate LNM assessment strategy (38, 39). And in our study, to reduce noise and redundant features, we used the LASSO algorithm to combine all features into a single feature, the marker panel, rather than a single predictor selection based on the strength of univariate regression analysis. The StrLNM signature identified in the current study was distributed significantly differently in LN+ and LN- tumors in the TCGA-BLCA training, test cohorts as well as GEO cohorts, indicating our StrLNM signature could be considered as a non-invasive factor for our preoperative prediction of lymph node metastasis. The nomogram and DCA curves were also used to confirm the predictive value of our StrLNM signature. The nomogram is a quantification tool, and people get an overall score by adding points to each nomogram feature, which calculated the contribution of features according to special algorithms. And in this study, our nomogram was pretty accurate in diagnosing lymph node metastases in the TCGA-BLCA training group (75.1%) and test group (65.4%) compared with 36.9% sensitivity for CT (9). Therefore, our StrLNM nomogram could effectively predict the LNM of patients before operation, and provide an effective auxiliary role for follow-up treatment. Furthermore, DCA curves revealed that our nomogram-based decision-making was superior to all patients with and without treatment, with a higher threshold probability and better clinical outcomes. All of these findings suggested that our developed nomogram had a high potential for clinical application.

Moreover, we found the expression levels of AZGP1, C11orf86, HSD17B2, SHH, CLIC3, and KCNK13 were significantly different between our own BLCA patients’ tumor and para-carcinoma tissues using qPCR, and the expression level of AZGP1, C11orf86, HSD17B2 and SHH in LN+ patients were significantly lower which indicates they all had a positive effect on patients’ prognosis. And in our StrLNM-Signature system, these genes contributed negatively to the StrLNM score. Since we had proved that LN+ patients had a higher StrLNM signature score, we believed AZGP1, C11orf86, HSD17B2 as well as SHH were potentially important biomarker candidates that were significantly associated with LNM and would contribute to the prediction of BLCA patients. On other hand, CLIC3, and KCNK13 involved the same result from a different angle. They were both highly expressed in LN+ patients indicating they both had a negative effect on prognosis while these genes contributed positively to our StrLNM score. As a matter of fact, Wang et al. had already reported the correlation between highly expressed HSD17B2 and better prognosis in bladder (40). And though there was little research in BLCA, there were a series of researches indicating the role of AZGP1 in various cancers. Kong B et al. found Zinc α2-glycoprotein (AZGP1, ZAG) acts as a tumor suppressor in pancreatic ductal adenocarcinoma, and its expression is lost due to histone deacetylation (41). The sonic hedgehog (SHH) signaling pathway, an evolutionarily conserved molecular cascade that is primarily involved in the development of the fetal central nervous system, was a well-documented trail along the bladder. And it was reported to be correlated with tumorigenesis, EMT, and BLCA stemness (42–44). Interestingly, as a canonical potassium channel, KCNK13 was rarely described in tumors. We all know that ion channels such as transient receptor potential (TRP) channels have been shown to play an important role in tumorigenesis and invasion in pan-cancer, and KCNK13 channels, another classic cation channels, deserve further study. These six genes depicted in Figure 8 showed a significant correlation in both LNM assessment and tumor/para-carcinoma tissue expression, indicating that these genes play a more important role in LNM and tumorigenesis and warrant more research.

FIGURE 8

Figure 8 Expression Level of StrLNM candidate genes in BLCA patients’ sample. (A) Expression level of AZGP1 between LN+ and LN- samples. (B) Expression level of C11orf86 between LN+ and LN- samples. (C) Expression level of CLIC3 between LN+ and LN- samples. (D) Expression level of HSD17B2 between LN+ and LN- samples. (E) Expression level of KCNK13 between LN+ and LN- samples. (F) Expression level of SHH between LN+ and LN- samples.

As for clinical use of our StrLNM signature, cystoscopy is necessary for the diagnosis of bladder cancer according to the guideline of European association (45). Also, histological evaluation of resected tissue was also recommended for ultimately diagnosis. Thus, in the process of preoperative cystoscopy, the patient’s BLCA tissue can be easily obtained, and the corresponding simple key gene expression level detection can be carried out, and then the score can be directly obtained for clinical reference.

In this study, we first introduced TME into our StrLNM signature, combined a large external dataset to validate the results, and validate our key genes using patient materials. Nevertheless, there was also a limitation in our study. Our research has not been fully integrated with techniques such as radiomics, and practical clinical applications are still a long way off. And at this stage, it is not available for us to carry out in vitro and in vivo experiments to further verify the functions of our key genes.

5 Conclusion

In conclusion, we created an StrLNM nomogram with StrLNM signature, which is an easily accessible and convenient tool to help preoperatively predict LN metastasis in BLCA patients.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

This study was approved by the Ethics Committee of the Affiliated Friendship Hospital of Capital Medical University.(NO.2021-P2-159). The patients/participants provided their written informed consent to participate in this study.

Author contributions

LY and RC made substantial contributions to the conception and design of the research. ZC, CQ, and GW integrated and analyzed the data. ZC wrote the paper. DS, YT, LY, and RC edited the manuscript and provided critical comments. All authors read and approved the final manuscript.

Funding

This work was supported by National Natural Science Foundation of China (Grant No. 81870509) to YT.

Acknowledgments

We acknowledge the Cancer Genome Atlas Research Network for the clinicopathological and genetic alteration data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.1099965/full#supplementary-material

References

1. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, et al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer (2015) 136(5):ES59–86. doi: 10.1002/ijc.29210

CrossRef Full Text | Google Scholar

2. He YT, Li DJ, Liang D, Zheng RS, Zhang SW, Zeng HM, et al. [Incidence and mortality of bladder cancer in China, 2014]. Zhonghua Zhong Liu Za Zhi. (2018) 40(9):647–52. doi: 10.3760/cma.j.issn.0253-3766.2018.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Liu S, Chen X, Lin T. Lymphatic metastasis of bladder cancer: Molecular mechanisms, diagnosis and targeted therapy. Cancer Lett (2021) 505:13–23. doi: 10.1016/j.canlet.2021.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Stein JP, Lieskovsky G, Cote R, Groshen S, Feng AC, Boyd S, et al. Radical cystectomy in the treatment of invasive bladder cancer: long-term results in 1,054 patients. J Clin Oncol (2001) 19(3):666–75. doi: 10.1200/JCO.2001.19.3.666

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Spradling K, Lotan Y, Shokeir A, Abol-Enein H, Mosbah A, Morgan JB, et al. Lymphovascular invasion is associated with oncologic outcomes following radical cystectomy for squamous cell carcinoma of the urinary bladder. Urol Oncol (2016) 34(9):417 e1–8. doi: 10.1016/j.urolonc.2016.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Cha EK, Sfakianos JP, Sukhu R, Yee AM, Sjoberg DD, Bochner BH. Poor prognosis of bladder cancer patients with occult lymph node metastases treated with neoadjuvant chemotherapy. BJU Int (2018) 122(4):627–32. doi: 10.1111/bju.14242

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Karl A, Carroll PR, Gschwend JE, Knuchel R, Montorsi F, Stief CG, et al. The impact of lymphadenectomy and lymph node metastasis on the outcomes of radical cystectomy for bladder cancer. Eur Urol (2009) 55(4):826–35. doi: 10.1016/j.eururo.2009.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Del Bene G, Calabro F, Giannarelli D, Plimack ER, Harshman LC, Yu EY, et al. Neoadjuvant vs. adjuvant chemotherapy in muscle invasive bladder cancer (MIBC): Analysis from the RISC database. Front Oncol (2018) 8:463. doi: 10.3389/fonc.2018.00463

PubMed Abstract | CrossRef Full Text | Google Scholar

9. McKibben MJ, Woods ME. Preoperative imaging for staging bladder cancer. Curr Urol Rep (2015) 16(4):22. doi: 10.1007/s11934-015-0496-8

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Girard A, Rouanne M, Taconet S, Radulescu C, Neuzillet Y, Girma A, et al. Integrated analysis of (18)F-FDG PET/CT improves preoperative lymph node staging for patients with invasive bladder cancer. Eur Radiol (2019) 29(8):4286–93. doi: 10.1007/s00330-018-5959-0

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wu S, Zheng J, Li Y, Yu H, Shi S, Xie W, et al. A radiomics nomogram for the preoperative prediction of lymph node metastasis in bladder cancer. Clin Cancer Res (2017) 23(22):6904–11. doi: 10.1158/1078-0432.CCR-17-1510

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ha HK, Koo PJ, Kim SJ. Diagnostic accuracy of f-18 FDG PET/CT for preoperative lymph node staging in newly diagnosed bladder cancer patients: A systematic review and meta-analysis. Oncology (2018) 95(1):31–8. doi: 10.1159/000488200

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ceci F, Bianchi L, Graziani T, Castellucci P, Pultrone C, Eugenio B, et al. 11C-choline PET/CT and bladder cancer: lymph node metastasis assessment with pathological specimens as reference standard. Clin Nucl Med (2015) 40(2):e124–8. doi: 10.1097/RLU.0000000000000604

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol (2016) 27(8):1482–92. doi: 10.1093/annonc/mdw168

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell (2012) 21(3):309–22. doi: 10.1016/j.ccr.2012.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Xiong Y, Huang F, Li X, Chen Z, Feng D, Jiang H, et al. CCL21/CCR7 interaction promotes cellular migration and invasion via modulation of the MEK/ERK1/2 signaling pathway and correlates with lymphatic metastatic spread and poor prognosis in urinary bladder cancer. Int J Oncol (2017) 51(1):75–90. doi: 10.3892/ijo.2017.4003

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Gajewski TF, Schreiber H, Fu YX. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol (2013) 14(10):1014–22. doi: 10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Birkhahn M, Mitra AP, Cote RJ. Molecular markers for bladder cancer: the road to a multimarker approach. Expert Rev Anticancer Ther (2007) 7(12):1717–27. doi: 10.1586/14737140.7.12.1717

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Cao R, Ma B, Wang G, Xiong Y, Tian Y, Yuan L. An epithelial-mesenchymal transition (EMT) preoperative nomogram for prediction of lymph node metastasis in bladder cancer (BLCA). Dis Markers (2020) 2020:8833972. doi: 10.1155/2020/8833972

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci (2012) 131(4):281–5. doi: 10.1007/s12064-012-0162-3

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Cao R, Yuan L, Ma B, Wang G, Qiu W, Tian Y. An EMT-related gene signature for the prognosis of human bladder cancer. J Cell Mol Med (2020) 24(1):605–17. doi: 10.1111/jcmm.14767

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Kim WJ, Kim EJ, Kim SK, Kim YJ, Ha YS, Jeong P, et al. Predictive value of progression-related gene classifier in primary non-muscle invasive bladder cancer. Mol Cancer (2010) 9:3. doi: 10.1186/1476-4598-9-3

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Riester M, Taylor JM, Feifer A, Koppie T, Rosenberg JE, Downey RJ, et al. Combination of a novel gene expression signature with a clinical nomogram improves the prediction of survival in high-risk bladder cancer. Clin Cancer Res (2012) 18(5):1323–33. doi: 10.1158/1078-0432.CCR-11-2271

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Chen C, He W, Huang J, Wang B, Li H, Cai Q, et al. LNMAT1 promotes lymphatic metastasis of bladder cancer via CCL2 dependent macrophage recruitment. Nat Commun (2018) 9(1):3826. doi: 10.1038/s41467-018-06152-x

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gautier L, Cope L, Bolstad BM, Irizarry RA. Affy–analysis of affymetrix GeneChip data at the probe level. Bioinformatics (2004) 20(3):307–15. doi: 10.1093/bioinformatics/btg405

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01

CrossRef Full Text | Google Scholar

27. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc (2008) 3(6):1101–8. doi: 10.1038/nprot.2008.73

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bland JM, Altman DG. Survival probabilities (the Kaplan-Meier method). BMJ (1998) 317(7172):1572. doi: 10.1136/bmj.317.7172.1572

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol (2008) 26(8):1364–70. doi: 10.1200/JCO.2007.12.9791

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wu T, Dai Y. Tumor microenvironment and therapeutic response. Cancer Lett (2017) 387:61–8. doi: 10.1016/j.canlet.2016.01.043

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Iesato A, Nucera C. Tumor microenvironment-associated pericyte populations may impact therapeutic response in thyroid cancer. Adv Exp Med Biol (2021) 1329:253–69. doi: 10.1007/978-3-030-73119-9_14

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Fu T, Dai LJ, Wu SY, Xiao Y, Ma D, Jiang YZ, et al. Spatial architecture of the immune microenvironment orchestrates tumor immunity and therapeutic response. J Hematol Oncol (2021) 14(1):98. doi: 10.1186/s13045-021-01103-4

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Cao R, Yuan L, Ma B, Wang G, Tian Y. Tumour microenvironment (TME) characterization identified prognosis and immunotherapy response in muscle-invasive bladder cancer (MIBC). Cancer Immunol Immunother (2021) 70(1):1–18. doi: 10.1007/s00262-020-02649-x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol (2020) 17(8):807–21. doi: 10.1038/s41423-020-0488-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Lai Y, Tang F, Huang Y, He C, Chen C, Zhao J, et al. The tumour microenvironment and metabolism in renal cell carcinoma targeted or immune therapy. J Cell Physiol (2021) 236(3):1616–27. doi: 10.1002/jcp.29969

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Toor SM, Sasidharan Nair V, Decock J, Elkord E. Immune checkpoints in the tumor microenvironment. Semin Cancer Biol (2020) 65:1–12. doi: 10.1016/j.semcancer.2019.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Liu J, Zheng Z, Zhang W, Wan M, Ma W, Wang R, et al. Dysregulation of tumor microenvironment promotes malignant progression and predicts risk of metastasis in bladder cancer. Ann Transl Med (2021) 9(18):1438. doi: 10.21037/atm-21-4023

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Kosiba M, Stolzenbach LF, Colla Ruvolo C, Nocera L, Mansour M, Tian Z, et al. Contemporary trends and efficacy of pelvic lymph node dissection at radical cystectomy for urothelial and variant histology carcinoma of the urinary bladder. Clin Genitourin Cancer (2022) 20(2):195e1–e8. doi: 10.1016/j.clgc.2021.10.010

CrossRef Full Text | Google Scholar

39. Bi L, Huang H, Fan X, Li K, Xu K, Jiang C, et al. Extended vs non-extended pelvic lymph node dissection and their influence on recurrence-free survival in patients undergoing radical cystectomy for bladder cancer: a systematic review and meta-analysis of comparative studies. BJU Int (2014) 113(5b):E39–48. doi: 10.1111/bju.12371

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang CT, Li CF, Wu WJ, Huang CN, Li CC, Li WM, et al. High expression of 17beta-hydroxysteroid dehydrogenase type 2 is associated with a better prognosis in urothelial carcinoma of the urinary tract. J Cancer (2016) 7(15):2221–30. doi: 10.7150/jca.16777

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Kong B, Michalski CW, Hong X, Valkovskaya N, Rieder S, Abiatari I, et al. AZGP1 is a tumor suppressor in pancreatic cancer inducing mesenchymal-to-epithelial transdifferentiation by inhibiting TGF-beta-mediated ERK signaling. Oncogene (2010) 29(37):5146–58. doi: 10.1038/onc.2010.258

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Islam SS, Mokhtari RB, Noman AS, Uddin M, Rahman MZ, Azadi MA, et al. Sonic hedgehog (Shh) signaling promotes tumorigenicity and stemness via activation of epithelial-to-mesenchymal transition (EMT) in bladder cancer. Mol Carcinog (2016) 55(5):537–51. doi: 10.1002/mc.22300

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Wang D, Xie D, Bi L, Wang Y, Zou C, Chen L, et al. Benzidine promotes the stemness of bladder cancer stem cells via activation of the sonic hedgehog pathway. Oncol Lett (2021) 21(2):146. doi: 10.3892/ol.2020.12407

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Syed IS, Pedram A, Farhat WA. Role of sonic hedgehog (Shh) signaling in bladder cancer stemness and tumorigenesis. Curr Urol Rep (2016) 17(2):11. doi: 10.1007/s11934-015-0568-9

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Witjes JA, Bruins HM, Cathomas R, Comperat EM, Cowan NC, Gakis G, et al. European Association of urology guidelines on muscle-invasive and metastatic bladder cancer: Summary of the 2020 guidelines. Eur Urol (2021) 79(1):82–104. doi: 10.1016/j.eururo.2020.03.055

PubMed Abstract | CrossRef Full Text | Google Scholar