{'fraction_non_alphanumeric': 0.03210674259847406, 'fraction_numerical': 0.014492753623188406, 'mean_word_length': 5.288682277318641, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SUMOylation Pattern Predicts Prognosis and Indicates Tumor Microenvironment Infiltration Characterization in Bladder Cancer

Background SUMOylation is an important component of post-translational protein modifications (PTMs), and bladder cancer (BCa) is the ninth most common cancer around the world. But the comprehensive role of SUMOylation in shaping tumor microenvironment (TME) and influencing tumor clinicopathological features and also the prognosis of patients remains unclear. Methods Using the data downloaded from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO), we comprehensively evaluated the SUMOylation patterns of 570 bladder cancer samples, and systematically correlated these SUMOylation patterns with TME immune cell infiltrating characteristics. The SUMO score was constructed to quantify SUMOylation patterns of individuals using principal component analysis (PCA) algorithms. Results Two distinct SUMOylation patterns and gene clusters were finally determined. Significant differences in the prognosis of patients were found among two different SUMOylation patterns and gene clusters, so were in the mRNA transcriptome and the landscape of TME immune cell infiltration. We also established a set of scoring system named SUMO score to quantify the SUMOylation pattern of individuals with BCa, which was discovered to be tightly connected with tumor clinicopathological characteristics and could predict the prognosis of patients with BCa. Moreover, SUMO score was a considerable predictive indicator for the survival outcome independent of tumor mutation burden (TMB) and low SUMO score was related to better response to immunotherapy using PD-1 blockade. We also found that there existed a significant relationship between sensitivity to commonly used chemotherapy drugs and SUMO score. Finally, a nomograph based on five features, namely, SUMO score, age, gender, T category, and M category was constructed to predict the survival probability of patients with BCa in 1, 3, and 5 years, respectively. Conclusions Our work demonstrated and overviewed the complicated regulation mechanisms of SUMOylation in bladder cancer, and better understanding and evaluating SUMOylation patterns could be helpful in guiding clinical therapeutic strategy and improving the prognosis of patients with BCa.


INTRODUCTION
Post-translational protein modification (PTM) is an important regulatory mechanism which regulates the interaction between protein and other molecules and changes the hydrophobicity, charge state, conformation and stability of proteins, thus affecting their functions in various biological processes (1). To date, more than 450 PTMs have been discovered, such as, acetylation, ubiquitination, phosphorylation, and methylation (2). SUMOylation, a specific kind of post-translational modification with the small ubiquitin-related modifier (SUMO), was identified in the mid-1990s (3) and found to play important roles in cell growth and migration, transcriptional regulation, immune response and tumorigenesis (4).
SUMOylation is a reversible and dynamic process mediated by several conserved proteins, namely, SUMO proteins (SUMO-1, SUMO-2, SUMO-3, and SUMO-4), SUMO E1 enzyme (SAE1-SAE2 heterodimer), SUMO E2 Enzyme (Ubc9), SUMO E3 Ligases, and SUMO Proteases (SENPs). A SUMOylation cycle consists of maturation, activation, conjugation, ligation, and de-modification as described below. Before being conjugated to target protein, SUMO protein should first be activated and matured by proteolytically cleaving the four amino acids at the C terminal to expose the glycine-glycine (GG) motif with the help of SUMO-specific proteases (SUPs), a family member of SENPs. Then the mature SUMO protein is activated by SUMO E1 enzyme and linked to the SAE2 subunit of SUMO E1 enzyme through the formation of a thioester bond in an ATP-dependent manner. Next the activated SUMO protein is transferred to Ubc9 by forming another thioester bond. Finally, with the assistance of a specific SUMO E3 ligase, Ubc9 catalyzes formation of an isopeptide bond between the Cterminal glycine of SUMO protein and a lysine (K) residue in the substrate. SUMOylation is a reversible process since the isopeptide can be hydrolyzed by SENPs, thus releasing SUMO protein to enter the next SUMOylation cycle (2,5,6).
Bladder cancer (BCa) is the ninth most common cancer around the world and the fourth most common cancer among men, with an estimated 550,000 new cases and approximately 200,000 related deaths in 2018 (7,8). It is also considered as the costliest highly prevalent cancer in the USA (9). Advanced age, male, tobacco smoking, exposure to aromatic amines and other chemicals are the main risk factors for the incidence of bladder cancer (8). Bladder cancer can be divided into non-muscle invasive bladder cancer (NMIBC) and muscle invasive bladder cancer (MIBC) according to the depth of tumor infiltration. NMIBC is confined to the mucosa (stage Ta, CIS) or submucosa (stage T1), for which transurethral resection of the bladder (TURB) and adjuvant intravesical chemotherapy instillations with bacillus Calmette-Guerin (BCG) or other drugs after TURB are conventional treatments (10). As for MIBC, which invades the muscle layer of bladder, is usually treated by radical cystectomy (RC) and some pre-or post-adjuvant therapy such as radiotherapy, chemotherapy and immunotherapy according to risk classification (11). The traditional chemotherapy follows a regimen of methotrexate, vinblastine, doxorubicin, and cisplatin (MVAC), which was proved to improve the prognosis of patients who accepted RC, however, this chemotherapy regimen has been replaced by gemcitabine plus cisplatin (GC) as its similar efficacy and better tolerability (12). Nowadays, due to deep studying and enhanced understanding of the tumor microenvironment (TME) and molecular profiling of bladder cancer, immunotherapies targeted at programmed cell death protein 1 (PD-1) and its ligand programmed cell death protein ligand 1 (PD−L1) have been developed over the past several years and expected to improve the survival outcome of patients with BCa in the future (12). However, only approximately 20% of patients have an objective response to immune checkpoint blockade (ICB), therefore novel therapeutic strategies for BCa need to be developed urgently.
SUMOylation has been proved to play critical roles in cell cycle progression by regulating the activity of cyclin-dependent kinases (CDKs), so aberrant SUMOylation can influence cell proliferation and lead to tumorigenesis (13). Many proteins involved in SUMOylation were found to be highly expressed in cancer tissues, namely, ovarian, colon and prostate cancer (14,15), suggesting SUMOylation was more activated in tumor and tightly linked to tumor growth. SUMOylation can also promote tumorigenesis via regulating oncogenes and tumor suppressors like c-Myc and FoxM1. Besides, like other PTMs, SUMOylation can occur in many tumor-associated proteins, thus producing new tumor antigens and providing novel therapeutic targets (16). Many previous studies have investigated the relationship between SUMOylation and BCa. SENP2 was reported to be downregulated in BCa and it could suppress migration and invasion of BCa cell in vitro through inhibiting the SUMOylation of transducin b-like 1 and transducin b-likerelated 1 (TBL1/TBLR1), thus inhibiting the nuclear translocation of b-catenin and the transcription of matrix metallopeptidase 13 (MMP13) (17). SENP2 was also found to regulate TGF-b signaling via the deSUMOylation of TGFb receptor I (TGF-bRI) in BCa, thus suppressing epithelialmesenchymal transition (EMT) in vivo and in vitro (18). Furthermore, Chen et al. also discovered the unique role of SUMOylation in inducing lymph angiogenesis and promoting tumor lymph node (LN) metastasis in BCa. The long noncoding RNA (lncRNA) LN-associated transcript 1 (ELNAT1) could induce Ubc9 overexpression, which would then mediate the SUMOylation of hnRNPA1 and promote packaging of ELNAT1 into the extracellular vesicles (EVs) secreted by bladder cancer cell. Finally, ELNAT1 was transported to human lymphatic endothelial cells (HLECs) by EVs and activated SOX18 transcription to trigger lymph angiogenesis (19). Given its specific roles in promoting tumorigenesis, metastasis and progression in BCa, SUMOylation could become a latent target for precise treatment in the future.
However, all the above studies focused on only one or two SUMOylation related proteins but no previous studies have investigated and summarized the comprehensive effect of SUMOylation. So, here in this article, we divided the patients with BCa into two SUMOylation patterns according to the expression levels of SUMOylation related genes and performed survival analyses, then we further explored the hallmark pathways and TME immune cell infiltration characteristics in distinct SUMOylation patterns and surprisingly found SUMOylation patterns were significantly connected with TME. Next, we discovered 1,934 SUMOylation cluster related differentially expressed genes (DEGs), classified the patients into two distinct genomic subgroups and explored the interaction between SUMOylation patterns and gene patterns. Moreover, we established a set of scoring system to quantify the SUMOylation pattern in individuals and explored the characteristics of SUMOylation in tumor clinicopathological characteristics, tumor somatic mutation, immunotherapy, and chemotherapy. Finally, we constructed a nomograph based on five features, namely, SUMO score, age, gender, T category, and M category to predict the survival probability of patients with BCa in 1, 3, and 5 years, respectively.

