A Comprehensive Bioinformatics Analysis of Notch Pathways in Bladder Cancer

Simple Summary The Notch pathway is important in embryology and numerous tumor diseases. However, its role in bladder cancer (BCa) has not been deeply investigated thus far. Gene expression data are available for BCa, and bioinformatics analysis can provide insights into a possible role of the Notch pathway in BCa development and prognosis. Using this information can help in better understanding the origin of BCa, finding novel biomarkers for prediction of disease progression, and potentially opening new avenues to improved treatment. Our analysis identified the Notch receptors NOTCH2/3 and their ligand DLL4 as potential drivers of BCa by direct interaction with basic cell functions and indirect by modulating the immune response. Abstract Background: A hallmark of Notch signaling is its variable role in tumor biology, ranging from tumor-suppressive to oncogenic effects. Until now, the mechanisms and functions of Notch pathways in bladder cancer (BCa) are still unclear. Methods: We used publicly available data from the GTEx and TCGA-BLCA databases to explore the role of the canonical Notch pathways in BCa on the basis of the RNA expression levels of Notch receptors, ligands, and downstream genes. For statistical analyses of cancer and non-cancerous samples, we used R software packages and public databases/webservers. Results: We found differential expression between control and BCa samples for all Notch receptors (NOTCH1, 2, 3, 4), the delta-like Notch ligands (DLL1, 3, 4), and the typical downstream gene hairy and enhancer of split 1 (HES1). NOTCH2/3 and DLL4 can significantly differentiate non-cancerous samples from cancers and were broadly altered in subgroups. High expression levels of NOTCH2/3 receptors correlated with worse overall survival (OS) and shorter disease-free survival (DFS). However, at long-term (>8 years) follow-up, NOTCH2 expression was associated with a better OS and DFS. Furthermore, the cases with the high levels of DLL4 were associated with worse OS but improved DFS. Pathway network analysis revealed that NOTCH2/3 in particular correlated with cell cycle, epithelial–mesenchymal transition (EMT), numbers of lymphocyte subtypes, and modulation of the immune system. Conclusions: NOTCH2/3 and DLL4 are potential drivers of Notch signaling in BCa, indicating that Notch and associated pathways play an essential role in the progression and prognosis of BCa through directly modulating immune cells or through interaction with cell cycle and EMT.


Introduction
The Notch signaling pathway is a highly conserved ligand-receptor signaling pathway involved in various aspects of cancer biology that plays an essential role in cancer stemness, angiogenesis, epithelial-mesenchymal transition (EMT), tumor immunity, and drug resistance [1][2][3]. In recent years, Notch pathways have been highlighted in the progression and prognosis of bladder cancer (BCa) [4]. However, the literature results are quite different.

