Identification and Validation of Molecular Subtype and Prognostic Signature for Bladder Cancer Based on Neutrophil Extracellular Traps

Abstract Neutrophil extracellular traps (NETs) could promote tumor growth and distant metastases. Molecular subtypes of bladder cancer were identified with consensus cluster analysis. A NETs-related prognostic signature was constructed with LASSO cox regression analysis. As a result, we identified three subtypes of bladder cancer, which had a distinct difference in prognosis, immune microenvironment, TIDE score, mRNAsi score and IC50 score. We also developed a prognostic signature based on 5 NETs-related genes, which had a good performance in clinical outcome prediction of bladder cancer. These results may provide more data about the vital role of NETs in bladder cancer.


Introduction
Bladder cancer (BLCA) is the most common and deadly cancer of urinary malignancies. There is an estimated of 573278 new cases of bladder cancer and 212,536 cases of bladder cancer-related death, ranking 3.0% of all new cancer cases and 2.1% of cancer-related deaths in 2020 globally [1]. Painless hematuria is the typical symptom of BLCA. BLCA has a high recurrence rate, so that patients still need to undergo cystoscopy frequently after surgery, resulting in serious economic and psychological burden [2][3][4]. The prognosis of BLCA patients is poor once metastasis and disease recurrence occur [5,6]. Recent studies suggested that immunotherapy based on checkpoint blockade was one of the most promising therapeutic strategies for BLCA patients with metastatic or advanced disease [7][8][9][10]. However, limited biomarkers could be used to predict immunotherapy benefits clinically for patients received immunotherapy. Thus, mining and identification of biomarkers predicting the prognosis and therapy response for bladder cancer is particularly important.
Neutrophil extracellular traps (NETs) are filamentous structures composed of extracellular DNA, granular proteins, and histones [11]. As key players in innate immune response and pathological progress, NETs exert a vital role in many diseases, including myocardial infarction [12], periodontitis [13], atherosclerosis [14], rheumatoid arthritis [15]. Tumor cells could recruit neutrophils and release NETs to tumor microenvironment [16]. As far as we know, the role of NETs in cancer is complicated and far from elucidated. On the one hand, NETs could promote tumor relapse by awaking dormant cancer cells [17]. By secreting matrix metalloproteinases, NETs could facilitate tumor growth and distant metastases [18]. On the other hand, NET components, including myeloperoxidase and histones, could directly kill tumor cells [18,19]. Recent study revealed that NET-related signature could predict the clinical outcomes and immunotherapy response in certain type of cancer [20]. Jing et al identified a prognostic signature for ccRCC using NETs-related genes [21]. Another study constructed a NETs-associated lncRNA signature predicting the clinical outcomes in LUAD [22]. These evidences indicated that clarification of the crosstalk between BLCA and NETs may help to devise potential approaches for the prognosis and treatment of BLCA.
In our study, we identified three NETs-related molecular subtypes and a prognostic signature for bladder cancer by using consensus cluster analysis and least absolute shrinkage and selection operator (LASSO) analysis, respectively. Moreover, we also explored the potential lncRNA-miRNA-mRNA regulatory axis of hub gene. Our study could provide more evidences about the role of NETs in BLCA.

