Wnt pathway-related three-mRNA clinical outcome signature in bladder urothelial carcinoma: computational biology and experimental analyses

The Wnt signaling pathway is core to the growth of bladder tumors. Epithelial-to-mesenchymal transition (EMT) is significant for bladder tumor metastasis. Nevertheless, the relationship between the Wnt signaling pathway, outcomes of bladder cancer (BLCA), and the specific mechanisms driving immune infiltration have not been studied. We obtained Wnt pathway-related gene mRNA and clinicopathological data from the Cancer Genome Atlas (TCGA). We obtained 34 genes that were greatly correlated with outcome using univariate Cox regression analysis and conducted a completely randomized data t-test to perform clinical staging. According to the single-sample gene set enrichment analysis (ssGSEA), the weighted correlation network analysis (WGCNA) was applied to identify relevant biological functions. Various subtypes were identified using consensus cluster analysis. Univariate Cox regression and least absolute shrinkage sum selection operator–Cox regression algorithm analysis were conducted on TCGA and Gene Expression Omnibus data to identify risk characteristics. The Kaplan–Meier method and receiver running feature curves were adopted to calculate overall survival. Single-sample gene set enrichment analysis (ssGSEA) was adopted for the assessment of the degree of immune infiltration. Then, we demonstrated the relationship between PPP2CB and EMT function in two cell lines. Thirty-four Wnt signaling pathway-related genes were risk factors for BLCA outcome, and their expression levels differed by clinical stage. The co-expression of WGCNA showed the relationship between the Wnt signaling pathway and biological functions and was closely associated with EMT. We divided BLCA patients into two subtypes using consensus clustering. Survival curves and clinical analysis showed that the Wnt pathway enriched group had worse outcomes. The Wnt signature showed the significance of the outcome for MAPK10, PPP2CB, and RAC3. Based on these genes, the degree of immune infiltration was evaluated. Cell function experiments suggested that PPP2CB drives the proliferation and migration of BLCA cells. We found that Wnt signaling pathway-related genes can be used as prognostic risk factors for BLCA, and the Wnt signaling pathway is a cancer-promoting signaling pathway associated with EMT. We identified three critical genes: MAPK10, RAC3, and PPP2CB. The genes in these three Wnt signaling pathways are associated with tumor cell EMT and immune cell infiltration. The most important finding was that these three genes were independent prognostic factors for BLCA.


Introduction
Bladder cancer (BLCA) is the tenth most usual cancer and ranks sixth among males. It is associated with high incidence and mortality rates. There are approximately 549,393 new cases and 199,922 deaths in the world annually [1]. The Wnt signaling pathway is a relatively classic signaling pathway related to cell differentiation, migration and cell polarity. The Wnt signaling pathway is fallen into three sub-pathways: the classic β-catenin-dependent path, the non-canonical Wnt/calcium path and the non-canonical planar cell polarity path. Abnormal Wnt signaling is related to some cancers, the most remarkable ones being colorectal, breast, lung, breast, oral, and cervical cancer and hematopoietic malignancies [2]. Most studies of the Wnt signaling pathway in BLCA only focused on the Wnt pathway as a downstream signaling pathway under the regulation of upstream molecules or medications to regulate the proliferation and differentiation of BLCA cells. These studies did not mention the influence of the Wnt signaling pathway on BLCA outcomes. Therefore, the current research determined the relationship between the Wnt signaling pathway and BLCA outcome and identified independent prognostic factors and the effects of Wnt pathway-related genes on BLCA outcome by using RNA-seq data downloaded from The Cancer Genome Atlas (TCGA). In addition, a prognostic model was set up on basis of recognized genes, the applicability and value of the model in BLCA were assessed; the generalizability of the model was studied using external verification using the Gene Expression Omnibus (GEO).