Enhanced SUMOylation Related Pathway and SUMOylation Patterns in Bladder Cancer
Having obtained 6 SUMOylation related gene sets and TCGA_BLCA dataset, we conducted gene sets enrichment analysis (GSEA) between bladder cancer tissues and normal tissues to investigate whether there existed enhanced or weakened SUMOylation functions in bladder cancer compared to the normal samples (23). Then we merged the transcriptional matrix of bladder cancer samples in TCGA_BLCA and GSE13507 dataset, and eliminated the batch effects by the "combat" algorithm using R package "sva". Following this, we extract the expression values of all genes enrolled in these 6 sumoylation related gene sets for all bladder cancer samples. Non-negative matrix factorization (NMF) algorithm was performed to identify the characteristics of SUMOylation patterns. Then we separately investigated the differential overall survival (Kaplan-Meier method survival curve and logrank test), enhanced sumoylation related functions [gene set variation analysis (GSVA)], and immune infiltrations [single sample gene set enrichment analysis (ssGSEA)] between SUMOylation patterns. Besides, the principal components analysis (PCA) was carried out to check the discrimination of the SUMOylation patterns.
Differentially Expressed Genes Between SUMOylation Patterns, DEGs Based Consensus Cluster, and Establishment of the SUMO Score R package "limma" was used to screen differentially expressed genes (DEGs) defined by the |log 2 fold change (FC)| >0 and adjusted p-value <0.001, then univariate Cox regression was conducted to check those DEGs with prognostic values. Following this, GO and KEGG gene set enrichment analysis were used to further explore the potential mechanisms or functions influenced by these DEGs with prognostic value. R package "ConsensusClusterPlus" was used to perform unsupervised consensus clustering. Similarly, prognosis differences, differentially enhanced or weakened functions and pathways were compared between DEGs based consensus clusters. Following this, the PCA algorithm was carried out to distinguish the molecular characteristics of these DEGs with prognostic value, and the SUMO score was calculated by the following formula: In the above formula. PC1 and PC2 represent two dimension expression patterns of these DEGs with prognostic value in all samples. Then the summary of these two scores can represent the SUMOylation levels of each sample to some extent.

Verification of the SUMO Score and External Dataset Validation in GSE69795 and GSE70691
Each sample has obtained a corresponding SUMO score according to its transcriptional matrix, then we firstly checked the best cut-off value of SUMO score in the TCGA_BLCA cohort to improve the survival predicting capability and obtain a best area under curve (AUC) by using R package "survminer". Then all samples from both TCGA_BLCA and GSE13507 cohorts were divided into high-/low-SUMO score groups. Then the survival analysis was performed to verify the efficacy of the SUMO score in the TCGA_BLCA cohort, GSE13507 cohort, and all samples.
Besides, we also performed survival analyses using external datasets GSE69795 and GSE70691 for further validation. Also, a gene set variation analysis (GSVA) was conducted to compare the differential active SUMOylation related functions between high SUMO score group and low SUMO score group. Immune cell infiltrations, immune related pathways, and immune related functions were quantified by the single sample gene set enrichment analysis (ssGSEA) algorithm. The correlation between SUMO score and ssGSEA results were conducted by the SPEARMAN correlation test. Having discovered the correlation between SUMO score and immune cell infiltrations, we wondered the correlation between SUMO score and tumor mutation burden (TMB). Thus, TMB of each sample was calculated by its somatic mutation profiles, and we further explored the correlation between SUMO score and TMB, then took these two factors to conduct joint prediction of the overall survival. Finally, the differences in the mutation atlas were compared between high SUMO score group and low SUMO score group by c 2 test.