Acquisition of Data
Target gene expression data of normal bladder samples were downloaded from the GTEx (https://gtexportal.org/home/datasets; accessed on 20 November 2020) database [22,23]. The clinical datasets and target gene expression for non-cancerous samples and bladder cancer samples (TCGA-BLCA) were downloaded from TCGA (https://portal.gdc.cancer.gov/; accessed on 18 September 2020) [24,25]. HTseq-count files were used to tabulate the number of uniquely mapping reads for each gene. We used statistical software R (Version 4.0.2. https://www.r-project.org/; accessed on 20 June 2020) with R studio (Version 1.2.5042) [26] and "DESeq2" R package to normalize the data and performed batch correction of the gene expression between sample in order to make features comparable. All the data used in this study are open for research, and we used the data following the guidelines in the appropriate Data Use Certification Agreement.

Analysis of Expression Levels in Different Subgroups
Cell type-specific expression levels were not available in the TCGA and GTEX datasets; however, we could infer the relative abundance of tumor-infiltrating lymphocytes (TILS) by using gene set variation analysis (GSVA) based on gene expression profile retrieved from the tumor tissue-TCGA-BLCA, allowing lymphocyte subtype analysis.

Acquisition of Data
Target gene expression data of normal bladder samples were downloaded from the GTEx (https://gtexportal.org/home/datasets; accessed on 20 November 2020) database [22,23]. The clinical datasets and target gene expression for non-cancerous samples and bladder cancer samples (TCGA-BLCA) were downloaded from TCGA (https://portal.gdc.cancer.gov/; accessed on 18 September 2020) [24,25]. HTseq-count files were used to tabulate the number of uniquely mapping reads for each gene. We used statistical software R (Version 4.0.2. https://www.r-project.org/; accessed on 20 June 2020) with R studio (Version 1.2.5042) [26] and "DESeq2" R package to normalize the data and performed batch correction of the gene expression between sample in order to make features comparable. All the data used in this study are open for research, and we used the data following the guidelines in the appropriate Data Use Certification Agreement.

Analysis of Expression Levels in Different Subgroups
Cell type-specific expression levels were not available in the TCGA and GTEX datasets; however, we could infer the relative abundance of tumor-infiltrating lymphocytes (TILS) by using gene set variation analysis (GSVA) based on gene expression profile retrieved from the tumor tissue-TCGA-BLCA, allowing lymphocyte subtype analysis.

Correlation Analysis and Evaluation of the Diagnostic Value
For correlation analysis and correlation matrix visualization, we used the statistical software R studio with the packages "performanceAnalytics", "Hmisc", "pROC", and "corrplot". We analyzed the correlation between different target genes. In addition, according to the tumor and immune system interaction database (TISIDB; http://cis.hku.hk/TISIDB; accessed on 18 March 2021) [45,46], correlations between the Notch members (Notch receptors and Notch ligands) with lymphocytes and the relationship with immunomodulators in BCa were analyzed. Criteria of correlation co-efficients were defined as follows: 0.00-0.19 (very weak), 0.20-0.39 (weak), 0.40-0.59 (moderate), 0.60-0.79 (strong), and 0.80-1.0 (very strong). A p-value < 0.05 was considered statistically significant [43]. We used receiver operating characteristic (ROC) curve analysis to evaluate the diagnostic value of the genes in question. Additionally, in order to identify the critical factors of Notch signaling, we used orthogonal partial least squares discriminant analysis (OPLS-DA; multivariate statistical analysis by R packages "ropls"). For quality criterion, we chose R 2 X > 0.4 in the PCA model, and R 2 Y (goodness of fit parameter) and Q 2 (predictive ability parameter) > 0.5 in OPLS-DA [47,48].

Overall Survival (OS) and Disease-Free Survival (DFS) Analyses
Clinical data were downloaded from TCGA-BLCA and analyzed with R studio (packages: "ComplexHeatmap", "clusterProfiler", "survival", "survMisc", "survminer", "RCol-orBrewer"). Kaplan-Meier estimates of OS and DFS were performed to evaluate the prognostic value of target genes. We used Cox regression to find independent factors, landmark analysis, and time-dependent covariates to evaluate the prognostic value in case of crossing over between groups. We also used GSCALite (http://bioinfo.life.hust.edu. cn/web/GSCALite/, accessed on 20 March 2021) [49,50], a web server for gene set cancer analysis, to estimate whether DNA methylation of Notch member genes had effects on the overall survival. A p-value < 0.05 was considered statistically significant. Data on cancer-specific survival are not available in the TCGA-BLCA dataset.

Gene-Gene Interaction (PPI) Network Analysis and Gene Set Enrichment Analysis (GSEA)
We used GeneMANIA (http://genemania.org/co-expression genes; accessed on 12 November 2020) [51,52] to construct the gene-gene interaction network (GGI) and GSEA software (Version 4.0.3) [53] for the functional analysis of the BCa and non-cancerous samples. We included biological process (BP) of Gene Ontology (GO) term enrichment analysis, canonical pathway analysis based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, REACTOME database analysis, and immunologic signature gene set analysis. Analysis results of GSEA, satisfying the threshold of normalized enrichment score |NES| > 1, nominal p-value (NOM p-value) < 0.05, and a false discovery rate (FDR) q-value < 0.10 were considered statistically significant.

Body Maps of the Target Genes and Oncomine Analysis
The expression levels of Notch-related genes were body mapped for various cancers and normal tissues using Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn; accessed on 12 November 2020) [57,58]. Additionally, we used Oncomine (https://www.oncomine.org; accessed on 13 November 2020) [59,60] for metanalysis of the target genes in previous studies. Immunochemistry (IHC) stainings of the proteins coded by the target genes in BCa and control tissues were accessed from the Human Protein Atlas (HPA, http://www/proteinatlas.org/; accessed on 13 December 2020) [61,62]. A p-value < 0.05 was considered statistically significant.

Other Statistical Analyses
All statistical analyses were performed using R statistical software. Univariate analysis was performed using ANOVA, t-test, Wilcoxon test, Tukey's HSD, and permutation tests. The false discovery rate (FDR) was used to conceptualize the rate of errors in null hypothesis testing, and Bonferroni/familywise error rate (FWER) was used to adjust p-values. Correlation analysis was performed by Pearson test or Spearman analysis. Multivariate analyses were performed using the principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA) model. Receiver operating characteristic (ROC) curve analysis was performed to evaluate the diagnostic value. Kaplan-Meier estimates of OS and DFS, landmark analysis, and time-dependent covariates were used to evaluate the prognostic value of target genes. We used Cox regression for survival analysis to find the vital factors. A p-value < 0.05/adjusted p-value < 0.05 were considered statistically significant.

Results
Four hundred and six tumor samples (n = 406) and 19 non-cancerous samples were available from TCGA-BLCA, and 21 normal bladder samples were extracted from GTEx. The normalized data passed the quality control and kept the differences in expression levels ( Figure S1). We defined non-cancerous bladder tissues from TCGA-BLCA and normal bladder samples extracted from GTEx as controls (n = 40) (Table S1). The demographic and patient characteristics including tumor classifications are shown in Table 1.

Gene Expression of Notch Pathway-Related Genes in BCa and Controls
Compared with controls, the expression levels of NOTCH1, NOTCH2, NOTCH4, DLL1, and DLL4 were significantly downregulated in the BCa samples. On the contrary, NOTCH3, DLL3, and HES1 were significantly overexpressed. Interestingly, in controls, NOTCH4 and DLL4 in non-cancerous samples from TCGA-BLCA were significantly lower than in normal tissues from GTEx. No significant differences were found for JAG1, JAG2, and HEY1 ( Figure 2, Table 2). In the Supplementary Materials, we included a complete set of figures illustrating the expression of Notch-related genes with respect to tumor TNM stage ( Figure S2), lymph node status ( Figure S3), histologic type (papillary vs. nonpapillary; Figure S4), and molecular subtypes ( Figure S5).  , and significance levels are marked as * p-value < 0.05, *** p-value < 0.001, **** p-value < 0.0001. (B) Heatmap of expression levels (color-coded). High expression marked in pink; low expression marked in blue.

Gene Expression in Different Tumor Stages
In tumor stage I, only two samples were available. Therefore, we omitted stage I from the analysis between different tumor stages. The expression of NOTCH2 was lower in stage II (S2) than in stage III (S3) and stage IV (S4); JAG2 was higher in S2 than in S4 ( Table 2).

Gene Expression in Patients Stratified for Lymph Nodal Metastasis
Compared to controls, we found that NOTCH2 was significantly downregulated in N0, N1, and N2. NOTCH3 and HES1 were significantly upregulated in N0, N1, and N2. NOTCH4, DLL1, and DLL3 were significantly altered in a range from N0, N1, N2, to N3. Furthermore, no significant differences were found between N1, N2, and N3, and HEY1 was also not significantly changed between subgroups ( Table 2). The NOTCH gene expression was up-or downregulated without correlation to lymph nodal metastasis status (data not shown).

Gene Expression in Papillary (PT) and Non-Papillary Tumors (NPT)
Compared to control samples, NOTCH3, DLL3, and HES1 were significantly increased in PT and NPT, whereas NOTCH1, NOTCH2, NOTCH4, and DLL1 were significantly decreased in PT and NPT.

Gene Expression in Patients Stratified for Race and Gender
We found significant differential expression of NOTCH2 and DLL4 between Asians (ASI, n = 43), Caucasians (CAU, n = 323), and Black or African Americans (AFA, n = 23; Table 2). NOTCH1, JAG1, and JAG2 were significantly upregulated in AFA compared to CAU, while NOTCH4 was significantly upregulated in ASI compared to AFA. No significant differences were found in the expression levels of NOTCH3, DLL1, DLL3, and HEY1.
The cohort included n = 299 male and n = 107 female patients (gender ratio of 3:1). The gender ratio was 2:1 (n = 14 male; n = 5 female) in the control group (Table S1). Only JAG2 and DLL3 showed gender specific differential expression with higher expression levels in female patients ( Table 2).

Gene Expression in Molecular Subtypes
We investigated the gene expression in neuronal (NET), basal (BT), basal squamous (BST), luminal (LT), luminal infiltrated (LIT), and luminal papillary (LPT). Compared to non-cancerous tissue controls from the TCGA-BLCA dataset (C-TCGA), NOTCH1 was downregulated only in LT; NOTCH2 was downregulated in LT, LIT, and LPT; and NOTCH4 was downregulated in BST. In contrast, NOTCH3 was upregulated in LT, LIT, LPT, and BST (see Table 3 for a complete list of all investigated target genes).

Gene Expression in Patients Stratified for Tumor-Suppressor TP53 Mutation
According to the TP53 mutation status, we separated the BCa samples into the cases with mutation (TP53 M ) and without mutation (TP53 WT ). Compared to C-TCGA, patients with TP53 mutation showed downregulation of NOTCH1, NOTCH2, and NOTCH4, while NOTCH3 was upregulated. Patients without mutation showed downregulated NOTCH2 and upregulation of NOTCH4 only. Comparison of the TP53 mutation status in BCa tissue revealed in TP53 WT higher expression of NOTCH1 and NOTCH4 but lower expression of NOTCH2 (see Table 3 for a complete list of all investigated target genes).

Correlation and Diagnostic Value of Notch-Related Genes, and Correlation between Notch
Pathway, Lymphocyte Subtypes, and Immunomodulators

Potential Diagnostic Value of Notch-Related Genes
As shown in Figure 3B1, the majority of BCa cases (n = 406) could be discriminated from controls (n = 40) by the OPLS-DA scores plot. Moreover, this OPLS-DA model met the quality criterion R 2 X > 0.4. However, the fitting of the data was R 2 Y = 0.317 (p < 0.01) and the predictive power was Q2 = 0.298 (p < 0.01). Nevertheless, the variable importance in projection (VIP) in OPLS-DA identified eight potential critical genes in descending order: DLL3, NOTCH4, NOTCH2, HES1, DLL1, NOTCH1, DLL4, and NOTCH3 ( Figure 3B2).

Notch Pathway Correlated with Lymphocyte Subtypes and Immunomodulating Genes
Correlation analysis was performed between the Notch receptors and ligands and the infiltration level/relative abundance of lymphocytes in BCa. NOTCH2 showed the strongest, mostly positive, correlation with lymphocytes. In contrast, NOTCH3 was negatively correlated with the majority of lymphocytes ( Figure 4A). We focused on relationships with ρ > 0.4 and p-value < 0.05.

Dependency of Overall Survival (OS) on the Expression Levels of Target Genes
Kaplan-Meier estimates showed that high expression levels of NOTCH3, JAG1, DLL4, and HEY1 significantly correlated with poor OS. In contrast, high expression of HES1 was associated with prolonged survival (Figure 6). Interestingly, we found a crossing in the survival curves of NOTCH2 ( Figure 6B). Subsequently, the landmark analysis revealed that high expression of NOTCH2 seems to be a risk factor in the early phase of follow-up. However, after 100 months (8.2 years, n = 12; high expression: n = 10; low expression: n = 2), NOTCH2 resembled an advantageous factor, with high expression being associated with a better prognosis ( Figure S7A). Cox regression for overall survival (OS) analysis revealed NOTCH3, JAG1, DLL4, and HES1 as independent predictive variables in a multivariate test (Table 4).   Table S6 for threshold (cut-off value) used to define the high and low levels of the expression. y-axis = overall survival probability, x-axis = follow up in days, high expression (orange line), low expression (green line); the number at risk in high-and low-expression groups are listed in the tables below the survival curve; significance levels are marked as * p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001, **** p-value < 0.0001. Cox regression analysis; hazard ratio = HR; 95% confidence interval = CI95%; * = p-value < 0.05; ** = p-value < 0.01.

Dependency of Disease-Free Survival (DFS) on the Expression Levels of Target Genes
Kaplan-Meier estimates of DFS results showed that patients with high expression of NOTCH3, JAG1, JAG2, and HEY1 had a shorter DFS than those with low expression levels. By contrast, high expression of DLL4 and HES1 was associated with better DFS (Figure 7). Intriguingly, we also found a crossing in the disease-free survival curves of NOTCH2. Landmark analysis showed that DFS was greater with higher NOTCH2 expression at up to 90 months (7.4 years, n = 29; high level, n = 17; low level, n = 12); thereafter, low expression of NOTCH2 was associated with greater DFS ( Figure S7B). Cox regression for DFS analysis revealed NOTCH3, JAG2, and HES1 as independent variables in multivariate testing (Table 5).

Combinations of the Independent Factors Correlated with OS and DFS, While Methylation of Notch Factors Did Not
Notch ligand binding leads to proteolysis of the Notch receptor and activation of distinct target programs, which can induce different cell fates [5,64]. Consequently, the combination of certain Notch receptors and particular ligands may be most relevant, and therefore we analyzed the correlation of several combinations of Notch members with OS and DFS on the basis of COX regression results and linear fitting models (risk score = coefficient1 * gene expression 1 +...+ coefficient N * gene expression N) [65].
Additionally, low values of NOTCH3 + JAG2 (cut-off = 10.4, p-value = 0.018) and NOTCH3 + JAG2 + HES1 (cut-off = 4.1, p-value = 0.015) were associated with prolonged DFS ( Figure S8). However, we did not find a significant influence of methylation status of the target genes on OS or DFS (data not shown).  Table S6. y-axis = DFS probability, x-axis = follow up days; high expression is marked with an orange line, low expression is marked with a green line. The number at risk in high-and low-expression groups are listed in the tables below the survival curve; significance levels are marked as * p-value (p) < 0.05, ** p-value < 0.01.   Table S6. y-axis = DFS probability, x-axis = follow up days; high expression is marked with an orange line, low expression is marked with a green line. The number at risk in high-and low-expression groups are listed in the tables below the survival curve; significance levels are marked as * p-value < 0.05, ** p-value < 0.01. Cox regression analysis; hazard ratio = HR; 95% confidence interval = CI95%; * = p-value < 0.05; ** = p-value < 0.01.
The functional analysis of the network describes the roles of genes in Notch-associated pathways and their functions. In Figure 8A, we indicate four major functions as colored circle sections. Overall, we found 23 genes directly related to the Notch signaling pathway (blue): 11 of them (NOTCH1-4, DLL1, DLL4, JAG1, JAG2, APH1A, APH1B, NCSTN) were involved in Notch receptor processing (orange), 5 of them (DTX1, DLL3, DLL1, JAG1, JAG2) in Notch binding (red), and 4 of them (DTX1, LENG, NEURL1, HEY1) in the regulation of the Notch signaling pathway (pink). The functions overlapped in part ( Figure 8A).
Additionally, the results also revealed some significant and interesting interactions with other functions including stem cell development, cell fate determination, and immune system development. Further information regarding results of the gene network analysis of function is supplied in Table S8.

Notch Pathway Interactions with Other Well-Known Cancer-Related Pathways
Using GSEA, we evaluated the relationship between the target genes and other wellknown cancer-related pathways. Furthermore, we also calculated the percentage of BCa cases in which the Notch pathway was either activated or inhibited by Notch genes. Analysis revealed that the Notch pathway plays a critical role in cell cycle and EMT, as indicated by high positive (activating) or negative (inhibiting) values for NOTCH4, NOTCH2, NOTCH3, JAG1, DLL1, and DLL4 ( Figure 8F). The results also indicated that the target genes can have dual roles in the same pathway, e.g., NOTCH4 predominantly activates EMT (EMT_A = 35%) but can also inactivate EMT (EMT_I = 3%). Furthermore, there was a high variability in the functional role of ligands, Notch receptors, and downstream genes with respect to the cancer-related signaling pathways analyzed ( Figure 8F).

Meta-Analysis of Notch Factors in BCa
We performed an Oncomine meta-analysis based on the expression level between infiltrating bladder urothelial carcinoma vs. non-cancerous tissue (https://www.oncomine. org/, accessed on 13 December 2020) [59,60]. The meta-analysis revealed that NOTCH3 and HES1 were also significantly upregulated in BCa in previous studies [69,70] (Figure 9A,C), while the expression level of DLL4 was not significantly altered and has not been widely evaluated before [50][51][52] (Figure 9B). Interestingly, in contrast to the present study, NOTCH1, NOTCH2, and JAG1 were significantly overexpressed in BCa in Dyrskjot's study, while the other examined genes did not show any significant difference [69]. The detail information on the other target genes is listed in Table S13.

Body Maps of Notch Factors in Normal and BCa Patients
To visualize the distribution and to compare the expression of target genes in normal and BCa patients, we constructed the expression body maps of the genes from the Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn, accessed on 12 November 2020) [57,58]. Figure 9 depicts NOTCH3, DLL4, and HES1 as representative genes of the Notch pathway. The gene expression of NOTCH3 and HES genes was upregulated in BCa 2.47 times and 1.48 times, respectively ( Figure 9A1,A2,C1,C2), while DLL4 was downregulated in BCa 0.66 times (Figure 9B1,B2). Detailed information on the gene expression of all target genes is summarized in Figure S9.

IHC of Notch-Related Proteins in Normal Bladder Tissue and BCa Samples
Immunohistochemistry (IHC) protein expression data from the Human Protein Atlas (HPA; https://www.proteinatlas.org/, accessed on 13 December 2020) [62] confirmed the expression and alteration of NOTCH3 ( Figure 9A3,A4) and DLL4 proteins ( Figure 9B3,B4) in BCa, No protein expression data were available for HES1, NOTCH4, DLL1, and DLL3, whereas the others can be found in Figure S9.

NOTCH2/3 and DLL4 Are Potential Drivers of Notch Signaling in BCa
In recent years, great effort has been put into studying the mechanisms and prognostic and predictive values of the Notch pathway in BCa, including NOTCH1, NOTCH2, NOTCH3, NOTCH4, DLL1, DLL3, DLL4, JAG1, JAG2, HES1, and HES1. In reviewing the literature, prior studies have noted and focused on the importance of NOTCH1 and NOTCH2. However, very little of the literature has addressed the role of ligands in BCa [4]. Several previous studies have documented expression level alterations in BCa and have also described the Notch pathway's binary action, either oncogenic or suppressive (Table 6). Nevertheless, the biological function of the Notch factors is still not fully understood in BCa.
Oncogene (↑) NOTCH1 knockdown led to cancer cell growth and significantly inhibited growth and proliferation [72].

Oncogene (↑)
High rate of NOTCH2 copy number gain and over-expression in BCa, and the activation of NOTCH2 promotes metastasis, resulting in poor survival [5].

Suppression (↓)
Deletion of the intracellular domain of NOTCH3 leads to negative function in BCa, i.e., turning it into potential tumor-suppressive gene [8].

NOTCH4
Suppression * (↓) NOTCH4 was significantly downregulated in BCa, especially in the T1 stage of bladder tumor. However, in patients diagnosed with muscle-invasive bladder tumor (MIBC), high NOTCH4 expression correlated with poor survival and more vascular invasion [73].

JAG2
Suppression NA Oncogene (↑) JAG2 was overexpressed in BCa and was significantly associated with the metastasis and recurrence [76].
Oncogene NA

DLL3
Suppression NA Oncogene * (↑) DLL3 was significantly upregulated in small cell components correlated with worse clinical outcomes in small cell bladder cancer (SCBC) [77]. Despite apparent contradictions among precious studies (Table 6), in line with our findings, Greife et al. also reported that HES1 was overexpressed in BCa [6]. In our study, comparing BCa samples to non-cancerous tissue, the expression levels of notch receptors NOTCH1, NOTCH2, and NOTCH4 were significantly downregulated in BCa, while NOTCH3 was overexpressed. The ligands DLL1 and DLL4 were significantly downregulated in BCa, while DLL3 was overexpressed. The expression of the downstream gene HES1 was significantly upregulated in BCa. Our results confirmed previous research (in Table 6, marked with *).
The Oncomine meta-analysis results also support our present findings. In Lee's studies, all the target genes were detected in BCa (Table S13) [71]. Moreover, the research genes in our current study were also detected in normal bladder tissue through the RNA sequencing (RNA-seq) by the National Center for Biotechnology Information (NCBI, USA; https://www.ncbi.nlm.nih.gov/pmc/, accessed on 13 December 2020), (BioProject: PR-JEB4337, n = 95, human; Figure S10) [79,80]. HPA data were also in line with our findings documenting the detection of NOTCH1, NOTCH2, NOTCH3, JAG1, and JAG2 by IHC in normal bladder tissue and urothelial cancer (https://www.proteinatlas.org, accessed on 13 December 2020) [61,62]. In summary, previous reports as well as our present study support the hypothesis that all Notch receptors coexist in the BCa but potentially have quite different functions [5].
Suggestive evidence from previous studies revealed that the incidence, mortality, and prognosis of BCa depend on histopathology, race, gender, and molecular and morphological features [28,29,38,81].
In addition, compared to non-cancerous control samples, JAG1 was significantly upregulated in BST; JAG2 was significantly upregulated in NET, BST, and LPT; and HEY1 showed significant overexpression in NET, BST, LPT, and LIT. Furthermore, JAG1, JAG2, and HEY1 were differentially expressed in BST and LT.
Taken together, our findings support the hypothesis of biological subtypes of BCa [82,83], and we also found that NOTCH2, NOTCH3, DLL4, and HES1 expression was broadly altered in BCa subgroups (Tables S14 and S15).
Therefore, it could conceivably be hypothesized that NOTCH2, NOTCH3, and DLL4 are major factors in BCa and potential drivers of the Notch pathway. Additionally, the association of different Notch genes with molecular subtypes could promote our understanding of the biological role of the Notch pathway in BCa, thereby expanding therapeutic options, molecular diagnostics, and risk stratification [84][85][86]. Furthermore, these findings could also improve our understanding of the pathogenesis of BCa [83][84][85][86][87].

Differentially Expressed Notch-Related Genes Discriminate BCa and Relate to BCa Prognosis
Hu et al. documented that Notch4 receptor protein was differentially expressed in BCa compared to adjacent non-cancerous tissue and used Notch4 as one of 10 different biomarkers to discriminate MIBC from non-cancerous tissues. [73]. However, the diagnostic value of other factors is not recognized yet. In our study, ROC analyses revealed that NOTCH1-4, DLL1, DLL3, DLL4, and HES1 were able to distinguish the patients with BCa from normal individuals (BCa, n = 406; control, n = 40) with an accuracy (AUC) ranging from 0.678 (NOTCH3) to 0.831 (NOTCH4), sensitivities from 57.9% (NOTCH3) to 76.3% (DLL1), and specificities from 65% (DLL4) to 90% (NOTCH1). The OPLS-DA model confirmed the results implicated by the ROC analyses, also identifying NOTCH1-4, DLL1, DLL3, DLL4, and HES1 as potentially predominant factors, distinguishing BCa patients from normal controls.
Rampias et al. found that NOTCH1 alterations (mutations and/or copy number losses) were associated with poor OS, and lower expression of HEY1 and HES1 correlated with poor OS or DFS [8]. Zhang et al. pointed out that patients with high expression of NOTCH3 had a poor OS [9], and Li et al. described that patients with high JAG2 expression had a significantly poorer prognosis than those with weak expression [76]. On the other hand, high levels of JAG1 were associated with prolonged OS [74] (Table 6).
Furthermore, in the present study, the Kaplan-Meier curve of OS and DFS showed that high expression levels of NOTCH3, JAG1, and HEY1 were significantly associated with a worse prognosis (Figures 6 and 7). In contrast, high expression of HES1 was associated with better OS and better DFS. The Kaplan-Meier survival analyses also showed that high expression of DLL4 significantly correlated with a poor OS, but a better DFS. Interestingly, overexpression of NOTCH2 correlated with poor prognosis up to about eight years of follow-up, while after this time point, high expression levels of NOTCH2 were associated with a better prognosis in both OS and DFS. Cox regression analysis revealed that both NOTCH3 and HES1 were independent predictors of OS and DSF in BCa.
Further analysis of gene combinations suggested that independent of the single gene's effect on OS/DFS, patients with high expression of combinations involving NOTCH3 were associated with a worse prognosis of OS. This implies that NOTCH3 potentially is the most critical gene for tumorigenesis and progression in BCa.
In addition, our findings are in line with HPA analysis, which also proposed HES1 as a favorable prognostic biomarker with high expression correlates with a better prognosis in urothelial cancer (HPA; https://www.proteinatlas.org/, accessed on 13 December 2020; Table S16) [62].
Intriguingly, our OS and DFS analyses revealed an unexpected crossing of the Kaplan-Meier curves of NOTCH2-dependent survival after 8 years and 7.4 years, respectively. No other marker showed such a curve progression (Figures 6 and 7). The interpretation of these unexpected results is difficult. We can think of four scenarios that could explain such a crossing: (i) interaction of the Notch2 receptor with different ligands may lead to different downstream events, of which the causalities are at date not fully understood; (ii) Notch2 receptor activity could be low in those patients, despite high NOTCH2 gene expression; (iii) crosstalk of the Notch 2 pathway with still unrecognized other pathways; or (iv) activation of non-canonical Notch pathways could cause the different results. In any case, much more specific molecular biologic experimental work is required to elucidate those complex events.
In summary, these results support the notion that (i) Notch pathway-related gene and/or protein expression can serve as diagnostic biomarkers for BCa, with NOTCH4 in particular being the most promising biomarker candidate; (ii) the Notch pathway may play a critical role in the progression of BCa, while NOTCH3 is the most critical tumorigenesis gene, and HES1 the most protective gene [12,15]; (iii) NOTCH2 and DLL4 have a dual role and function in BCa and the interactions between different receptors, whereas ligands and downstream genes may explain the variable results [72,73]; (iv) all in all, according to the diagnostic, prognostic, and OPLS-DA analyses, NOTCH2, NOTCH3, DLL4, and HES1 are the most essential factors of Notch signaling in BCa, associated with clinical value (Table S17).

Notch Pathway Modulates the Development and Progression of BCa via Immune System
Notch signaling plays a crucial role in T cell development, and potentially plays a dual role in activating or inhibiting T cells [88]. In addition, Notch signaling could induce macrophage phenotype polarization to tumoricidal M1-type macrophage [89], and Notch2 signaling in CD8+ T cells is required for generating potent antitumor cytotoxic T cells [90]. Furthermore, Notch2 function seems to be essential for the proper development of the immune system under normal conditions [91], and gain-of-function mutations play a critical role in tumor immunity [92]. However, inactivation of the Notch2 pathway in CD8+ cytotoxic T cells has also been reported to impair tumor immune response [93]. Numerous studies have explored the role of Notch signaling in anti-tumor immune response in a variety of tumor entities, but thus far none have focused on BCa.
In the present study, GGI network showed that the Notch pathway potentially modulates the immune system in BCa. Additionally, GESA immune signature enrichment analysis also supported the fact that CD4, CD8, and dendritic cell (DC)-associated biological processes are linked to BCa. Using TISIDB, we made a systematic analysis that found the Notch pathway to be significantly associated with the immune system. Our findings confirmed previous studies in the literature, and the results showed that NOTCH2 was broadly and positively associated with the lymphocytes and immunomodulating genes, while NOTCH3 was generally and negatively associated with lymphocytes and immunomodulating genes. On the basis of the coefficient values, we found a moderate positive relationship between the expression level of NOTCH2 with Tcm CD8+ cells, Th2 cells, Mem B cells, Treg cells, and NKT cells. NOTCH2 seems to be the most crucial Notch factor and plays an essential role in modulating the progression of BCa through the immune system.
Moreover, the results also unveiled the Notch pathway modulating the immune system potentially via lymphocytes and immunomodulating genes [63]. Cancer genotypes determine tumor immunophenotypes and tumor escape mechanisms [63]. While gene expression studies cannot uncover molecular mechanisms, our results further strongly support the view that the Notch pathway is involved in regulation of the tumor immunity in BCa.

Notch Pathway Potentially Regulates or Crosstalks with Other Cancer Pathways
The Notch signaling pathway performs crucial roles in maintaining cell fate/proliferation, regulating the stem cell maintenance and differentiation and epithelial-mesenchymal transition (EMT) [94][95][96]. GGI network analysis revealed some exciting functions of Notch factors in BCa, including stem cell development and cell fate determination, indicating peculiar roles for DTX1, DLL3, DLL1, JAG1, and JAG2 in terms of Notch binding, while others are not involved. Furthermore, the canonical pathway and biological process enrichment from GESA also verified that "cell cycle", "DNA damage telomere stress-induced senescence", and "mitotic cell cycle checkpoint" are pathways associated with cell proliferation or cell fate and play a critical role in BCa. Abdou et al. described the fact that NOTCH1 over-expression in stromal cell potentially activates Oct-4 (one of the crucial regulators of cell self-renewal) and is associated with poor prognosis and liability for recurrence in BCa [96][97][98]. Furthermore, loss-of-function mutations in NOTCH1 and NOTCH2 promote an epithelial-mesenchymal transition (EMT) in BCa, stabilizing the mesenchymal phenotype [7]. However, the question remains as to how the other Notch factors contribute to BCa fate. Our analysis of the relationship between target genes and well-known cancer-related pathways showed that NOTCH4, NOTCH2, NOTCH3, DLL4, DLL1, and JAG1 can either stimulate or inhibit critical cancer-related pathways, e.g., EMT, DNA damage response, cell cycle, apoptosis pathways, TSC/mTOR, RAS/MAPK, and PI3K/AKT ( Figure 8F). Thereby, our findings corroborate the notion of the dual role of the Notch pathway depending on a delicate balance between its activating and inhibiting factors.
The tumor suppressor p53 plays an essential role in many malignancies, including promoting progression of BCa [99][100][101][102]. Dueñas et al. found that mutations in the TP53 gene and epigenetically mediated increased p53 protein expression favor the progression and recurrence of non-muscle invasive bladder cancer (NMIBC). Furthermore, mutations in TP53 and expression of splice variants were frequently observed in advanced stages and were associated with a worse prognosis [102]. However, the relationship between the TP53 and the Notch pathway was rarely reported in BCa until now. Consistent with the literature, canonical pathway enrichment based on KEGG and the Reactome Pathway Knowledgebase also showed significant enrichment in "P53 signaling pathway" and "TP53 regulates transcription of cell cycle genes". Interestingly, our results revealed that NOTCH1, NOTCH2, NOTCH4, DLL1, DLL4, and HES1 expression was significantly different between the patients with TP53 mutation and the cases without.
Taken together, these findings suggest that the Notch pathway plays a crucial role in regulating the progression of BCa through modulation of cell cycle and stemness, directly or indirectly interacting with well-known cancer-related pathways such as p53, TSC/mTOR, RAS/MAPK, and PI3K/AKT. Rather than a single pathway, carcinogenesis is often dependent on multiple pathways [14]. Hence, a better understanding of the related pathways and the molecular crosstalk is pivotal [103]. Consequently, instead of targeting a single pathway, targeting the crosstalk network could be a more meaningful strategy.

Limitations of the Study
Our study allowed us to integrate information of expression levels, mutations, and immune response coupled with pathological and medical information to uncover potential biomarkers and alterations of the Notch pathway in BCa. Our results promote our understanding of the complex influence of the Notch pathway on BCa and may help to improve clinical decision making by proposing certain members of the Notch pathway as potential prognostic biomarkers.
However, our study also has inherent limitations similar to the other bioinformatics approaches due to the restricted availability of cell type-specific data and the still very limited protein expression data. In addition, up until now, functional data as the activation status of proteins/enzymes and ligand-receptor interactions are not recorded in the TCGA and GTEx databases. Bioinformatics try to circumvent this problem by introducing network, gene-gene, and protein-protein interaction analyses on the basis of target gene expression overlaps and predicted interactions. However, specific and extensive experimental in vitro research is needed to confirm the hypothesized functions and to elucidate the complex interactions between the different members of the Notch pathway and the interactions with other cancer-related pathways. Therefore, our study cannot directly answer the urging question of potential therapeutical consequences.

Conclusions
In conclusion, our bioinformatics analysis provided further evidence for the involvement of the Notch pathway in carcinogenesis, progression, metastasis, and the modulation of the tumor immunity in BCa. The preponderance of the results indicates that NOTCH2, NOTCH3, and DLL4 are potential drivers of Notch signaling in BCa, playing the most crucial role in the progression and development of BCa via modulating the cell cycle and stemness, as well as direct or indirect interactions with well-known cancer-related pathways such as p53, TSC/mTOR, RAS/MAPK, and PI3K/AKT. Therefore, the Notch pathway is potentially a candidate target for diagnosis and therapies in BCa. Nevertheless, further work is required to perform the novel strategies, not only concern on the role of single gene/protein but also focusing on function of multigene/protein combinations in BCa. Our in silico study is limited to uncover associations and correlations but cannot provide mechanistic insights. Therefore, validation of the cancer biological functions in future in-depth experimental research is needed. Meanwhile, cytological studies are ongoing in our laboratory, seeking to confirm these results. Our findings can help in understanding the function of the core members of the Notch pathway in BCa and can provide further support for the hypothesis of the integral cancer pathway analysis approach.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cancers13123089/s1: File S1: List of full gene names and description of "GSE20727 CTRL vs. ROS INH and DNFB allergen-treated DC DN". Figure S1: Data before and after normalization. Figure S2: Expression of Notch-related genes in relation to tumor TNM stage. Figure S3: Lymph node status. Figure S4: Histologic type (papillary vs. non-papillary). Figure S5: Molecular subtypes. Figure S6: ROC analysis of JAG1, JAG2, and HEY1. Figure S7: Landmark analysis of Notch2 in OS and DFS. Figure S8: Analysis of combination in OS and DFS analysis. Figure S9: Body map construction and Oncomine meta-analysis and IHC of Notch pathway. Figure S10: The expression level of target genes in normal bladder tissue. Table S1: Demographics of controls in present study based on GTEx and TCGA-BLCA. Table S2: Diagnostic values of the components of Notch pathway. Table S3: Correlation between the components of Notch pathway and lymphocytes. Table S4: Correlation between the components of Notch pathway and immunoinhibitors. Table S5: Correlation between the components of Notch pathway and immunostimulators. Table S6: OS/DFS analysis of target genes and combinations in BCa. Table S7: Gene network analysis on interactions. Table S8: Gene network analysis results of function. Table S9: GSEA-canonical pathway enrichment based on KEGG. Table S10: GSEA canonical pathway analysis based on Reactome. Table S11: GSEA biological process of GO analysis. Table S12: GSEA immunologic signature gene set enrichment. Table S13: Meta-analysis of Notch BCa from Oncomine database. Table S14: Venn diagram data of subgroup analysis based on Notch receptors. Table S15: Venn diagram data of subgroup analysis based on Notch ligands and downstream genes. Table S16: Overall survival analysis of target genes in BCa based on HPA. Table S17: Venn diagram data of diagnostic, prognostic, and OPLS-DA analysis.

Data Availability Statement:
The data presented in this study are openly available in the Cancer Genome Atlas Program (TCGA), the Genotype-Tissue Expression (GTEx), and the Human Protein Atlas (HPA). Furthermore, the other web servers provide easy access to publicly available cancer OMICS data (TCGA, GTEx, et al.) and allow users to identify biomarkers or to perform in silico validation of potential genes, including GEPIA, GeneMANIA, Oncomine, GSCALite, UALCAN, and TISIDB for analyzing RNA sequencing expression and other data projects via a standard processing pipeline, and following the guidelines in the Data Use Certification Agreement.