Data source
NETs-related genes were defined as those genes correlated with neutrophils and NETosis. From the study of Şenbabao glu et al [23], we obtained the gene set of neutrophils. The NETosis-related gene set was a summary of the research progress of NETs in immunity and various diseases [24]. To summarize, we converged a total of 69 genes as the NETs-related genes for further study (Supplementary table 1). We could also see similar study from Yi Zhang et al [25]. From TCGA database, we collected the mRNA level profile, single nucleotide variants (SNV) and copy number variation (CNV) data of bladder cancer cases and matched clinical characters. Moreover, GSE13507 and GSE32894 dataset of bladder cancer were employed as test dataset for prognostic signature validation. Using student T-test, we identified the differentially expressed NETs-related genes with the "limma" package (https://bioconductor.org/ packages/release/bioc/html/limma.html) [26] in R.

Construction and validation of prognostic signature
To identify NETs-related genes with significant prognosis in bladder cancer, univariate cox regression analysis was conducted. This was followed by LASSO analysis to determine the coefficient of candidate genes and construct a NETs-related prognostic signature with 10-fold cross-validation. In this case, the risk score could be calculated with Ri(CoefiÁExpi) (Coef is the coefficient and Exp is the gene expression). We also employed ROC curve to compare prognostic model with other pathological features in prognostic prediction in BLCA using "survivalROC" package (https://cran.r-project.org/web/packages/ survivalROC/index.html). The C-index was calculated with "survival" package. These analyses were performed in TCGA dataset (training cohort) and GSE13507/GSE32894 dataset (test cohort).

Go and KEGG pathway analysis
As shown in Figure 2(A,B), these differentially expressed NETs-related genes were enriched in glycosaminoglycan binding, cytokine binding, chemokine binding, leukocyte migration, response to molecule of bacterial origin, positive regulation of interleukin-8 production, and neutrophil mediated immunity. Moreover, KEGG pathways analysis suggested the involvement of these differentially expressed NETsrelated genes in neutrophil extracellular traps formation, Toll-like receptor signaling pathway, PIK-ATK signaling pathway, PD-L1 expression and PD-L1 checkpoint pathway in cancer, and IL-17 signaling pathway (Figure 2(C,D)).

Consensus clustering identified three subtypes of bladder cancer
Based on the consensus CDF values ( Figure  3(A)) and delta area (Figure 3(B)), k ¼ 3 was suggested as the optimum number and three subtypes of bladder cancer were identified (Figure 3(C)). The expression pattern of NETsrelated genes was presented in Figure 3(D). Further survival analysis indicated that cluster 3 was correlated with favorable clinical outcome and cluster 1 had a poor overall survival in BLCA (Figure 3(E), p ¼ 0.0026). Previous study had suggested the crosstalk between NETs and tumor microenvironment [16]. In our study, the CIBERSORT score of Macrophage M0, Macrophage M1, resting Myeloid dendritic cell, activated Myeloid dendritic cell, Neutrophil, naïve B cell, plasma B cell, naïve CD4 þ T cell, resting CD4þ memory T cell, and Tregs in these three clusters was significantly different ( Figure  4(A)). As presented in Figure 4(B,C), the level of immune checkpoints and human leukocyte antigen (HLA) genes were higher in cluster 1 versus in cluster 2/3 in bladder cancer (p < 0.05). mRNAsi score was suggested to be correlated with the prognosis of cancer [40,41]. Our result revealed that the mRNAsi score in cluster 3 is higher than cluster 1 and cluster 2 ( Figure 4(D0, p ¼ 3.8e-14). TIDE score was good biomarker evaluating immunotherapy response [42]. Interestingly, the data suggested that cluster 3 of bladder cancer was correlated with lower TIDE score (Figure 4(E), p ¼ 4.2e-22). As chemotherapy and molecular targeted therapy were vital therapeutic strategies for bladder cancer, we also detected the IC50 score of common drugs of chemotherapy and molecular targeted therapy. Overall, the IC50 score of Vinorelbine, Mitomycin, Gemcitabine, Paclitaxel, Lapatinib, and Erlotinib in cluster 2 was lower than that in cluster 1/3 in bladder cancer ( Figure 5(A-F), all p < 0.05).

Development and validation of a prognostic signature for BLCA based on NETs-related genes
A total of five NET-related genes (CD93, CRISPLD2, KCNJ15, IRAK4, and MAPK3) were significantly correlated with the prognosis of bladder cancer (Figure 6(A)). Using these potential biomarkers, LASSO analysis was performed to  Supplementary table 2) and construct a NETsrelated prognostic signature. Interestingly, all these 5 genes were included in this prognostic signature based on the coefficient and partial likelihood deviance ( Figure 6(B,C)). After calculating the risk score of bladder cancer cases, we also compared the difference of high and low risk score in overall survival. As expected, high risk patients were correlated with poor clinical outcome in bladder cancer ( Figure 6(D), p < 0.001). Figure 6(E) showed the riskscore, survival status, and gene expression of prognostic signature in TCGA cohort. Further ROC analysis suggested that the AUC of 1-year, 3-year and 5-year ROC curve were 0.652, 0.640, and 0.670, respectively ( Figure 6(F)), demonstrating a good performance of this prognostic signature in prognosis prediction of bladder cancer. Interestingly, we also obtained similar results in GSE13507 ( Figure 6(G-I)) and GSE32894 dataset ( Figure 6(J-L)). We also constructed ROC curve and C-index analysis, suggested that risk score and clinical stage could serve as a good biomarker for prognosis prediction of bladder cancer in TCGA cohort (Figure 7(A,B)), GSE13507 (Figure 7(C,D)) and GSE32894 dataset ( Figure  7(E,F)). Interestingly, the result of univariate and multivariate Cox analysis also suggested NETsrelated risk score as an independent risk factor for bladder cancer (Figure 7(G,H)). We then analyzed the correlation between NETs-related risk score and immune cell infiltration. As shown in Figure  8(A-F), the level of M2 Macrophage, activated CD4þ memory T cell increased while the level of activated dendritic cell, resting dendritic cell, plasma cell, and follicular helper T cells decreased, as the riskscore increased.

Prognostic signature genes analysis
In order to clarify the role of prognostic signature genes in the progression of BLCA, we then analyzed their correlation with clinical characters. Significant difference in the expression of KCNJ15, CRISPLD2 and CD93 was obtained between different groups based on pT stage (Figure 9(A)), pN stage (Figure 9(B)), pM stage (Figure 9(C)) and pathological stage ( Figure  9(D)) analysis. Thus, these three genes were selected for further analysis. We found that all of three genes were involved in the regulation of some of famous cancer-related pathways ( Figure  9(E)). Drug sensitivity revealed that high CRISPLD2 expression and low CD93 were  associated with drug resistance of CTRP ( Figure  9(F)). As CD93 had been studied in some types of cancer [43,44], we selected CRISPLD2 for further analysis. To clarify the mechanism of CRISPLD2 in bladder cancer, we then explored CRISPLD2-associated lncRNA-miRNA-mRNA regulatory axis. As shown in Figure 10(A), four miRNAs (hsa-miR-6835-3p, hsa-miR-370-5p, hsa-miR-671-5p, and hsa-miR-493-5p) were suggested as the targets of CRISPLD2. However, only miR-370-5p was differentially expressed in bladder cancer and significantly correlated with the clinical outcome (Figure 10(B,C)). Thus, we selected miR-370-5p as the most promising target of CRISPLD2. As for the lncRNA targets of miR-370-5p, despite 12 lncRNAs were suggested as potential targets of miR-370-5p (Figure 10(D)), only AP4B1-AS1 was differentially expressed in bladder cancer and significantly correlated with the clinical outcome ( Figure 10(E,F)). Thus, we selected lncRNA AP4B1-AS1 as the most promising target of miR-370-5p. And we identified lncRNA AP4B1-AS1/miR-370-5p/CRISPLD2 regulatory axis for the progression of BLCA.

Discussion
NETs, the extracellular release of decondensed chromatin, exert a vital function in host defense against foreign pathogens [45]. The role of NETs in cancer is complex and it is far from fully being clarified. In most cases, NETs are suggested as pro-tumorigenic factors in cancer [45]. NET deposition is positively correlated with the tumor cell proliferation, immunosuppression and thrombosis [46,47]. Moreover, NETs can favore metastasis by facilitating ETM pathway [48]. NET components, including myeloperoxidase and histones, could directly kill tumor cells [18,19]. Though certain studies have been performed to clarify the role of NETs in BLCA, their function in BLCA was far from being clarified. Thus, our study is performed.
We first determined three subtypes of bladder cancer using consensus clustering analysis. Further analysis revealed that cluster 3 was correlated with a favorable outcome and cluster 1 had a poor overall survival in BLCA. mRNAsi score was suggested to be correlated with the prognosis of cancer [40,41]. Hang et al performed a bioinformatics study, revealing high mRNAsi score correlated with a favorable clinical outcome in colorectal cancer [49]. Another study suggested that gastric cancer patients with high mRNAsi had a better overall survival [50]. Our result revealed that the mRNAsi score in cluster 3 is higher than cluster 1/2 and BLCA patients in cluster 3 was correlated with a favorable clinical outcome, which was similar with previous study. TIDE score was a good biomarker evaluating immunotherapy response [42]. In our study, cluster 3 of bladder cancer was correlated with lower TIDE score, suggesting that bladder cancer patients in cluster 3 may be more sensitive to immunotherapy. Interestingly, bladder cancer patients in cluster 2 may be more sensitive to chemotherapy and target therapy, as the IC50 score of Vinorelbine, Mitomycin, Gemcitabine, Paclitaxel, Lapatinib, and Erlotinib in cluster 2 was lower than that bladder cancer patients in cluster 1/3. Using 5 NETs-related genes (CD93, CRISPLD2, KCNJ15, IRAK4, and MAPK3), we then constructed a prognostic signature for bladder cancer. As expected, this NETs-related prognostic signature had a good performance in predicting the prognosis of bladder cancer in TCGA cohort and GSE13507 and GSE32894 cohort. These results further highlighted the vital function of NETs in the prognosis of cancers. In fact, previous study had identified some NETsrelated signatures for certain types of cancer. Bioinformatics study had found that NETs score could act as a prognostic marker and be correlated with poor clinical outcome in most types of cancer [25]. Moreover, NETs-related lncRNA signature could predict the prognosis in lung cancer [51]. Based on NETs-related genes, Naifei ei al. also developed an innovative prognostic symbol, which could predict the prognosis and immunotherapy response in squamous cell carcinoma [52].
We also identified lncRNA AP4B1-AS1/miR-370-5p/CRISPLD2 regulatory axis that may be involved in the progression of bladder cancer. CRISPLD2 had anti-inflammatory property and was referred as a therapeutic agent in sepsis [53]. In most types of the cancer, CRISPLD2 showed positive correlation with HMGB1, which was correlated with DNA damage repair and genomic stability maintenance [53,54]. Interestingly, previous study revealed that miR-370-5p could regulate Figure 9. NETs-relative prognostic signature genes analysis. Gene expression in different groups based on pT stage(A), pN stage(B), pM stage(C) and pathological stage(D) analysis. The role of prognostic signature genes in famous cancer-associated pathways (E) and drug sensitivity (F). Ã p < 0.05, ÃÃ p < 0.01, ÃÃÃ p < 0.001. the proliferation, migration and invasion in bladder cancer by interacting with LINC01232 and PIM3 [55]. Moreover, miR-370-5P was associated with tumor progression in hepatocellular carcinoma by sponging Circular RNA circUBE2J2 [56]. However, limited study had been performed to clarify the role of AP4B1-AS1 in cancer. In our further study, we will focus on the validation of the role of lncRNA AP4B1-AS1/miR-370-5p/CRISPLD2 regulatory axis in bladder cancer by in vivo and in vitro studies.
The current study also had some limitations. The prognostic signature should be verified using more dataset. Moreover, in vivo and in vitro studies should be conducted to verify the role of lncRNA AP4B1-AS1/miR-370-5p/CRISPLD2 regulatory axis in bladder cancer.

Conclusion
Our study identified three molecular subtypes and a prognostic signature for bladder cancer based on NETs. These results may provide more data about the vital role of NETs in bladder cancer.

Declarations of interest
The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the article.

Data availability statement
The analyzed data sets generated during the study are available from the corresponding author on reasonable request.