Prediction of the Response to Chemotherapy/Immunotherapy and External Dataset Validation in IMvigor-210 Cohort
Having checked the effectiveness of the SUMO score, we wondered about the correlation between SUMO score and drug sensitivity to either immunotherapy or chemotherapy. Thus, we conducted three different algorithms to predict the response of each sample to the immunotherapy: TIDE (24), TCIA (22), and submap (25). All the predicted response to immunotherapy was compared by Wilcoxon or c 2 test. Besides, we calculated the SUMO score of each patient in IMvigor-210 cohort to externally validate the predicted response to immunotherapy. Responses to anti-PD-L1 immunotherapy in IMvigor-210 cohort were also compared between high-or low-SUMO score groups. As for the drug sensitivity to the chemotherapy, we used the "ProPhetic" package to separately predict the drug sensitivity to two commonly used chemotherapy regimens in bladder cancer (26), namely, GC scheme (gemcitabine plus cisplatin) and MVAC scheme (methotrexate, vinblastine, doxorubicin, and cisplatin). Furthermore, we compared the differential drug sensitivity to chemotherapy between high-/low-SUMO score groups.

Subgroup Analysis and Construction of SUMO Score Based Nomogram
Detailed clinicopathological characteristics of both TCGA_BLCA cohort and GSE 13507 cohort were collected and sorted. All samples were divided into several sub-groups containing age (>65 or ≤65), gender (male or female), grade (high grade or low grade), T stage (Ta or T1-2 or T3-4), N stage (N0 or N1-3 or Nx), M stage (M0 or M1 or Mx), and survival status (dead or alive), all of which were used to check the efficacy of SUMO score in different subgroups. Following this, we performed both univariate Cox regression and multivariate Cox regression to discover the potential independent prognostic factors between age, gender, grade, T stage, N stage, M stage, and SUMO scores. Finally, we constructed a nomogram including the independent prognostic factors, calculated the corresponding concordance index (C-index), and plotted the calibration curves and receiver operating characteristic (ROC) curve to verify the effectiveness of the nomogram.

Statistical Analysis
All the data cleaning, statistical analysis, and pictures plotting were conducted by R program version 4.1.1. A p-value <0.05 was defined as the statistically significant.

RESULTS
The Landscape of SUMOylation and the Establishment of SUMOylation Patterns in Bladder Cancer PTM, which includes acetylation, ubiquitination, phosphorylation, methylation, SUMOylation and so on ( Figure 1A), was found to play vital roles in many biological processes such as substrates transportation (27), signaling transduction (28), autophagy regulation (29), immune responses (30) and many others. As mentioned above, SUMOylation, as an important component of various PTMs, is a reversible process regulated by several highly conserved SUMO family proteins, namely, SUMO-1, SUMO-2, SUMO-3, SUMO-4, SUMO E1-activating enzyme consisting of two subunits named SUMO activating enzyme subunit 1 (SAE1) and SAE2, SUMO−conjugating enzyme E2 (SUMO E2 or Ubc9), SUMO E3 ligase and sentrin/SUMO-specific proteases (SENPs) (2,5,6), and mainly participates in nuclear function and intracellular signaling like chromatin organization, DNA damage response and repair, DNA methylation, DNA replication and immune response ( Figure 1B) (31). Here in Supplementary Table 1, we summarized genes involved in protein SUMOylation modification and SUMOylation related biological processes mentioned in Figure 1B. One dataset (GSE13507) with complete survival and clinical information was downloaded from the GEO database. Next, we used Combat R packages to eliminate the heterogeneity between the GSE13507 and the TCGA-BLCA cohort and enrolled them into a new metacohort. Before processing, we could easily distinguish the two datasets by principal component analysis (PCA) ( Figure 1C), while the two datasets merged together well after processing ( Figure 1D). The characteristics of these patients are shown in Table 1. As shown in Figure 2A, we found significant differences in the activation of SUMOylation related pathways between tumor and normal adjacent tumor tissue, and all these pathways were more activated in tumor than in normal adjacent tumor tissue. To further investigate the interaction between SUMOylation and tumor characteristics, we used Nonnegative Matrix Factorization (NMF) clustering to divide the patients in the meta-cohort into two SUMOylation patterns, termed as SUMO Cluster C1 and SUMO Cluster C2, based on choosing k = 2 as the optimal k value after calculating the cophenetic correlation coefficients of NMF ( Figure 2B and Supplementary Figure 1). Using PCA analysis, we could easily distinguish the two SUMOylation patterns ( Figure 2C), which revealed a remarkable difference in transcriptional profile of the SUMOylation related genes between the two different SUMOylation patterns. SUMO cluster C2 was characterized by higher expression level of SUMOylation related genes than SUMO cluster C1, indicating enhanced SUMO modification in SUMO cluster C2 ( Figure 2D). The clinical landscape of the two distinct SUMOylation patterns is shown in Figure 2E and there seemed to exist no significant differences in several clinical characteristics such as gender, age and TNM categories between the two SUMOylation patterns. However, a survival analysis for the two SUMO clusters revealed the particularly prominent survival advantage in patients in SUMO cluster C1 ( Figure 2G).