Materials and approaches
Achievement of BLCA expression profiles from TCGA datasets UCSC Xena provided RNA-seq data (level 3) and clinical data of BLCA samples. The Genome Tissue Expression (GTEx) [3] tool was adopted to collect expression levels of genes discussed in normal tissues. Then, the "sva" package of R software [4] was used to normalize RNA expression profiles and remove batch effects. GSE13507 [5] was also obtained from GEO, the platform for which was GPL6102. The GSE13507 cohort contained 165 primary BLCA samples, 58 normal-looking bladder mucosa samples surrounding cancer, 23 recurrent non-muscle invasive tumor tissues, and ten normal bladder mucosae.

Analysis of clinical phenotype of the Wnt signaling path
We downloaded the complete list of Wnt path genes from KEGG [6]. The list included several Wnt pathwayrelated genes (Additional file 1: Table S1). To identify the Wnt signaling pathway genes associated with BLCA outcome, we downloaded the relevant genes of the Wnt signaling pathway and performed univariate Cox regression analysis applying TCGA (P < 0.05). To characterize Wnt-related gene expression patterns in BLCA, RNA sequencing data in TCGA samples were analyzed. Subsequently, based on the median value of Wnt pathway activation in the sample, all patients were fallen into two groups (activated and inhibited), and KEGG was used for difference analysis to create a volcano map.

ssGSEA and WGCNA analysis for the analysis of Wnt-related biological functions
Based on TCGA BLAD samples, single-sample gene set enrichment analysis (ssGSEA) [7] analysis was employed for the evaluation of 6935 Gene Oncology (GO) [8] biological functions of each sample. Five Wnt pathways were recognized applying weight gene co-expression network analysis (WGCNA) [9] clinical phenotype, and then module assignment analysis was performed based on WGCNA. We obtained the correlation diagrams for the ssGSEA scores of GO functions and the Wnt pathway. We found epithelial-mesenchymal transition (EMT)related biological functions were the GO terms most associated with the Wnt pathway. Use WGCNA (version: 1.61) for co-expression network analysis. First, we select the soft threshold for building the network. The soft threshold gives the adjacency matrix a continuous value between 0 and 1, so that the established network meets the power law distribution and is near the actual state of the biological network. Second, we use the blockwise modules function to build a ladderless network, and then identify the gene co-expression modules that cluster genes with similar expression patterns by performing modulus assignment analysis. The definition of module is made by cutting the clustering tree into branches and assigning them different colors with the dynamic tree cutting algorithm. Then the Eigengene modulus (ME) of every modulus is calculated. ME represents the expression level of each module. Then, the association between the ME and the clinical characteristics of every module is calculated. At the end, the importance of genes in the EMT and immune cell infiltration. The most important finding was that these three genes were independent prognostic factors for BLCA. computation module represents the association between genes and samples. The soft threshold for network construction was 5. Finally, we conducted in-depth research on the correlation between EMT and Wnt pathway genes.

