Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Anoikis-related gene signature as novel prognostic biomarker to predict immunotherapy with bladder urothelial carcinoma

Anoikis-related gene signature as novel prognostic biomarker to predict immunotherapy with... IntroductionBLCA as a prevalent tumor of the urinary system, ranking in the top 10 of all malignancies. Although the survival proportion can reach 50 % for early-stage patients with bladder cancer, the survival rate for distant metastases is still less than 10 %. Recently, the rise of immunotherapy offers new options for BLCA treatment, but low response rates remain a current problem [1, 2]. So far, an adequate survival prediction model is still missing for BLCA. Therefore, it is significant to search for possible mechanisms of metastasis and explore robust multiple biomarkers to predict the clinical prognosis. Meanwhile, the identification of the new model is also important for effectively distinguishing risk groups and promoting immune precision therapy to improve clinical outcomes.Anoikis is regarded as a programmed form of cell death triggered by cells detached from the extracellular matrix (ECM), which plays a great role in the growth of our body, tissue self-equilibrium, disease progression, and tumor metastasis. ​Anoikis can protect the organism by preventing detached cells from re-adhering to other substrates to undergo abnormal proliferation. However, tumor cells are insensitive to anoikis, which adhere to substrates to survive and proliferate without ECM.Previous research results have demonstrated that anoikis-associated genes are important in tumor development and metastasis, such as ovarian cancer, melanoma, prostate cancer, and breast cancer. CBX2 is regarded as a contributing factor to anoikis resistance and spread in serous ovarian cancer [3]. Overexpression of BRN2 can induce anoikis resistance in melanoma [4]. TMPRSS4 promotes anti-anoikis, improving the growth of circulating tumor cells and facilitating initial metastasis of prostate cancer [5]. HBXIP promotes resistance to breast cancer anoikis by inhibiting the activation of JNK1 [6].Although anoikis has been shown to be related to the prognosis of several tumors, prognostic indicators related to anoikis are rarely mentioned in BLCA. Therefore, we investigated the prognostic relationship between anoikis genes and BLCA. In our work, our team constructed a robust and reliable genetic signature related to anoikis and explored its clinical prognostic and significance of immunotherapy in patients with BLCA.Materials and methodsData acquisition and processingThe transcriptional and survival materials of BLCA were obtained from GEO and TCGA databases and excluded patients of surviving time fewer than 30 days. The TCGA contained 412 tumor samples treated as a training cohort. The external independent dataset GSE13507 and GSE32894 were used as the test and validation cohort. For data processing, we used the mean value of multiple probes for expression when the gene was tagged with multiple probes. We also downloaded the external independent cohort of immunotherapy IMvigor210 online (http://research-pub.gene.com/IMvigor210CoreBiologies/). We presented the relative information of patients from the 4 datasets in Supplementary Table S1.Construction of the anoikis-associated signatureWe obtained genes associated with anoikis in the GeneCards database and set the correlation coefficient to be greater than 0.4 (Supplementary Table S2). We also extracted anoikis-associated differentially expressed genes (DEGs) set as p<0.05 and | log2 FC |>1.322. Then, univariate Cox hazard regression was used for choosing anoikis-associated prognostic value genes with a threshold of p<0.05. The R package ‘veen’ was conducted to find overlapped parts between DEGs and prognostic genes. The network was visualized with the R ‘circlize’ package. To better simulate the survival of the training queue, the selection operator and least absolute shrinkage regression were performed to build the risk model according to the ‘glmnet’ package. ​The choice of indicator lambda was based on cross-validation with k=10. High- and low-risk clusters were then grouped by the median of the model risk score.Risk score=∑k=1ncoef(genek)*expr(genek)$$Risk\hspace{0.17em}score=\sum\limits _{k=1}^{n}coef({gene}^{k})\ast expr({gene}^{k})$$The coef (genek) stood for a parameter of genes related to survival and expr (genek) was the expression of genes.Evaluation of ARGS in BLCAWe performed an independent prognostic analysis of the training, test, and validation groups in the risk model using the survival package. Besides, ROC was used to test the constructed model’s efficiency. Then multivariate and univariate Cox regressions were also applied to known independent risk factors in BLCA. We developed a nomogram diagram that included multiple clinicopathological characteristics and model risk scores to better simulate the prognostic value of patients with BLCA. We also drew a calibration curve to assess whether the predictions were in good agreement with reality.Enrichment analysis of pathways and functionsKyoto Encyclopedia of Genes and Genomes and Gene Ontology analysis on 21 overlapped genes were performed with package clusterProfiler. The GSEA analysis was then carried out using the software GSEA (https://www.gsea-msigdb.org/gsea/login.jsp:kegg.v.7.5.2). Meanwhile, we evaluated the degree of infiltration of T and macrophage M2 cells in ARGS using the CIBERSORT and XCELL algorithms.ARGS evaluation and response to immune checkpoint inhibitors (ICIs)Initially, the state of immune checkpoint activation was compared between different risk groups with the ggpubr package. Then immunophenoscore (IPS) was applied to observe the response rate of patients to ICIs, which was got in the TCIA atlas (https://tcia.at). An increased IPS usually reveals a better response to immunotherapy. We also analyzed cancer-associated CAF and immune exclusion data using the TIDE dataset (http://tide.dfci.harvard.edu/login/) to explore immune evasion mechanisms. A higher score on CAF and exclusion means more proper to escape from immunotherapy. The external independent cohort of immunotherapy IMvigor210 was further used to validate our model.Validation of the expression level of ARGS in BLCAWe performed survival analysis on the Kaplan-Meier Plotter database (Kaplan-Meier plotter (kmplot.com)) and cancer and adjacent tissue were obtained from the Second Hospital of Tianjin Medical University. The participants were aware of and signed an informed consent form. Our work was endorsed by the Ethics Committee of the Second Hospital of Tianjin Medical University. TRIzol reagents (Invitrogen, Shanghai, China) were utilized to separate and purify total RNA, and BIOG cDNA synthesis kits (BioDai, Changzhou, China) were applied to synthesize DNA. Then, qRT-PCR was applied with FastStart Universal SYBR Green Master (Roche (ROX), Mannheim, Germany). (GAPDH-F: GGAAGGTGAAGGTCGGAGTCA, GAPDH-R: GTCATTGATGGCAACAATATATCCACT, MYC-F: TGGAAAACCAGCCTCCCG, MYC-R: TTCTCCTCCTCGTCGCAGTA; RAC3-F: GGAGATTGGTGGGCTCTGTG, RAC3-R: TTCTTCCCCGGCTTCTTCAC; S100A7-F: GCTTCCCAGCTCTGGCTTTT, S100A7-R: GCAGGCTTGGCTTCTCAATC).ResultsDifferential expression level of anoikis genes in BLCATo show our research process more clearly, we made a flowchart and showed it as follows (Figure 1). We downloaded the BLCA expression matrix and information on clinicopathological parameters from the TCGA dataset and performed a prognostic analysis based on 434 anoikis related genes. 101 genes were significantly correlated with prognosis (Figure 2A and Supplementary Table S3). Then differential analysis (|Log2FC|>1.322, p<0.05) based on 434 genes was also performed. We screened out 56 differential genes (DEGs) (Figure 2B and Supplementary Table S4) and 21 of them were overlapped DEGs (Figure 2C). Furthermore, most of the 21 genes were positively related to each other (Figure 2D). KEGG and GO analyses were also conducted to probe the biological processes of 21 genes. The results showed that they were mainly enriched in cell junction, adhesion, metastasis, and proliferation (Supplementary Figure S1).Figure 1:Research flow chart.Figure 2:Construction of anoikis-related gene signature in BLCA patients. The forest plot of 101 anoikis-related genes with prognostic relevance (A). The volcano plot of 56 anoikis-related genes with different expression (|Log2FC|>1.322, p<0.05) (B). 21 overlap genes in DEGs and prognostic genes are represented in a venn diagram (C). 21 gene correlation network graph (D).Construction of the cellular anoikis-related signature in BLCANext, we built an ARGS for survival forecasting. To eliminate overlapping of the ARGS model, lasso regression was carried out on 21 screened genes above, and 15 ARGS in BLCA were targeted according to the first level value of Log(λ) (Figure 3A and B). Hazard scores were obtained by Multivariate cox regression based on the following formula: hazard score=CAV1 × (−0.0865) + ITGA5 × (−0.0426) + GSPG4 × (0.1723) + MYC × (0.0961) + NTF3 × (−0.4377) + ZEB2 × (−0.2730) + CRYAB × (0.0297) + ZEB1 × (0.3589) + F10 × (0.2523) + SFRP1 × (−0.0146) + TAGLN × (−0.0884) + CCDC80 × (0.1880) + ID2× (−0.0982) + RAC3 × (0.2285) + S100A7 × (0.0413).Figure 3:Prognostic evaluation of 15 anoikis-related gene models. 10-Fold counter verification and variable selection for LASSO regression (A). The overview of the LASSO parameters of 15 anoikis-related genes (B). Survival analysis of TCGA (C) and GEO (D). ROC of the hazard score and clinicopathological parameters (E). 5 years ROC curve of the prognostic model (F). The distribution of hazard score (G), survival outcomes (I), and the expression profile of 15 genes with a heatmap (K) from the TCGA database. The distribution of ARGS hazard score (H), survival outcomes (J), and the expression profile of 15 genes with a heatmap (L) from the GEO database. Prognostic value of ARGS in tumor stage (M-N). Tumor grade between the two groups (O). Prognostic value in the GSE32894 independent cohort (P).According to the results of the hazard score formula, patients were grouped into different risk groups depending on median hazard scores. Patients in TCGA as the training group and patients in GSE32894 as the test group. We evaluated survival probability (Figure 3C and D). The outcomes indicated that the high-risk cluster showed a poorer prognosis. In addition, the sensitivity and specificity of ARGS were evaluated with ROC and the risk score showed a primary predictive value compared to other clinical factors (Figure 3E). Meanwhile, the proportion under ROC for OS 1, 3, and 5 years was 0.735, 0.701, and 0.688 (Figure 3F). Then the layout of the hazard scores for ARGS, patient prognosis status, and expression profiles of 15 genes in different risk groups were also explored (Figure 3G–L). The results also showed that the high-risk group represented a worse prognosis. Furthermore, the TNM stage also showed the same results, which further proved the precision of the ARGS model (Figure 3M and N). We also analyzed the association between ARGS and tumor grade and found that high-grade tumors had higher risk scores, indicating that higher risk scores were related to the malignant progression of BLCA (Figure 3O). Finally, we again validated the survival value of our model in the external bladder cancer queue GSE32894 (Figure 3P).Construction of nomogramThen multivariate and univariate Cox regression were applied to explore whether ARG hazard scores could be considered as a specific risk indicator. The univariate regression demonstrated that TNM stage, hazard scores, and age were all related to survival (Figure 4A). Next, multivariate analysis revealed that hazard scores were still associated with survival with the highest hazard ratio (HR=3.239, 95%CI:2.291–4.578, p<0.001) (Figure 4B). To probe the survival value of ARGS, we constructed a nomogram in combination with other clinical indicators and found that the survival percentage of patients at 1, 3, and 5 years was 0.909, 0.711, and 0.606 under the risk model. In the 5-year calibration curve, we noticed that 1,3, and 5-year prediction curves had less deviation and therefore could be better used to predict the patient’s prognosis (Figure 4C and D).Figure 4:Nomogram of the ARGS risk model. Multi- and univariate regression (A-B) to assess independent risk indicators of ARGS. The nomogram includes the hazard score and other clinical parameters to forecast 1, 3, and 5-year survival rates (C). 1, 3, and 5-year overall survive calibration curves (D).Analysis of the biological function of ARGSOur team also used GSEA software to investigate biological activities in the model of patients. We found that CHEMOKINE signaling, HEDGEHOG signaling, WNT signaling, MAPK signaling, VEGF signaling, TGF BETA signaling, chemokine signaling, P53 signaling, and CYTOKINE signaling were significantly enriched in the high-risk cluster (Figure 5A and Supplementary Figure S2). The transmission of the cellular signaling pathway is based on secreted proteins that are closely related to disease development and play an essential role in biological responses and homeostasis [7]. Additionally, we also noticed that several secreted proteins were expressed in greater amounts in the high-risk cluster compared to the low-risk cluster (Figure 5B). Importantly, some secreted proteins can induce an immunosuppressive microenvironment, such as the CXCL family, TNFRSF11B, and IL6 [8], [9], [10]. In addition, the WNT pathway can induce some tumors with non-T cell phenotypes and reduce the prognosis of immunotherapy [11]. TGF‐β signaling pathway modifies TME to suppress T cell infiltration and limit antitumor immunity [12]. Therefore, we speculated that the high-risk group tends to be an immunosuppressive phenotype. Next, we probed the relevance between risk scores and T cells or TAM infiltration using CIBERSORT and XCELL algorithms. ​We observed a favorable correlation between TAM and ARGS scores, while T cell infiltration was negatively correlated with ARGS scores (Figure 5C).Figure 5:Analysis of the biological function and immune status of ARGS. Analysis of multiple pathways in different risk groups (A). Comparison of multiple secreted proteins between two risk clusters (B). The relationship between ARGS risk scores and inflation levels of T cells and macrophageM2 (C).The significance of immune checkpoints in ARGSImmune checkpoints are critically important for modulating the TME, according to previous research [13]. Hence, we need to better understand the relationship between ARGS and immune checkpoints. The result of the barplot showed that most immune checkpoints were highly expressed in the ARGS high-risk group, including CTLA4, CD274\PD-L1, CD276, and TIGIT (Figure 6A). This suggested that high-risk ARGS cluster might be involved in the immune resistance of BLCA patients.Figure 6:Distinct expression of the immune checkpoints in two risk groups (A). (*: p<0.05, **: p<0.01, and ***: p<0.001).The study of ARGS in immunotherapy responseThe practice of ICIs in uroepithelial tumors has been successful but objective response rates are not as good as expected [2]. The efficacy of immunotherapy is influenced by many factors, especially CAF. CAF can express suppressive immune checkpoint signaling and immunosuppressive factors to further inhibit the ability of various immune effector cells [14]. We discovered that CAF was enriched in the high-risk ARGS cluster and showed immune exclusion (Figure 7A and B). Next, we further analyzed the predictive potential of ARGS in the direction of immunotherapy. As illustrated in Figure 7C and D, the low-risk ARGS cluster had a higher IPS and anti-CTLA4 treatment was relatively sensitive to the low-risk cluster. Based on these results, the low-risk cluster would benefit from immunotherapy. Then, we utilized the immunotherapy cohort IMvigor210 to test the forecast power of immunotherapy. We noticed that the ARGS low-risk cluster inclined to a better prognosis and the high-risk cluster showed significantly poorer results and lower response rates (Figure 7E and F). Next, we explored the prognostic value of the model combined with the immune checkpoints. The results revealed that low-risk ARGS cluster with a higher CTLA4 expression level had longer survival time and a better prognosis (Figure 7G). The above results suggested that patients with higher CTLA4 expression in the low-risk ARGS group would have a better clinical benefit when receiving the corresponding immunotherapy.Figure 7:Evaluation of immunotherapy of the ARGS model. The distribution of CAF scores (A), exclusion scores (B), and IPS scores (C–D) in the two groups. Prognosis and potential predictivity of the immunotherapy response of the ARGS model in the IMvigor cohort. CR/PR: Complete/partial response, PD/SD: Progressive/stable disease (E–G).External test of anoikis-associated genes as a latent biomarkerWe utilized the Kaplan-Meier Plotter database to verify the clinical prognosis of MYC, RAC3, and S100A7. The plots showed that high expression of MYC (HR=1.81 (1.35–2.43), p=5.8e-05) (Figure 8A), RAC3 (HR=1.57 (1.17–2.11), p=0.0024) (Figure 8C) and S100A7 (HR=1.5 (1.1–2.05), p=0.01) (Figure 8E) was associated with a worse prognosis. We then carried out in vitro experiments to test the expression trends of the above genes in BLCA. RT-qPCR also showed that MYC (Figure 8B), RAC3 (Figure 8D), and S100A7 (Figure 8F) had a greater amount of expression in BLCA tissues compared to adjacent tissues.Figure 8:External test of anoikis-related genes as latent biomarkers. KM survival analysis of MYC (A), RAC3 (C), and S100A7 (E). Vitro RT-qPCR of MYC (B), RAC3 (D) and S100A7 (F).DiscussionBladder urothelial carcinoma is a common urological carcinoma with a trend of increasing mortality and morbidity. With the treatment of conventional therapy, such as surgery, radiation, chemotherapy, and bacillus calmette guerin (BCG), the survival time for BLCA patients only has 14 months [15, 16]. Although the use of ICIs for urothelial cancers has also achieved great success, the response rate is only about 20–23 % [2]. Therefore, we need to develop a reliable prognostic indicator and achieve precise treatment for BLCA. Anoikis, a form of programmed cell death, can regulate a variety of tumor biological behaviors. Silencing of ITGA5 in vitro can further promote anoikis of hepatocellular carcinoma (HCC) and suppress the development of tumors [17]. MYH9 also promotes the escape of anoikis by inducing CTNNB1 transcription in gastric cancer [18]. Overexpressed ZNF32 can also promote resistance to anoikis in HCC by preserving redox homeostasis and promoting Src/FAK signaling [19]. Therefore, anoikis-related genes show promise as new targets for tumor therapy.To obtain a reliable ARGS for BLCA, TCGA was used as the training cohort, and the external independent database GEO as the test and validation cohort. To avoid overfitting, we performed a LASSO regression to generate a 15-gene prognostic model that had significant prognostic outcomes in the TCGA and GEO databases. To test the performance of the ARGS, our team built a norman plot combining age, grade, stage, and risk scores. The calibration plots was then used to suggest good discriminative and goodness-of-fit predictions for the prognosis. Zhang’s article built a model for bladder cancer [20]. The areas under the ROC of Zhang’s signature were 0.667, 0.632, and 0.637 for 1, 3, and 5 years, all of which were lower than the areas under the ROC curves of our model. The results demonstrated that our model was more accurate in assessing the survival ratio of patients with 1, 3, and 5 years.The ARGS model we constructed was significantly associated with the survival outcome of BLCA. ARGS contained 15 genes, including CAV1, ITGA5, CSPG4, MYC, NTF3, ZEB2, CRYAB, ZEB1, F10, SFRP1, TAGLN, CCDC80, ID2, RAC3, S100A7. All genes included in the ARGS model have been reported to be associated with tumors. CAV1 is involved in CD26-mediated angiogenesis and metastasis in colorectal cancer [21]. ITGA5, a member of the integrin family, facilitates communication among several cells and cells with ECM. Especially ITGA5 can promote the development of colorectal cancer through its own O-GlcNAcylation [22]. CSPG4 has been investigated to regulate tumor growth and invasion via RTK signaling, MAPK resistance, and integrin signaling by promoting focal adhesion kinase activation [23]. MYC is involved in the CAF oncogenic process in the tumor microenvironment and contributes to tumor proliferation by hindering the PI3K/Akt/β-Catenin signaling [24]. NTF3 as a new target promotes apoptosis in hepatocellular cancer through the signaling of P38 MAPK and JNK signaling [25]. In lung carcinoma, the PAX6-ZEB2 pathway increases cisplatin resistance and migration via PI3K/AKT signaling [26]. CRYAB promotes prostate cancer progression and metastasis as a downstream target of MiRNA-671-5p targeting NFIA [27]. ZEB1 can promote the warburg effect through the activation of PFKM in HCC [28]. Hypermethylation of the F10 promoter in gliomas is an independent prognostic factor and further increases tumor aggressiveness [29]. SFRP1 can affect the proliferation of tumor with Wnt pathway [30]. In colorectal cancer, TAGLN mediates ovarian cancer progression by regulating stromal stiffness through the RhoA/ROCK pathway [31]. CCDC80 plays an important role in melanoma migration via the FAK/CCDC80/E-cadherin axis [32]. ID2 has been investigated to accelerate the aggressiveness of head and neck squamous cell carcinoma via the ID2-MMP1 axis [33]. RAC3 leads to colorectal cancer (CRC) progression and stemming of cancer cells by acting on its downstream target CTFR [34]. S100A7 facilitates esophageal squamous cell carcinoma (ESCC) growth and migration through the incorporation of intracellular JAB1 and increases tumor angiogenesis through the promotion of p-FAK and p-ErK signaling [35]. To explore the biological signaling pathways of ARGS, we performed a GSEA analysis and identified multiple signaling pathways that were enriched in the high-risk group, including hedgehog, MAPK, p53, TGF β, VEGF, and WNT signaling. The above regulation paths have been reported to contribute to BLCA growth. Propofol suppresses the development and stem cellularity of BLCA by inhibiting hedgehog signaling [36]. Mir-186-5p acts on Ras-related protein RAB27 A/B to inhibit BLCA development and EMT process through PI3K/MAPK signaling pathway inactivation [37]. As a cancer suppressor gene, UNC5B-mediated dephosphorylation of p53 in bladder cancer can arrest the cell cycle and inhibit tumor proliferation, providing a therapeutic target for BLCA patients [38]. LncRNA PLAC2 can increase miR-663 which inhibits TGF-β1 to reduce BLCA invasion and metastasis [39]. Silencing RIPK4 will significantly reduce the expression level of VEGF-A to suppress tumor development [40]. TMEM88 can inhibit proliferation and invasion outcomes in BLCA via reducing Wnt/β-catenin signaling [41].Immunotherapy has potency as a new option for BLCA treatment [42]. TME contains a wide range of immune and non-immune components, creating a chronic inflammatory, pro-tumor environment, and immunosuppression [43]. Meng’s article also reported that stromal cells had a great effect on the modification of TME. Stromal cells can suppress the immune process by increasing the secretion of the immunosuppressive factor IL10 to create an immunosuppressive microenvironment [44]. In the risk model, IL10 was substantially enriched in the high-risk ARGS group, indicating that the high-risk ARGS group might have a suppressive immune microenvironment. TAM can shape TME and affect the efficacy of immunotherapy [45]. Higher TAM can predict unfavorable survival outcomes and evaluate the sensitivity of PD-L1 immunotherapy in patients with BLCA [46, 47]. CAF, as the main component of TME, increases BLCA aggressiveness through paracrine IL-6 signaling [48]. Furthermore, CAF can also induce the formation of an immunosuppressive tumor microenvironment [49]. Meng’s article illustrated that the immune-depleted subgroup had a fibroblast-regulated cellular matrix enriched with WNT/TGF-β and LAG3 checkpoint [44]. Our TIDE, GSEA, and checkpoint analysis outcomes also suggested that the high-risk ARGS cluster was concentrated fibroblasts, WNT/TGF-β, and a high level of LAG3 according to the characteristics of the immune-depleted subgroup. The above findings clarified that the high-risk group was an immunosuppressed state.We also explored whether ARGS could identify the sensitivity of patients with BLCA to immunotherapy. The immune checkpoint (ICI) is one of the reliable indicators for evaluating immunotherapy. Our results showed that PD-L1, TIGIT, CD276, and CTLA4 had a high expression level in the high-risk cluster. Additionally, the high ARGS risk group also had a stronger immune exclusion. The TCIA results also showed a better response to ICIs in the low-risk cluster. The above results suggested that the high ARGS risk cluster would not benefit from immunotherapy. We also validated our findings in the immune cohort IMvigor210 with the same results. Overall, we have identified a subtype of BLCA with a better prognosis and sensitivity to immunotherapy. These results will facilitate the improvement of clinical outcomes for BLCA patients and enable precision therapy.ConclusionsIn general, we constructed and validated an anoikis-associated genetic signature based on 15 anoikis-associated genes as a predictor of immunotherapy with a significant prognosis for BLCA patients. Furthermore, ARGS is involved in the regulation of TME and contributes to the development of the immunosuppressive microenvironment. Last, we also explored the combined role of ARGS scores and immune checkpoint genes in patients with BLCA. Based on the above findings, ​we believe that the combination of ARGS scores with specific immune checkpoints has the potential to be a predictive indicator of ICIs response, leading to more precise treatment and increases patient benefit. Therefore, identifying the anoikis-associated gene signature can help us risk stratify and achieve precision therapy, providing new insights to improve the response of BLCA to immunotherapy. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png ONCOLOGIE de Gruyter

Anoikis-related gene signature as novel prognostic biomarker to predict immunotherapy with bladder urothelial carcinoma

Download PDF
18 pages

Loading next page...
 
/lp/de-gruyter/anoikis-related-gene-signature-as-novel-prognostic-biomarker-to-LvMtq1gwFQ

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
de Gruyter
Copyright
© 2023 the author(s), published by De Gruyter, Berlin/Boston
eISSN
1765-2839
DOI
10.1515/oncologie-2023-0049
Publisher site
See Article on Publisher Site

Abstract

IntroductionBLCA as a prevalent tumor of the urinary system, ranking in the top 10 of all malignancies. Although the survival proportion can reach 50 % for early-stage patients with bladder cancer, the survival rate for distant metastases is still less than 10 %. Recently, the rise of immunotherapy offers new options for BLCA treatment, but low response rates remain a current problem [1, 2]. So far, an adequate survival prediction model is still missing for BLCA. Therefore, it is significant to search for possible mechanisms of metastasis and explore robust multiple biomarkers to predict the clinical prognosis. Meanwhile, the identification of the new model is also important for effectively distinguishing risk groups and promoting immune precision therapy to improve clinical outcomes.Anoikis is regarded as a programmed form of cell death triggered by cells detached from the extracellular matrix (ECM), which plays a great role in the growth of our body, tissue self-equilibrium, disease progression, and tumor metastasis. ​Anoikis can protect the organism by preventing detached cells from re-adhering to other substrates to undergo abnormal proliferation. However, tumor cells are insensitive to anoikis, which adhere to substrates to survive and proliferate without ECM.Previous research results have demonstrated that anoikis-associated genes are important in tumor development and metastasis, such as ovarian cancer, melanoma, prostate cancer, and breast cancer. CBX2 is regarded as a contributing factor to anoikis resistance and spread in serous ovarian cancer [3]. Overexpression of BRN2 can induce anoikis resistance in melanoma [4]. TMPRSS4 promotes anti-anoikis, improving the growth of circulating tumor cells and facilitating initial metastasis of prostate cancer [5]. HBXIP promotes resistance to breast cancer anoikis by inhibiting the activation of JNK1 [6].Although anoikis has been shown to be related to the prognosis of several tumors, prognostic indicators related to anoikis are rarely mentioned in BLCA. Therefore, we investigated the prognostic relationship between anoikis genes and BLCA. In our work, our team constructed a robust and reliable genetic signature related to anoikis and explored its clinical prognostic and significance of immunotherapy in patients with BLCA.Materials and methodsData acquisition and processingThe transcriptional and survival materials of BLCA were obtained from GEO and TCGA databases and excluded patients of surviving time fewer than 30 days. The TCGA contained 412 tumor samples treated as a training cohort. The external independent dataset GSE13507 and GSE32894 were used as the test and validation cohort. For data processing, we used the mean value of multiple probes for expression when the gene was tagged with multiple probes. We also downloaded the external independent cohort of immunotherapy IMvigor210 online (http://research-pub.gene.com/IMvigor210CoreBiologies/). We presented the relative information of patients from the 4 datasets in Supplementary Table S1.Construction of the anoikis-associated signatureWe obtained genes associated with anoikis in the GeneCards database and set the correlation coefficient to be greater than 0.4 (Supplementary Table S2). We also extracted anoikis-associated differentially expressed genes (DEGs) set as p<0.05 and | log2 FC |>1.322. Then, univariate Cox hazard regression was used for choosing anoikis-associated prognostic value genes with a threshold of p<0.05. The R package ‘veen’ was conducted to find overlapped parts between DEGs and prognostic genes. The network was visualized with the R ‘circlize’ package. To better simulate the survival of the training queue, the selection operator and least absolute shrinkage regression were performed to build the risk model according to the ‘glmnet’ package. ​The choice of indicator lambda was based on cross-validation with k=10. High- and low-risk clusters were then grouped by the median of the model risk score.Risk score=∑k=1ncoef(genek)*expr(genek)$$Risk\hspace{0.17em}score=\sum\limits _{k=1}^{n}coef({gene}^{k})\ast expr({gene}^{k})$$The coef (genek) stood for a parameter of genes related to survival and expr (genek) was the expression of genes.Evaluation of ARGS in BLCAWe performed an independent prognostic analysis of the training, test, and validation groups in the risk model using the survival package. Besides, ROC was used to test the constructed model’s efficiency. Then multivariate and univariate Cox regressions were also applied to known independent risk factors in BLCA. We developed a nomogram diagram that included multiple clinicopathological characteristics and model risk scores to better simulate the prognostic value of patients with BLCA. We also drew a calibration curve to assess whether the predictions were in good agreement with reality.Enrichment analysis of pathways and functionsKyoto Encyclopedia of Genes and Genomes and Gene Ontology analysis on 21 overlapped genes were performed with package clusterProfiler. The GSEA analysis was then carried out using the software GSEA (https://www.gsea-msigdb.org/gsea/login.jsp:kegg.v.7.5.2). Meanwhile, we evaluated the degree of infiltration of T and macrophage M2 cells in ARGS using the CIBERSORT and XCELL algorithms.ARGS evaluation and response to immune checkpoint inhibitors (ICIs)Initially, the state of immune checkpoint activation was compared between different risk groups with the ggpubr package. Then immunophenoscore (IPS) was applied to observe the response rate of patients to ICIs, which was got in the TCIA atlas (https://tcia.at). An increased IPS usually reveals a better response to immunotherapy. We also analyzed cancer-associated CAF and immune exclusion data using the TIDE dataset (http://tide.dfci.harvard.edu/login/) to explore immune evasion mechanisms. A higher score on CAF and exclusion means more proper to escape from immunotherapy. The external independent cohort of immunotherapy IMvigor210 was further used to validate our model.Validation of the expression level of ARGS in BLCAWe performed survival analysis on the Kaplan-Meier Plotter database (Kaplan-Meier plotter (kmplot.com)) and cancer and adjacent tissue were obtained from the Second Hospital of Tianjin Medical University. The participants were aware of and signed an informed consent form. Our work was endorsed by the Ethics Committee of the Second Hospital of Tianjin Medical University. TRIzol reagents (Invitrogen, Shanghai, China) were utilized to separate and purify total RNA, and BIOG cDNA synthesis kits (BioDai, Changzhou, China) were applied to synthesize DNA. Then, qRT-PCR was applied with FastStart Universal SYBR Green Master (Roche (ROX), Mannheim, Germany). (GAPDH-F: GGAAGGTGAAGGTCGGAGTCA, GAPDH-R: GTCATTGATGGCAACAATATATCCACT, MYC-F: TGGAAAACCAGCCTCCCG, MYC-R: TTCTCCTCCTCGTCGCAGTA; RAC3-F: GGAGATTGGTGGGCTCTGTG, RAC3-R: TTCTTCCCCGGCTTCTTCAC; S100A7-F: GCTTCCCAGCTCTGGCTTTT, S100A7-R: GCAGGCTTGGCTTCTCAATC).ResultsDifferential expression level of anoikis genes in BLCATo show our research process more clearly, we made a flowchart and showed it as follows (Figure 1). We downloaded the BLCA expression matrix and information on clinicopathological parameters from the TCGA dataset and performed a prognostic analysis based on 434 anoikis related genes. 101 genes were significantly correlated with prognosis (Figure 2A and Supplementary Table S3). Then differential analysis (|Log2FC|>1.322, p<0.05) based on 434 genes was also performed. We screened out 56 differential genes (DEGs) (Figure 2B and Supplementary Table S4) and 21 of them were overlapped DEGs (Figure 2C). Furthermore, most of the 21 genes were positively related to each other (Figure 2D). KEGG and GO analyses were also conducted to probe the biological processes of 21 genes. The results showed that they were mainly enriched in cell junction, adhesion, metastasis, and proliferation (Supplementary Figure S1).Figure 1:Research flow chart.Figure 2:Construction of anoikis-related gene signature in BLCA patients. The forest plot of 101 anoikis-related genes with prognostic relevance (A). The volcano plot of 56 anoikis-related genes with different expression (|Log2FC|>1.322, p<0.05) (B). 21 overlap genes in DEGs and prognostic genes are represented in a venn diagram (C). 21 gene correlation network graph (D).Construction of the cellular anoikis-related signature in BLCANext, we built an ARGS for survival forecasting. To eliminate overlapping of the ARGS model, lasso regression was carried out on 21 screened genes above, and 15 ARGS in BLCA were targeted according to the first level value of Log(λ) (Figure 3A and B). Hazard scores were obtained by Multivariate cox regression based on the following formula: hazard score=CAV1 × (−0.0865) + ITGA5 × (−0.0426) + GSPG4 × (0.1723) + MYC × (0.0961) + NTF3 × (−0.4377) + ZEB2 × (−0.2730) + CRYAB × (0.0297) + ZEB1 × (0.3589) + F10 × (0.2523) + SFRP1 × (−0.0146) + TAGLN × (−0.0884) + CCDC80 × (0.1880) + ID2× (−0.0982) + RAC3 × (0.2285) + S100A7 × (0.0413).Figure 3:Prognostic evaluation of 15 anoikis-related gene models. 10-Fold counter verification and variable selection for LASSO regression (A). The overview of the LASSO parameters of 15 anoikis-related genes (B). Survival analysis of TCGA (C) and GEO (D). ROC of the hazard score and clinicopathological parameters (E). 5 years ROC curve of the prognostic model (F). The distribution of hazard score (G), survival outcomes (I), and the expression profile of 15 genes with a heatmap (K) from the TCGA database. The distribution of ARGS hazard score (H), survival outcomes (J), and the expression profile of 15 genes with a heatmap (L) from the GEO database. Prognostic value of ARGS in tumor stage (M-N). Tumor grade between the two groups (O). Prognostic value in the GSE32894 independent cohort (P).According to the results of the hazard score formula, patients were grouped into different risk groups depending on median hazard scores. Patients in TCGA as the training group and patients in GSE32894 as the test group. We evaluated survival probability (Figure 3C and D). The outcomes indicated that the high-risk cluster showed a poorer prognosis. In addition, the sensitivity and specificity of ARGS were evaluated with ROC and the risk score showed a primary predictive value compared to other clinical factors (Figure 3E). Meanwhile, the proportion under ROC for OS 1, 3, and 5 years was 0.735, 0.701, and 0.688 (Figure 3F). Then the layout of the hazard scores for ARGS, patient prognosis status, and expression profiles of 15 genes in different risk groups were also explored (Figure 3G–L). The results also showed that the high-risk group represented a worse prognosis. Furthermore, the TNM stage also showed the same results, which further proved the precision of the ARGS model (Figure 3M and N). We also analyzed the association between ARGS and tumor grade and found that high-grade tumors had higher risk scores, indicating that higher risk scores were related to the malignant progression of BLCA (Figure 3O). Finally, we again validated the survival value of our model in the external bladder cancer queue GSE32894 (Figure 3P).Construction of nomogramThen multivariate and univariate Cox regression were applied to explore whether ARG hazard scores could be considered as a specific risk indicator. The univariate regression demonstrated that TNM stage, hazard scores, and age were all related to survival (Figure 4A). Next, multivariate analysis revealed that hazard scores were still associated with survival with the highest hazard ratio (HR=3.239, 95%CI:2.291–4.578, p<0.001) (Figure 4B). To probe the survival value of ARGS, we constructed a nomogram in combination with other clinical indicators and found that the survival percentage of patients at 1, 3, and 5 years was 0.909, 0.711, and 0.606 under the risk model. In the 5-year calibration curve, we noticed that 1,3, and 5-year prediction curves had less deviation and therefore could be better used to predict the patient’s prognosis (Figure 4C and D).Figure 4:Nomogram of the ARGS risk model. Multi- and univariate regression (A-B) to assess independent risk indicators of ARGS. The nomogram includes the hazard score and other clinical parameters to forecast 1, 3, and 5-year survival rates (C). 1, 3, and 5-year overall survive calibration curves (D).Analysis of the biological function of ARGSOur team also used GSEA software to investigate biological activities in the model of patients. We found that CHEMOKINE signaling, HEDGEHOG signaling, WNT signaling, MAPK signaling, VEGF signaling, TGF BETA signaling, chemokine signaling, P53 signaling, and CYTOKINE signaling were significantly enriched in the high-risk cluster (Figure 5A and Supplementary Figure S2). The transmission of the cellular signaling pathway is based on secreted proteins that are closely related to disease development and play an essential role in biological responses and homeostasis [7]. Additionally, we also noticed that several secreted proteins were expressed in greater amounts in the high-risk cluster compared to the low-risk cluster (Figure 5B). Importantly, some secreted proteins can induce an immunosuppressive microenvironment, such as the CXCL family, TNFRSF11B, and IL6 [8], [9], [10]. In addition, the WNT pathway can induce some tumors with non-T cell phenotypes and reduce the prognosis of immunotherapy [11]. TGF‐β signaling pathway modifies TME to suppress T cell infiltration and limit antitumor immunity [12]. Therefore, we speculated that the high-risk group tends to be an immunosuppressive phenotype. Next, we probed the relevance between risk scores and T cells or TAM infiltration using CIBERSORT and XCELL algorithms. ​We observed a favorable correlation between TAM and ARGS scores, while T cell infiltration was negatively correlated with ARGS scores (Figure 5C).Figure 5:Analysis of the biological function and immune status of ARGS. Analysis of multiple pathways in different risk groups (A). Comparison of multiple secreted proteins between two risk clusters (B). The relationship between ARGS risk scores and inflation levels of T cells and macrophageM2 (C).The significance of immune checkpoints in ARGSImmune checkpoints are critically important for modulating the TME, according to previous research [13]. Hence, we need to better understand the relationship between ARGS and immune checkpoints. The result of the barplot showed that most immune checkpoints were highly expressed in the ARGS high-risk group, including CTLA4, CD274\PD-L1, CD276, and TIGIT (Figure 6A). This suggested that high-risk ARGS cluster might be involved in the immune resistance of BLCA patients.Figure 6:Distinct expression of the immune checkpoints in two risk groups (A). (*: p<0.05, **: p<0.01, and ***: p<0.001).The study of ARGS in immunotherapy responseThe practice of ICIs in uroepithelial tumors has been successful but objective response rates are not as good as expected [2]. The efficacy of immunotherapy is influenced by many factors, especially CAF. CAF can express suppressive immune checkpoint signaling and immunosuppressive factors to further inhibit the ability of various immune effector cells [14]. We discovered that CAF was enriched in the high-risk ARGS cluster and showed immune exclusion (Figure 7A and B). Next, we further analyzed the predictive potential of ARGS in the direction of immunotherapy. As illustrated in Figure 7C and D, the low-risk ARGS cluster had a higher IPS and anti-CTLA4 treatment was relatively sensitive to the low-risk cluster. Based on these results, the low-risk cluster would benefit from immunotherapy. Then, we utilized the immunotherapy cohort IMvigor210 to test the forecast power of immunotherapy. We noticed that the ARGS low-risk cluster inclined to a better prognosis and the high-risk cluster showed significantly poorer results and lower response rates (Figure 7E and F). Next, we explored the prognostic value of the model combined with the immune checkpoints. The results revealed that low-risk ARGS cluster with a higher CTLA4 expression level had longer survival time and a better prognosis (Figure 7G). The above results suggested that patients with higher CTLA4 expression in the low-risk ARGS group would have a better clinical benefit when receiving the corresponding immunotherapy.Figure 7:Evaluation of immunotherapy of the ARGS model. The distribution of CAF scores (A), exclusion scores (B), and IPS scores (C–D) in the two groups. Prognosis and potential predictivity of the immunotherapy response of the ARGS model in the IMvigor cohort. CR/PR: Complete/partial response, PD/SD: Progressive/stable disease (E–G).External test of anoikis-associated genes as a latent biomarkerWe utilized the Kaplan-Meier Plotter database to verify the clinical prognosis of MYC, RAC3, and S100A7. The plots showed that high expression of MYC (HR=1.81 (1.35–2.43), p=5.8e-05) (Figure 8A), RAC3 (HR=1.57 (1.17–2.11), p=0.0024) (Figure 8C) and S100A7 (HR=1.5 (1.1–2.05), p=0.01) (Figure 8E) was associated with a worse prognosis. We then carried out in vitro experiments to test the expression trends of the above genes in BLCA. RT-qPCR also showed that MYC (Figure 8B), RAC3 (Figure 8D), and S100A7 (Figure 8F) had a greater amount of expression in BLCA tissues compared to adjacent tissues.Figure 8:External test of anoikis-related genes as latent biomarkers. KM survival analysis of MYC (A), RAC3 (C), and S100A7 (E). Vitro RT-qPCR of MYC (B), RAC3 (D) and S100A7 (F).DiscussionBladder urothelial carcinoma is a common urological carcinoma with a trend of increasing mortality and morbidity. With the treatment of conventional therapy, such as surgery, radiation, chemotherapy, and bacillus calmette guerin (BCG), the survival time for BLCA patients only has 14 months [15, 16]. Although the use of ICIs for urothelial cancers has also achieved great success, the response rate is only about 20–23 % [2]. Therefore, we need to develop a reliable prognostic indicator and achieve precise treatment for BLCA. Anoikis, a form of programmed cell death, can regulate a variety of tumor biological behaviors. Silencing of ITGA5 in vitro can further promote anoikis of hepatocellular carcinoma (HCC) and suppress the development of tumors [17]. MYH9 also promotes the escape of anoikis by inducing CTNNB1 transcription in gastric cancer [18]. Overexpressed ZNF32 can also promote resistance to anoikis in HCC by preserving redox homeostasis and promoting Src/FAK signaling [19]. Therefore, anoikis-related genes show promise as new targets for tumor therapy.To obtain a reliable ARGS for BLCA, TCGA was used as the training cohort, and the external independent database GEO as the test and validation cohort. To avoid overfitting, we performed a LASSO regression to generate a 15-gene prognostic model that had significant prognostic outcomes in the TCGA and GEO databases. To test the performance of the ARGS, our team built a norman plot combining age, grade, stage, and risk scores. The calibration plots was then used to suggest good discriminative and goodness-of-fit predictions for the prognosis. Zhang’s article built a model for bladder cancer [20]. The areas under the ROC of Zhang’s signature were 0.667, 0.632, and 0.637 for 1, 3, and 5 years, all of which were lower than the areas under the ROC curves of our model. The results demonstrated that our model was more accurate in assessing the survival ratio of patients with 1, 3, and 5 years.The ARGS model we constructed was significantly associated with the survival outcome of BLCA. ARGS contained 15 genes, including CAV1, ITGA5, CSPG4, MYC, NTF3, ZEB2, CRYAB, ZEB1, F10, SFRP1, TAGLN, CCDC80, ID2, RAC3, S100A7. All genes included in the ARGS model have been reported to be associated with tumors. CAV1 is involved in CD26-mediated angiogenesis and metastasis in colorectal cancer [21]. ITGA5, a member of the integrin family, facilitates communication among several cells and cells with ECM. Especially ITGA5 can promote the development of colorectal cancer through its own O-GlcNAcylation [22]. CSPG4 has been investigated to regulate tumor growth and invasion via RTK signaling, MAPK resistance, and integrin signaling by promoting focal adhesion kinase activation [23]. MYC is involved in the CAF oncogenic process in the tumor microenvironment and contributes to tumor proliferation by hindering the PI3K/Akt/β-Catenin signaling [24]. NTF3 as a new target promotes apoptosis in hepatocellular cancer through the signaling of P38 MAPK and JNK signaling [25]. In lung carcinoma, the PAX6-ZEB2 pathway increases cisplatin resistance and migration via PI3K/AKT signaling [26]. CRYAB promotes prostate cancer progression and metastasis as a downstream target of MiRNA-671-5p targeting NFIA [27]. ZEB1 can promote the warburg effect through the activation of PFKM in HCC [28]. Hypermethylation of the F10 promoter in gliomas is an independent prognostic factor and further increases tumor aggressiveness [29]. SFRP1 can affect the proliferation of tumor with Wnt pathway [30]. In colorectal cancer, TAGLN mediates ovarian cancer progression by regulating stromal stiffness through the RhoA/ROCK pathway [31]. CCDC80 plays an important role in melanoma migration via the FAK/CCDC80/E-cadherin axis [32]. ID2 has been investigated to accelerate the aggressiveness of head and neck squamous cell carcinoma via the ID2-MMP1 axis [33]. RAC3 leads to colorectal cancer (CRC) progression and stemming of cancer cells by acting on its downstream target CTFR [34]. S100A7 facilitates esophageal squamous cell carcinoma (ESCC) growth and migration through the incorporation of intracellular JAB1 and increases tumor angiogenesis through the promotion of p-FAK and p-ErK signaling [35]. To explore the biological signaling pathways of ARGS, we performed a GSEA analysis and identified multiple signaling pathways that were enriched in the high-risk group, including hedgehog, MAPK, p53, TGF β, VEGF, and WNT signaling. The above regulation paths have been reported to contribute to BLCA growth. Propofol suppresses the development and stem cellularity of BLCA by inhibiting hedgehog signaling [36]. Mir-186-5p acts on Ras-related protein RAB27 A/B to inhibit BLCA development and EMT process through PI3K/MAPK signaling pathway inactivation [37]. As a cancer suppressor gene, UNC5B-mediated dephosphorylation of p53 in bladder cancer can arrest the cell cycle and inhibit tumor proliferation, providing a therapeutic target for BLCA patients [38]. LncRNA PLAC2 can increase miR-663 which inhibits TGF-β1 to reduce BLCA invasion and metastasis [39]. Silencing RIPK4 will significantly reduce the expression level of VEGF-A to suppress tumor development [40]. TMEM88 can inhibit proliferation and invasion outcomes in BLCA via reducing Wnt/β-catenin signaling [41].Immunotherapy has potency as a new option for BLCA treatment [42]. TME contains a wide range of immune and non-immune components, creating a chronic inflammatory, pro-tumor environment, and immunosuppression [43]. Meng’s article also reported that stromal cells had a great effect on the modification of TME. Stromal cells can suppress the immune process by increasing the secretion of the immunosuppressive factor IL10 to create an immunosuppressive microenvironment [44]. In the risk model, IL10 was substantially enriched in the high-risk ARGS group, indicating that the high-risk ARGS group might have a suppressive immune microenvironment. TAM can shape TME and affect the efficacy of immunotherapy [45]. Higher TAM can predict unfavorable survival outcomes and evaluate the sensitivity of PD-L1 immunotherapy in patients with BLCA [46, 47]. CAF, as the main component of TME, increases BLCA aggressiveness through paracrine IL-6 signaling [48]. Furthermore, CAF can also induce the formation of an immunosuppressive tumor microenvironment [49]. Meng’s article illustrated that the immune-depleted subgroup had a fibroblast-regulated cellular matrix enriched with WNT/TGF-β and LAG3 checkpoint [44]. Our TIDE, GSEA, and checkpoint analysis outcomes also suggested that the high-risk ARGS cluster was concentrated fibroblasts, WNT/TGF-β, and a high level of LAG3 according to the characteristics of the immune-depleted subgroup. The above findings clarified that the high-risk group was an immunosuppressed state.We also explored whether ARGS could identify the sensitivity of patients with BLCA to immunotherapy. The immune checkpoint (ICI) is one of the reliable indicators for evaluating immunotherapy. Our results showed that PD-L1, TIGIT, CD276, and CTLA4 had a high expression level in the high-risk cluster. Additionally, the high ARGS risk group also had a stronger immune exclusion. The TCIA results also showed a better response to ICIs in the low-risk cluster. The above results suggested that the high ARGS risk cluster would not benefit from immunotherapy. We also validated our findings in the immune cohort IMvigor210 with the same results. Overall, we have identified a subtype of BLCA with a better prognosis and sensitivity to immunotherapy. These results will facilitate the improvement of clinical outcomes for BLCA patients and enable precision therapy.ConclusionsIn general, we constructed and validated an anoikis-associated genetic signature based on 15 anoikis-associated genes as a predictor of immunotherapy with a significant prognosis for BLCA patients. Furthermore, ARGS is involved in the regulation of TME and contributes to the development of the immunosuppressive microenvironment. Last, we also explored the combined role of ARGS scores and immune checkpoint genes in patients with BLCA. Based on the above findings, ​we believe that the combination of ARGS scores with specific immune checkpoints has the potential to be a predictive indicator of ICIs response, leading to more precise treatment and increases patient benefit. Therefore, identifying the anoikis-associated gene signature can help us risk stratify and achieve precision therapy, providing new insights to improve the response of BLCA to immunotherapy.

Journal

ONCOLOGIE – de Gruyter

Published: May 1, 2023

Keywords: anoikis; bladder carcinoma; gene signature; immunotherapy; risk model

There are no references for this article.