Hallmark Pathways and TME Cell Infiltration Characteristics in Distinct SUMOylation Patterns
In order to further explore the latent biological behavior differences behind the distinct SUMOylation patterns, we performed gene set variation analysis (GSVA) enrichment analysis based on KEGG datasets. As shown in Figure 2F, SUMO cluster C1 was dramatically enriched in pathways associated with metabolism, such as glycerophospholipid metabolism, retinol metabolism, ether lipid metabolism, and drug metabolism mediated by cytochrome P450, while SUMO cluster C2 showed enrichment in pathways related to cell proliferation and tumorigenesis, such as cell cycle, DNA mismatch repair, base excision repair, and DNA replication. Then we analyzed the TME immune cell infiltration and were surprised to find that SUMOylation patterns were closely linked to nearly all kinds of immune cells ( Figure 2H). SUMO cluster C1 was significantly enriched in few immune cells such as natural killer cell (NK cell), eosinophil, monocyte, and Type 17 T helper cell (Th17), while SUMO cluster C2 exhibited remarkable enrichment in most immune cells such as activated B cell, activated CD4 + T cell, activated CD8 + T cell, activated dendric cell (DC), myeloid derived suppressor cell (MDSC), natural killer T cell (NKT cell), Type 1 T helper cell (Th1), Type 2 T helper cell (Th2), and regulatory T cell (Treg). However, patients in SUMO cluster C2 did not show a corresponding survival advantage ( Figure 2G). MDSC is an important component of TME and can perform immune suppressive function via various mechanisms, namely, producing nitric oxide (NO), reactive oxygen species (ROS) and mostly antiinflammatory cytokines such as IL-10 and TGF-b, eliminating key nutrition factors which are vital to T cell proliferation, and preventing CD8 + T cell from entering the tumor region (32,33). Treg is another type of immunosuppressive cell in TME which can interfere the presentation of tumor associated antigens and attenuating cytotoxic T lymphocytes (CTL)-mediated cytotoxicity by influencing the release of cytolytic granules, which was demonstrated to contribute to immune evasion and tumorigenesis (34,35). As SUMO cluster C2 was remarkably enriched in MDSC and Treg, we speculated these two immunosuppressive cell types could play a crucial role in the poorer clinical outcome of patients with this SUMOylation pattern.

Generation of SUMOylation Gene Signatures and Functional Annotation
To further investigate the latent biological behavior of each SUMOylation pattern, we used the limma R package to discover SUMO cluster related DEGs. A total of 6,140 DEGs were identified between the two distinct SUMOylation patterns and finally we screened out 1,934 genes with prognostic value using univariate COX regression for further analysis. Then we performed GO and KEGG enrichment analyses for the DEGs by using the clusterProfiler R package. The results of the GO enrichment   analysis showed enrichment in biological process (BP), namely, neutrophil mediated immunity, cell cycle and antigen presentation, cellular component (CC), namely, chromosome and mitochondrial, and molecular function (MF), namely, ubiquitin-like protein ligase binding ( Figure 3A). The genes in the KEGG analysis also exhibited enrichment in pathways related to immunity, cancer, DNA repair and replication, which was consistent with previous results ( Figure 3B). The above results demonstrated again that SUMOylation was tightly linked with nuclear function, intracellular signaling and immune response. Then we performed unsupervised clustering analyses based on the 1,934 SUMO cluster related DEGs to dig out the potential regulation mechanism. We chose k = 2 as the optimal k value and successfully classified the patients into two distinct genomic subgroups using the unsupervised clustering algorithm ( Figure 3C and Supplementary Figures 2A-H). Then the cumulative distribution function (CDF) curve and scree plot were used to verify the rationality of the grouping (Supplementary Figures 2I, J). The track plot showed the details of grouping (Supplementary Figure  2K). We named these patterns as gene cluster A and gene cluster B, respectively. Similarly, gene cluster B was characterized by higher expression level of SUMOylation related genes than gene cluster A, indicating enhanced SUMO modification in gene cluster B ( Figure 3D). Thus, it was not surprising to find patients in gene cluster A had a better clinical outcome than those in gene cluster B ( Figure 3E). We also performed GSVA enrichment analysis using KEGG datasets between the two distinct gene clusters ( Figure 3F). To our surprise, gene cluster A was also remarkably enriched in pathways associated with metabolism, such as glycerophospholipid metabolism, fatty acid metabolism, retinol metabolism, and drug metabolism mediated by cytochrome P450, while gene cluster B exhibited enrichment in pathways linked to cell cycle and immune response to infection, which was approximately in accordance with the results in SUMOylation patterns. Interestingly, we also found peroxisome proliferator-activated receptors (PPAR) signaling was dramatically activated in gene cluster A. Previous studies have indicated that the activation of PPARg signaling can reduce the level of insulin-like growth factor-1 (IGF-1) in the blood, which was proven to promote tumorigenesis (36), thus interfering the growth of tumor in pancreas, colon, liver, and prostate (37). However, in bladder cancer, studies have shown that PPARg signaling could enhance the risk of developing bladder cancer (38,39), and promote the progression and metastasis of bladder cancer via providing a tumor-promoting microenvironment (37,40). But the results of our analyses showed the activation of PPARg was related to better prognosis in patients with BCa. Overall, the role of PPARg signaling in bladder cancer is still controversial and further studies are needed in the future.