Consensus clustering and principal component analysis
Using the R software package "ConsensusClusterPlus" [10] from TCGA, based on the difference in expressing Wnt signaling pathway-associated genes in BLCA patients, we identified two subtypes (50 iterations, 80% resampling rate Pearson correlation, http:// www. bioco nduct or. org/). We used the cumulative distribution function and consensus matrix to calculate the appropriate number of subtypes. Then, we drew Kaplan-Meier overall survival curves and the classification of clinical stages based on these two subtypes.

Establishment and verification of risk scoring signature and line graph model
The LASSO-Cox algorithm [11] in the "glmnet" R software package was employed to identify the overall survival clinical outcome signature. We verified the prognostic prediction value of the signature conducted in TCGA-BLAD by using Sangerbox (http:// sange rbox. com/). Sangerbox is an R-based comprehensive analysis tool for bioinformatics.

Analysis of three essential genes in the immune microenvironment
As an algorithm, CIBERSORT [12] analyzes the ratio of cells in the body tissue gene expression matrix. As a genetic signature matrix, LM22 defines 22 immune cell subtypes that can be downloaded from the CIBERSORT website portal (https:// ciber sort. stanf ord. edu/). We used LM22 Matrix and cipher algorithm to analyze the ratio of CD8+ t cells. The result showed that the number of samples was P < 0.05. Subsequently, the "ESTIMATE" R software package [13] was employed to measure tumor unit (referring to the ratio of cancer cells in tumor), immune mark (referring to the infiltration of immune cells in tumor) and matrix mark (recording the presence of tumor cells).We evaluated the immune microenvironment panorama of the factors in signatures to calculate correlations among the Wnt pathway, the immune microenvironment, and immune responses.

Cell source and transfection
The

Cell proliferation assay
We use BeyoClick ™ EdU-488 Cell Proliferation Kit (Beyotime Biotechnology, China) to measure cell proliferation. We incubated the transfected cells according to the instructions and obtained images on a fluorescence microscope (Olympus Corporation, Japan), and used ImageJ software (National Institutes of Health, Bethesda, Maryland) Perform a cell count.

Wound healing assay
We seeded the cells in a 6-well plate. After 48 h of transfection, when the cells reached confluence, we used a pipette tip to cut the cells to form a mechanical wound.
Use an optical microscope to image at 0 and 24 h, and compare cell migration by calculating the size of the gap in each field of view.

Cell migration assay
We put the 8 μm Transwell chamber into a 24-well plate (Corning Costar, Corning, New York, USA), and preliminarily added 600 μl of medium containing 10%-FBS to each bottom well, and inoculated each upper chamber at the same time 8000 cells in suspension in 200 μl FBS-free medium. After incubation at 5% CO2 and 37 °C for 24 h, the upper layer film was removed in the cells, while the cells migrated to the lower side is fixed in 4% paraformaldehyde, and stained with 1.0% crystal violet. Use a microscope to obtain an image.

Immunoblotting assay
We extract total protein from the cells by RIPA lysis buffer with 1% PMSF, and then determine the protein concentration by the bicinchoninic acid assay (Beyotime Institute of Biotechnology Co., Ltd., Beijing, China) and the chemiluminescence system (Bio-Rad), California, USA).

Statistics analysis
We applied the R packages "glmnet" and "survival" to perform multivariate, Kaplan-Meier analysis and LASSO-Cox regression to evaluate risk signature, and use the R package "survivalROC" to perform Roc curve analysis [14]. R software (version4.0.3; https:// www.r-proje ct. org/) was used to perform all statistical analyses. It was considered that P < 0.05 was greatly significant. GraphPad Prism was also applied for data analysis.

The relationship between Wnt signaling pathway-related genes and phenotypic characteristics of BLCA
On basis of the univariate Cox regression model, we found that 34 factors in the Wnt pathway were greatly related to OS. A forest map of risk scores displayed that 32 genes were risk elements and two genes (MAPK10 and WNT7B) were protective factors (Fig. 1A). The t-test of completely randomized data displayed the variation in the expression levels of these 34 genes in the clinical stages of BLCA and World Health Organization stages II and III (Fig. 1B). We used KEGG to analyze the metabolic pathways and signal transduction pathways of tumors significantly enriched in the Wnt pathway activation and inhibition groups. We found that these tumor-related pathways were also highly activated (Fig. 1C). Therefore, we believe that the Wnt signal pathway promotes cancer.

The association between the Wnt signaling pathway and biological functions in BLCA
We used ssGSEA and WGCNA analysis to enrich the biological functions associated with the Wnt signaling pathway and tumors in TCGA. Three modules were obtained (dark turquoise, yellow-green and dark green) (Fig. 1D). Dark turquoise and yellow-green negatively correlated with Wnt pathways, and dark green positively correlated with Wnt pathways. In these three modules, we found that EMT-related biological GO terms were associated with the Wnt pathway (Table 1). EMT positively correlated with the activation of the Wnt signaling pathway. Therefore, we speculated that the high activation of Wnt signaling in BLCA would induce EMT in tumor cells, increasing tumor migration and invasion and affecting the outcome in terms of computational biology. Other biological functions were closely correlated with the activation or inhibition of Wnt signaling pathways, including angiogenesis, vascular endothelial cell migration, G protein-coupled receptor pathways, and the mitochondrial respiratory chain ( Table 1). According to these outcomes, the Wnt signaling pathway is closely related to EMT. We further explored the correlation between 34 Wnt pathway-related genes that affect BLCA outcome and EMT. We calculated the influence of these 34 genes on EMT (Fig. 1E) and found that PPP2CB, LEF1, CTNNB1 were positively correlated with EMT. RAC3 was significantly negatively related to EMT. Figure 1 illustrates univariate Cox analysis of 34 genes in the Wnt pathway that had significant regression coefficients and were associated with BLCA outcome and EMT. Based on these 34 genes, we classified the Wnt signaling pathway in patients with BLCA. The consensus clustering of Wnt-associated genes divides BLCA specimens into two clusters. Using the similarity with Wnt-related gene expression, we chose the value of k = 2 ( Fig. 2A).

The subgroup of Wnt signal expression in BLCA
On basis of the expression of these 34 genes, the BLCA samples from TCGA data set were fallen into two subgroups (C1 and C2). We compared the expression differences of these genes in these subgroups (Fig. 2B).
The high expression group was the C1 group, and the low expression group was the C2 group. The C1 group was the Wnt signaling pathway enrichment group, and the C2 group was the Wnt signaling pathway deficiency group. The C1 group had poor outcomes, and the C2 group had better outcomes. The risk ratio of the C1 group was 1.80 of the C2 group, and the confidence interval was 1.06-3.05 (P = 0.028) (Fig. 2C). According to the database, the proportions of clinical stages in the two groups differed. The proportion of stage II BLCA (WHO) in group C1 was 28%; the proportion of stage II BLCA (WHO) in group C2 was 42%; the proportion of stage IV BLCA (WHO) in group C1 was 36%, and the proportion of BLCA in group C2 Phase III (WHO) accounted for 24%. The difference between the two groups of Phase III (WHO) was small (Fig. 2D). These findings suggest that the BLCA outcome in Wnt signaling pathway enrichment group is worse. The outcome of stage IV patients was worse, the outcome of stage II patients was better, and no great variation was found between stage III patients and the Wnt signaling deficiency group.

The prognostic prediction model associated with the Wnt pathway in BLCA
Having found that the Wnt signaling pathway was related to BLCA outcome, we tried to establish a set of prognostic risk scores associated with the Wnt pathway based on BLCA. We found 34 Wnt pathway-related genes. After applying the LASSO-Cox algorithm to these 34 genes, we identified MAPK10, RAC3, and PPP2CB as independent prognostic genes (Fig. 3A, B). Considering the role of these three genes in BLCA outcome, a risk mark signature was set up for the integration of the effects of these genes using LASSO-Cox regression coefficients (Fig. 3C). Risk score = − 0.001 * MAPK10 + 0.002 * RAC3 + 0.001 * PPP2CB. Based on the median risk mark, BLCA patients were fallen into high-risk and low-risk groups in TCGA. The survival rate in the low-risk group was greatly higher than that of the high-risk group. According to the time-varying ROC curve, the region under the nomogram curve (AUC) used for the prediction of the 3-year survival rate was 0.62, while the 5-year survival rate predicted by the AUC was 0.6 ( Fig. 3D).
Universality is an essential indicator for evaluating prognostic models, and it is necessary to verify signatures externally using data from other sources. GEO data sets applying various chip platforms were selected for the confirmation of the prognostic predictive power of Wntassociated genes and the signature. We found that these gene signatures had prognostic significance in GSE13507 (Fig. 4A-C).

The association between the three genes and the immune microenvironment of BLCA
Various immune infiltration and anti-tumor effects in the immune tumor microenvironment are caused by immune cells and immune-related pathways. The immune infiltration status of the transcriptome of BLCA data was assessed by CIBERSORT algorithm. Eight immunerelated genes in 22 immune cells were used to assess immune infiltration in BLCA. We analyzed the expression of MAPK10, RAC3, and PPP2CB in BLCA, the risk score, and the infiltration of immune cells (Fig. 5). By comparing with the MAPK10 high expression group, the MAPK10 low expression group had more M2 macrophages. Compared with the PPP2CB high expression group, the PPP2CB low expression group had fewer M2 macrophages (Fig. 5A).
The ESTIMATE algorithm was applied to assess the association between immune infiltration status and tumor purity, and the calculation of tumor stage, survival status, tumor purity, ESTIMATE mark, stroma mark, and the immune mark of TCGA was made (Fig. 5B). The heat map displayed that when the expression of MAPK10 is low, the purity of the tumor is lower, the immune score is higher, the stromal cell score is higher, and the expression of eight immune-associated genes is higher. When the expression of PPP2CB is low, the expression of the eight immune-related genes is low.

Functional characteristics of PPP2CB
To further validate our bioinformatics analysis, we selected a risk gene (PPP2CB) and performed cell function experiments by silencing its expression using siRNA. To date, no studies have reported the function of this gene in BLCA cells. We performed functional verification to determine the possible effects of this gene on the BLCA cell lines T24 and UM-UC3. After downregulating the expression of PPP2CB, the proliferation and migration of BLCA cells were significantly inhibited (Figs. 6A, B; 7A). Using western blotting, it was observed that silencing PPP2CB led to downregulation of protein levels of N-cadherin and vimentin, suggesting that this gene promotes EMT in BLCA cells (Fig. 7B). These results provide us with a part of the experimental basis for verifying bioinformatics analysis and provide further insights for the further study of molecular functions and mechanisms.

Discussion
Several biological phenomena were controlled by the Wnt signal transduction cascade throughout all animals' growth and adult life, and aberrant Wnt signaling emphasizes various pathologies in humans [15]. In this study, we found that Wnt signaling is related to the occurrence, growth, metastasis, and recurrence of BLCA, and the Wnt signaling pathway is a prognostic risk factor for BLCA.
In BLCA, many upstream and downstream molecular members of Wnt signaling participate in various biological functions. For example, Wnt signaling regulate (See figure on next page.) Fig. 3 A,  angiogenesis [16][17][18][19], induces vascular endothelial cell migration [20,21], and participates in G-protein-coupled receptor signal transmission [22][23][24], regulation of mitochondrial respiratory chain [25][26][27][28], and EMT. Based on the WGCNA analysis, we found EMT-related GO terms strongly correlated with Wnt pathways. Therefore, we focused on the regulation mechanism of EMT and Wnt. EMT is a necessary process in morphogenesis. Epithelial cells lose cell-to-cell contact, polarity, and other epithelial cells' characteristics and obtain the unique characteristics of mesenchymal cells (including enhanced motility). The characteristics of EMT are the loss of E-cadherin on the plasma membrane, the increase of vimentin and fibronectin, and the growing EMT occurs in collectively migrating cells [29]. These findings suggest that EMT is closely related to tumor metastasis and recurrence.
EMT is regulated by signal transduction pathways including TGF-β, Notch, NF-κB, Wnt, and receptor tyrosine kinase [30]. The regulation of EMT by the Wnt signaling pathway has been determined including Wnt5a [31], FOXP3 [32], SERPINH1 [33], CUL4B [34], and SUFU [35], which can be used in BLCA, gastric cancer, and pancreatic cancer. Activating the Wnt signaling pathway induces EMT. The Wnt pathway also inhibits EMT, and EFEMP2 may inhibit BLCA progression through the Wnt/β-catenin pathway to prevent EMT [36]. The downregulation of E-cadherin causes the nuclear translocation of β-catenin and activation of canonical Wnt signaling [37]. As a marker gene of EMT, the gene slug also induces nuclear translocation of β-catenin [38]. Studies have shown a close association between the Wnt signaling pathway and EMT and two-way modulation between the two. Because BLCA often recurs and has high invasiveness and metastasis, we explored the correlation between the Wnt signaling pathway and EMT and the impact of the Wnt pathway on BLCA outcome. We found that patients with highly activated Wnt signaling pathways had a more unsatisfactory outcome, suggesting that the Wnt signaling pathway is cancer-promoting. We also obtained five important prognostic independent risk factors for BLCA based on the Wnt pathway (MAPK10, PPP2CB, LEF1, CTNNB1, and RAC3). These Wnt signaling pathway genes are closely associated with EMT. PPP2CB, LEF1, CTNNB1 positively correlated with EMT; RAC3 significantly negatively correlated with EMT; MAPK10 positively regulated and inhibited EMT.
Wnt signal transduction and its crosstalk with various immune cells have both positive and negative effects on tumor progression. For example, it can renew white blood cells or it can promote immune tolerance and limit the anti-tumor response.
The interaction between inflammatory cells and cancer cells has been researched well. Signal pathways such as TGF-β, NF-κB and Notch have been confirmed to affect the occurrence of EMT. Changes in the tumor microenvironment are also closely related to tumor proliferation and metastasis, affecting the development of cancer and the prognosis of patients. Pro-inflammatory factors and immune-related factors also play different roles in the process of tumor metastasis. For example, mitochondrial reactive oxygen species (ROS) can further affect tumor metastasis by regulating the inflammatory response [39][40][41]. Thus some natural antioxidants or anti-inflammatory drugs may also be used to inhibit tumor metastasis [39,42,43]. In colorectal cancer, including the spectrum from normal colorectal adenoma to cancer, infiltrating macrophages display high levels of Wnt2 and Wnt5a.  19:409 This finding suggests that the activation of paracrine Wnt by macrophages leads to cancer progression [44].
The current research showed that the Wnt/β-catenin pathway-related factors of the Wnt classic pathway andB cell transcription coactivator 9 (BCL9) were related to immune cell infiltration. The expression of BCL9/BCL9L negatively related to the infiltration of CD8+ T cells in triple negative breast cancer, and BCL9/BCL9L inhibited the infiltration of CD8+ T cells in the tumor microenvironment [45]. The Wnt pathway promotes tumor growth by inhibiting immune cell infiltration. Wnt ligands secreted by tumor cells stimulate tumor-related macrophages to polarize to the M2 subtype through the canonical Wnt signaling pathway, leading to tumor development and migration [46]. This finding was consistent with our results, showing that when PPP2CB is highly expressed, the infiltration of M2-type macrophages is more substantial. M2 type macrophages have lost the functions of tumor antigen peptide presentation and processing. When the number of M2 type macrophages increases, immunity weakens, and outcomes worsen, suggesting that PPP2CB is an oncogene. When the expression of MAPK10 is low, the infiltration of M2-type macrophages is more substantial, suggesting that MAPK10 is a tumor suppressor gene.
In summary, as the central transportation hub in tumor tissue, PPP2CB affects tumor cell migration by regulating EMT, and it also regulates the tumor immune microenvironment and affects tumor immunity by interacting with immune cells.

Conclusion
We set up a three-gene signature on basis of Wnt signaling-related genes, EMT, and immune infiltration in BLCA. These genes affect tumor proliferation as members of the Wnt pathway family and affect tumor metastasis by regulating EMT, even participating in immune infiltration. We also identified PPP2CB as a risk factor for the activity, proliferation, and migration of BLCA cells and established relevant immune infiltration assessments to predict their survival rate accurately. We expect that our findings will deepen the understanding of the molecular mechanisms taking part in the occurrence and development of BLCA and provide a unique method for discovering predictive biomarkers and selecting targeted therapies.
Abbreviations CDF: Cumulative distribution function; ME: Module Eigengene; GS: Gene significance; GTEx: Genome Tissue Expression; ROC: Receiver Operating Fig. 7 A Migration of T24 and UM-UC-3 after silencing the expressions of PPP2CB. The analysis of data was made by t-test, which was showed as the mean ± standard deviation. B The detection of PPP2CB, vimentin and N-cadherin proteins in T24 and UM-UC-3 was made after silencing the expressions of PPP2CB. *P < 0.05; **P < 0.01; ***P < 0.001