Establishment of SUMO Score and Characteristics of SUMOylation in Tumor Mutation Burden and TME Cell Infiltration
The above analyses have revealed a landscape of the nonnegligible role SUMOylation played in BCa based on the patient population. Given the individual heterogeneity and complexity of SUMOylation, we next explored the latent SUMOylation pattern in individuals. Based on the expression levels of these prognostic DEGs, we constructed a set of scoring system named SUMO score to quantify the SUMOylation pattern of individual patients with BCa. Significant differences were observed in the genes participating in SUMOylation related to immune response, DNA replication, and DNA damage response between patients with high and low SUMO score ( Figure 4A). Patients with high SUMO score showed higher expression levels in the transcriptional profile of the SUMOylation related genes. The Sankey diagram was used to visualize the attribute changes of individual patients ( Figure 4B).
Not surprisingly, significant differences in SUMO score were also found between different gene cluster ( Figure 4C) and SUMO cluster ( Figure 4D). Gene cluster B exhibited higher SUMO score and so did SUMO cluster C2. Moreover, we also observed a significant survival advantage in patients with low SUMO score in all of the TCGA-BLCA cohort (Supplementary Figure 3A), GSE13507 cohort (Supplementary Figure 3B) and the meta cohort ( Figure 4E), which was consistent with patients in SUMO cluster C1 ( Figure 2G) and gene cluster A ( Figure 3E). We also performed survival analyses in external datasets GSE69795 (Supplementary Figure 4A) and GSE70691 (Supplementary Figure 4B) for further verification. Although lacking statistical significance, obvious survival differences were observed between patients in high and low SUMO score groups. Tumor mutation burden (TMB) has been considered a potential biomarker to predict the response to treatment with ICBs. High TMB could predict high ICB efficacy in bladder cancer (41,42). A multiple cancer research has also revealed patients with higher somatic TMB were associated with better overall survival (43). Consistent with previous studies, we also found patients with high TMB had better clinical outcome compared to those with low TMB ( Figure 4F). However, we did not find a significant correlation between TMB and SUMO score ( Figure 4G). Then, we further analyzed the distribution differences of somatic mutation between low and high SUMO score groups in the TCGA-BLCA cohort using maftools R package. As shown in Figures 5A, B and Table 2, in general, there were no obvious distribution differences of TMB in most mutated genes between low and high SUMO score groups, but for some popular genes in bladder cancer studies such as TP53, RB1, and FGFR3, the mutation frequency difference between high and low SUMO score groups was more than 10%. TP53 and RB1 were important tumor suppressive genes and studies have found they were regulated by PTMs such as SUMOylation, which could lead to abnormal proliferation and tumorigenesis (44). FGFR3 gene aberrations were common in bladder cancer and FGFR3 signaling was reported to promote epithelial-tomesenchymal transition (EMT) (45). But no previous studies have reported the relationship between SUMOylation and FGFR3, so further investigations were needed in this field. Nevertheless, we found the combination of TMB and SUMO score could excellently predict the clinical outcome of the patients with BCa ( Figure 4H). Patients with high TMB and low SUMO score had the greatest survival advantage while those with low TMB and high SUMO score had the worst. All these (H) Survival analyses for four groups grouped according to tumor mutation burden and SUMO score in the TCGA-BLCA cohort using Kaplan-Meier curves, namely, 62 cases in high TMB and high SUMO score group, 78 cases in high TMB and low SUMO score group, 138 cases in low TMB and high SUMO score group and 124 cases in low TMB and low SUMO score group. The high TMB and low SUMO score group showed significantly better overall survival than the other three groups (p <0.001, Log-rank test). (I) Correlations between SUMO score and the known immune cells in meta-cohort using Spearman analysis. Negative correlation was marked with blue and positive correlation with red. The number represented the correlation coefficient and the cross meant lack of statistical significance. results revealed SUMO score was a considerable predictive biomarker independent of TMB for the prognosis of patients with BCa. We also observed that a remarkable relationship existed between SUMO score and immune cell infiltration in the TME ( Figure 4I). SUMO score was positively correlated to the enrichment of immune cells such as activated B cell, activated CD4 + T cell, activated CD8 + T cell, DC, MDSC, NKT cell, Th1, Th2, and Treg, which was in accordance with the results in SUMOylation patterns ( Figure 2H) and indicated SUMO score could finely evaluate TME cell-infiltration characterization in individuals.

Characteristics of SUMOylation in Immunotherapy and Chemotherapy
Overall, patients with low SUMO score had better response to immunotherapies (Supplementary Figures 3C, D). We also found patients in the high SUMO group had higher Tumor Immune Dysfunction and Exclusion (TIDE) score ( Figure 5C), which could predict the response to ICB immunotherapy (24). Then we further explored whether the SUMO score had a predictive significance for the outcome of immunotherapies with different kinds of ICBs. As shown in Figure 5D, we evaluated the interaction between SUMO score and the response to anti-PD-1 and anti-CTLA-4 immunotherapy. We found low SUMO score was related to better response to anti-PD-1 immunotherapy. After Bonferroni correction, there still existed a remarkable correlation between low SUMO score and response to anti-PD-1 immunotherapy although lacked statistical significance. Then we divided the patients into four subgroups according to the use of anti-CTLA-4 and anti-PD-1 immunotherapies: CTLA-4 negative PD-1 negative ( Figure 5E), CTLA-4 negative PD-1 positive ( Figure 5F), CTLA-4 positive PD-1 negative ( Figure 5G) and CTLA-4 positive PD-1 positive ( Figure 5H). In CTLA-4 positive PD-1 negative and CTLA-4 negative PD-1 negative subgroups, low SUMO score was related to better immunotherapy response and prognosis, while in CTLA-4 positive PD-1 positive and CTLA-4 negative PD-1 positive subgroups the results were exactly opposite, which further proved that SUMOylation had a tighter relationship with immunotherapy targeted at PD-1 compared to other ICBs and patients with low SUMO score would benefit more from anti-PD-1 immunotherapy.
Anti-PD-L1 immunotherapy has been proven effective for patients with metastatic urothelial carcinoma in a multicenter, single-arm phase 2 trial using atezolizumab (IMvigor 210, NCT02108652) (20). Using the data acquired from IMvigor 210 cohort, we further verified the interaction between SUMO score and response to immunotherapies. Patients could be divided into three different immune phenotypes according to their response to immunotherapies: the immune-inflamed phenotype, the immuneexcluded phenotype and the immune-desert phenotype (46). The immune-inflamed phenotype was characterized by enrichment of various immune cells like CD4 + T cells, CD8 + T cells, and monocytic cells in the tumor parenchyma, as well as many proinflammatory and effector cytokines. Clinical responses to immunotherapies usually occurred in patients with this immune phenotype. The immune-excluded phenotype was also characterized by abundance in various kinds of immune cells, however, not in the tumor parenchyma, but in the stroma. So clinical responses to immunotherapies were uncommon in patients with this phenotype. The immune-desert phenotype was characterized by lack of T cells in both the tumor parenchyma and stroma, so it is not surprising that patients with this phenotype rarely responded to immunotherapies (46). Therefore, we further investigated the relationship between SUMO score and immune phenotypes. As shown in Figure 5I, the proportion of the immune-excluded phenotype was considerable between the high and low SUMO score groups, but the high SUMO score group was significantly enriched in the immune-desert phenotype and the low SUMO score group was dramatically enriched in the immune-inflamed phenotype. Similarly, we also found patients with the immune-desert phenotype had the highest SUMO score and patients with the immune-desert phenotype had the lowest ( Figure 5J). As immune phenotype was important to predict the response to immunotherapies, we next explored whether SUMO score was related to response to anti-PD-L1 immunotherapy. Not surprisingly, patients in the low SUMO score group showed significant therapeutic advantages and clinical response to immunotherapy ( Figures 5K, L). When combining complete response and partial response, as well as stable disease and progressive disease to create a binary outcome, the differences of SUMO score among patients with different main responses were more evidenced (Figures 5M, N). All these results revealed SUMO score could serve as a latent indicator for predicting the response to immunotherapies.
Next, we would like to investigate the relationship between SUMO score and response to chemotherapy. We screened out several commonly used drugs for intravesical chemotherapy in NMIBC such as doxorubicin, mitomycin C, and gemcitabine (10) and for adjuvant treatment after surgery in MIBC such as cisplatin, methotrexate, and vinblastine (11). We found that high SUMO score was related with low half maximal inhibitory concentration (IC50) in cisplatin, doxorubicin, mitomycin C, vinblastine and gemcitabine ( Figures 6A, B, D and Supplementary Figures 3E, F), which means higher sensitivity to chemotherapy, but methotrexate was just the reverse ( Figure 6C). We also found significant relationship between SUMO score and other drugs (Supplementary Figures 5G-N), some of which were newly developed and have not been put into clinical use to treat bladder cancer, thus we could screen out appropriate drugs according to the SUMOylation patterns of patients. In summary, the above results showed a unique role of SUMO score to predict the efficacy of immunotherapy and chemotherapy and could be used to guide clinical treatment.

Characteristics of Clinical Traits in SUMOylation-Related Phenotypes and the Construction of a Nomogram Using SUMO Score to Predict Prognosis in Clinical Scenarios
Then we lucubrated the interaction between SUMO score and clinical signatures and found SUMO score was significantly The violin diagram showed the differences of response index between high and low SUMO score groups among four subgroups. (E) If no immunotherapy was conducted, the high SUMO score resulted in a poor prognosis compared to low SUMO score (P <0.001, Wilcoxon test). (F) If only anti-PD1 immunotherapy was used, the low SUMO score resulted in a poor prognosis compared to high SUMO score (P = 0.034, Wilcoxon test). (G) If only anti-CTLA4 immunotherapy was used, the higher SUMO score group tended to get a poorer therapeutic efficacy compared to low SUMO score group (P <0.001, Wilcoxon test). (H) When anti-PD1 and anti-CTLA4 immunotherapy methods were simultaneously adopted, the high SUMO score group might get significantly better prognosis compared to low SUMO score group (P = 0.003, Wilcoxon test). related to age ( Figures 6E, F), gender (Figures 6G, H), grade (Supplementary Figures 5G, H), T category ( Figures 6I, J), M category (Supplementary Figures 5C, D) and final survival status (Supplementary Figures 5A, B). However, the distribution difference of SUMO score in N category did not show a statistical significance (Supplementary Figures 5E, F). In addition, we performed subgroup analyses and found SUMO score was a good predictor of survival in all subgroups divided by age (>65 and ≤65) ( Figures 6K, L), gender (female and male) ( Figures 6M, N), grade (high and low grade) ( Supplementary  Figures 5O, P) and T category (Ta, T1-2, T3-4) (Figures 6O-Q). But for subgroups divided by M and N categories, SUMO score was only significantly related to patients with M0 and N0 category (Supplementary Figures 5I, L).
Given the dramatic significance of SUMO score in predicting the prognosis of patients with BCa, we would like to establish a clinical predictive model to predict the survival probability of patients with BCa at 1, 3, and 5 years based on the clinical characteristics of the patients in the meta-cohort. First, we screened out five clinical features, namely, SUMO score, age, gender, T category and M category using univariate and multivariate COX regression model (Supplementary Figures  5Q, R). Then we constructed a nomograph based on these five clinical features which were easily accessible ( Figure 6R). The Cindexes for 1, 3, and 5 years were 0.823, 0.791, and 0.785 respectively, which indicated that the nomogram had a good predictive value (Figures 6S, T).

DISCUSSION
Nowadays, more and more studies have focused on the interaction between SUMOylation and tumorigenesis, metastasis, and progression (13), including bladder cancer (16)(17)(18)(19). As described above, many SUMOylation related proteins were overexpressed in tumor tissues, so drugs targeted at these aberrantly expressed components in SUMOylation process would be a latent treatment for cancer in the future. For example, ginkgolic acid, which could inhibit SAE1 subunit of SUMO E1-activating enzyme, was reported to suppress the growth of NOTCH1-activated breast epithelial cells and induce apoptosis (47). Treatment targeted at SUMOylation could not only directly impair the growth and survival of tumor, but also contribute to overcome the tumor heterogeneity and provide benefit to immunotherapy (16).
In the last decade, immunotherapies with ICBs targeted at PD-1 and PD-L1 were developed rapidly in BCa. Clinical trials using atezolizumab (20), nivolumab (48) and pembrolizumab (49) have all reported promising results, and immunotherapy has gradually become a new choice for those with high grade and metastatic disease (11). However, only approximately 20% of the patients could respond to ICB treatment because of the molecular heterogeneity of BCa, so it is important for us to find out biomarkers of genomic signatures to predict and select the patients who might have favorable outcomes during immunotherapies. Similarly, some biomarkers could also help in stratifying patients to predict the response to chemotherapy. For example, excision repair cross complementing 1 (ERCC1) has been considered as a prognostic biomarker to predict the sensitivity to cisplatin-based chemotherapy in patients with BCa (50). Therefore, it must be intriguing and urgent to establish a panel of validated biomarkers for response to chemotherapy and immunotherapy in the future.
Although previous studies have inquired the critical roles SUMOylation played in BCa, however, most of them paid more attention to limited SUMOylation related proteins and did not go deep into the comprehensive effect of SUMOylation in BCa and the interaction between TME cell infiltration and SUMOylation, which was necessary to guide more effective immunotherapy, chemotherapy and therapies targeted at SUMOylation.
Here in this article, we first investigated the differences in activated pathways related to SUMOylation between tumor and normal adjacent tumor tissue. Consistent with previous studies, we found all these pathways were more activated in tumor than in normal adjacent tumor tissue and the differences were significant. Then we combined the GEO dataset GSE13507 and TGCA-BLCA cohort into a new meta-cohort, and divided the patients with BCa in the meta-cohort into two SUMOylation patterns termed SUMO cluster C1 and C2 using NMF clustering according to the expression level of these SUMOylation related genes. Surprisingly, although we did not find significant differences in clinical characteristics, remarkable differences were discovered in the overall survival of patients and TME immune cell infiltration between the two SUMO clusters. It was interesting that SUMO cluster C2 was enriched in almost all kinds of immune cells but did not show a corresponding  (53). Studies also found SENP1 regulated the development and function of MDSC and the loss of SENP1 could increase SUMOylation of CD45, promote MDSC development and lead to tumorigenesis (54). Therefore, it is conceivable that the TME immune cell infiltration is tightly related to SUMOylation patterns. We also performed a GSVA enrichment analysis to dig out hallmark pathways and biological processes behind SUMOylation patterns, and we found SUMO cluster C2 showed enrichment in pathways related to cell proliferation and tumorigenesis, which was in accordance with the poor prognosis. Next, we used univariate COX regression to screen out 1,934 SUMOylation related DEGs with prognostic value and classified the patients into two distinct genomic subtypes named gene cluster A and B using unsupervised clustering based on the expression level of these genes. We found patients in gene cluster A had better survival outcome and the gene clusters were tightly connected with SUMOylation patterns. We also performed a GSVA enrichment analysis between the two distinct gene clusters gene and found cluster A was also remarkably enriched in pathways associated with metabolism, while gene cluster B exhibited enrichment in pathways linked to cell cycle and immune response to infection, which was approximately in accordance with the results in SUMOylation patterns. Interestingly, we found PPAR signaling was significantly activated in gene cluster A, but the previous studies have reported PPARg signaling could increase the incidence and progression of bladder cancer, which might lead to worse prognosis (37,38). So, it remains unclear what specific role PPARg signaling plays in the development and prognosis of bladder cancer.
Furthermore, given the individual heterogeneity and complexity of SUMOylation, it was necessary for us to explore the latent SUMOylation pattern in individuals. Therefore, we constructed a set of scoring system named SUMO score to quantify the SUMOylation pattern of individual patients with bladder cancer. We found significant differences in SUMO score between different gene clusters and SUMO clusters. SUMO cluster C2 and gene cluster B showed higher SUMO score. Similarly, patients with low SUMO score exhibited a significant survival advantage. As TMB has been demonstrated as an effective biomarker to predict the response to treatment with ICBs in BCa (43), we further explored the relationship between TMB and SUMO score. However, we did not find significant correlation between TMB and SUMO score and there were also no obvious distribution differences of TMB in most mutated genes between low and high SUMO score groups. But we found remarkable survival differences when combining SUMO score and TMB to divide the patients into four distinct subgroups, which indicated that SUMO score might be a potent predictive biomarker independent of TMB for the prognosis of patients with BCa. We also observed there existed a remarkable relationship between SUMO score and immune cell infiltration in the TME, which was consistent with the results in SUMOylation patterns and implied SUMO score could finely evaluate TME cellinfiltration characterization in individuals.
Many patients have benefited from immunotherapy using ICBs such as PD-1, PD-L1 and CTLA-4 blockade, but many more patients did not see pronounced clinical response to immunotherapeutic intervention (55). In patients with BCa, PD-1/PD-L1 blockade has demonstrated significant benefit in patients with unresectable and metastatic BC in the second-line setting, either as monotherapy or in combination with chemotherapy or CTLA-4 checkpoint inhibition (11). The results of the phase II trial using the PD-1 inhibitor pembrolizumab reported a complete pathological remission (pT0) in 42% and pathological response (<pT2) in 54% of patients (49), whereas another single-arm phase II trial with atezolizumab showed a pathologic complete response rate of 31% (20). These results suggested that the response rate was still low and it was important to screen out patients who were appropriate for immunotherapy. Our results found lower SUMO score was connected with better response to immunotherapy using PD-1 blockade. Next, using the data acquired from IMvigor 210 (NCT02108652) cohort, we further verified the interaction between (E) The proportion of patients with different age in low or high SUMO score groups. Age ≤65/age >65: 45%/55% in the low SUMO score groups and 34%/66% in the high SUMO score groups. (F) Differences in SUMO score between age >65 and age ≤65 subgroups in meta-cohort (P = 0.024, Wilcoxon test). (G) The proportion of patients with different gender in low or high SUMO score groups. Female/male: 20%/80% in the low SUMO score groups and 28%/72% in the high SUMO score groups. (H) Differences in SUMO score between female and male subgroups in meta-cohort (P = 0.011, Wilcoxon test). (F) Differences in SUMO score between age >65 and age ≤65 subgroups in meta-cohort (P = 0.024, Wilcoxon test). (I) The proportion of patients with different T category in low or high SUMO score groups. T1-2/T3-4/Ta: 52%/40%/8% in the low SUMO score groups and 34%/64%/1% in the high SUMO score groups. (J) Differences in SUMO score among T1-2, T3-4 and Ta subgroups in meta-cohort (Kruskal-Wallis H test). (K-S) Survival analyses for high and low SUMO score groups in subgroups of metacohort using Kaplan-Meier curves. (R) Nomogram for patients predicting survival outcome for 1, 3, and 5 years. Gender, M category, SUMO score, T category and age were marked as "points". Total points by adding the five points could predict survival probability. One patient who was male and had tumor with T2M0 category was randomly selected for analysis. After adding the score of each item, a total score of 297 and the corresponding survival probability for 1, 3, and 5 years respectively were obtained. (S) Calibration curves of nomograms in terms of the agreement between predicted and observed 1, 3, and 5 years of outcomes in meta-cohort. (T) Receiver operating characteristic (ROC) curve. The area under curve (AUC) for 1, 3, and 5 years was 0.823, 0.791, and 0.785 respectively, which showed a favorable ability of discrimination.
SUMO score and response to immunotherapy. Not surprisingly, we also found low SUMO score was connected better response to anti-PD-L1 treatment. Besides, we discovered the immune phenotype was also tightly linked to SUMO score and low SUMO score group was dramatically enriched in the immune-inflamed phenotype, which indicated a better response to immunotherapy. Therefore, we showed that SUMOylation patterns played a nonnegligible role in distinguishing different TME and could be an effective biomarker to select patients who were appropriate for immunotherapy.
Moreover, we also explored the relationship between SUMO score and response to chemotherapy. Mitomycin C (MMC), interferon (IFN), adriamycin, epirubicin, and gemcitabine (GEM) were commonly used drugs for intravesical chemotherapy after TURB apart from Bacille Calmette-Gueŕin (BCG) to prevent recurrence in patients with NMIBC (56). The traditional chemotherapy for MIBC followed a regimen of methotrexate, vinblastine, doxorubicin, and cisplatin (MVAC), and was gradually replaced by a new plan including gemcitabine plus cisplatin (GC) (12). Therefore, we evaluated the relationship between SUMO score and sensitivity to commonly used chemotherapy drugs and found IC50 for all these drugs exhibited a significant difference between high and low SUMO score group, which indicated a unique role of SUMO score to predict the efficacy of chemotherapy and guide clinical treatment.
Finally, we lucubrated the interaction between SUMO score and clinical signatures and found SUMO score was significantly related to age, gender, grade, T category, M category, and final survival status. We also performed subgroup survival analyses and found SUMO score was a good predictor of survival. Then we constructed a nomograph based on five features, namely, SUMO score, age, gender, T category, and M category to predict the survival probability of patients with BCa in 1, 3, and 5 years, respectively. The calibration curve and AUC of ROC curve showed a favorable ability of discrimination.
In general, our study provided a comprehensive insight into the interaction between SUMOylation, tumor clinicopathological characteristics, TMB, TME immune cell infiltration, and response to chemotherapy and immunotherapy. We demonstrated that different SUMOylation patterns could contribute to distinguishing the landscape of TME immune cell infiltration and clinical characteristics among patients, which was further verified using SUMO score within individuals. Moreover, SUMO score could also function as a predictive indicator for the survival of patients independent of TMB. Finally, we also evaluated the ability of SUMO score to predict the response to immunotherapy using ICBs and chemotherapy, which might help in improving therapeutic strategies, screening patients eligible for immunotherapy or chemotherapy and guiding individual precision therapy in the future.
However, we also realize that there still exist several shortcomings and limitations in our study. First, as many new SUMOylation related proteins were discovered recently, the SUMOylation related genes which we collected and used for analyses might not be comprehensive enough, which would bring some bias into our studies. Second, the current omics data only provide the level of mRNA but the SUMOylation process relies on proteins, which will bring in some inaccuracies.
Third, the number of clinical samples is limited and our study is lack of verification from other clinical data sets apart from the public data which will be helpful to further confirm our conclusions, and whether SUMOylation has a similar role in other types of cancer has not been verified. Therefore, we are prepared to collect some clinical samples to further verify our conclusions, and assess the role of SUMOylation in other urinary system tumors. Finally, the specific mechanisms in the interaction between SUMOylation patterns and TMB immune cell infiltration remain unclear and need further studies.
In conclusion, our work demonstrated and interpreted the complicated regulation mechanisms of SUMOylation in bladder cancer. The differences in SUMOylation patterns in population or individuals could significantly influence the heterogeneity in tumor clinicopathological features and TME, thus influencing the response to immunotherapy and chemotherapy. Therefore, better understanding and evaluating SUMOylation patterns could be beneficial in selecting appropriate patients, guiding clinical therapeutic strategy and improving the prognosis of patients with BCa.

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.

AUTHOR CONTRIBUTIONS
Q-DX and J-XS analyzed the data, wrote the manuscript and drew the figures. S-GW, ZL, and JH designed the study. YX, JX, C-QL, J-ZX, YA, and M-YX contributed to the critical revision of the manuscript. Q-DX, J-XS, JH, and S-GW contributed equally to this work. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication. analyses for low (13 cases) and high (25 cases) SUMO score patient groups in GSE69795 cohort using Kaplan-Meier curves (P = 0.183, Log-rank test). (B) Survival analyses for low (12 cases) and high (37 cases) SUMO score patient groups in GSE70691 cohort using Kaplan-Meier curves (P = 0.074, Log-rank test).
Supplementary Figure 5 | Survival analyses for patients with high and low SUMO score in several subgroups. (A) The proportion of patients with different final survival status in low or high SUMO score groups. Alive/dead: 68%/32% in the low SUMO score groups and 45%/55% in the high SUMO score groups. (B) Differences in SUMO score between alive and dead subgroups in meta-cohort (P < 0.001, Wilcoxon test). (C) Differences in SUMO score among M0, M1 and MX groups in meta-cohort (Kruskal-Wallis H test). (D) The proportion of patients with M0, M1 and MX stage tumor in low or high SUMO score groups. € Differences in SUMO score among N0, N1-N3 and NX groups in meta-cohort (Kruskal-Wallis H test). (F) The proportion of patients with N0, N1-N3 and NX stage tumor in low or high SUMO score groups. (G) The proportion of patients with different tumor grade in low or high SUMO score groups. High grade/low grade: 62%/38% in the low SUMO score groups and 91%/9% in the high SUMO score groups. (H) Differences in SUMO score between high and low tumor grade groups in meta-cohort (P < 0.001, Wilcoxon test). (I-P) Survival analyses for high and low SUMO score groups in subgroups of meta-cohort using Kaplan-Meier curves.