Hub gene associated with prognosis in bladder cancer is a novel therapeutic target

Data selection

On 23 October 2022, gene expression profiles (including 412 tumor tissues and 19 normal tissues) and copy number variation (CNV) dataset (including 417 cases of bladder cancer samples) with complete follow-up information of bladder cancer patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Besides, the GSE13507 cohort from the GEO database (https://www.ncbi.nlm.nih.gov/) was used to validate the prognostic of XPO1.

Analysis of copy number variation data

The CNV data from 417 bladder cancer samples were annotated using the Genome Research Consortium Human build 38 (GRCh38) as a reference genome. The whole-genome amplification and deletion in the sample were identified using the pipeline Genomic ldentification of Significant Targets in Cancer (GISTIC2.0) (version 2.0.23) (Mermel et al., 2011) from GenePattern (Reich et al., 2006) (https://cloud.genepattern.org/). The levels derived from GISTIC2.0 indicate the copy-number level per gene: −2 or Deep Deletion indicates a deep loss, -1 or Shallow Deletion indicates a shallow loss, 0 is diploid, 1 or Gain indicates a low-level gain, 2 or Amplification indicate a high-level amplification (Cerami et al., 2012; Gao et al., 2013).

Identification of differentially expressed genes (DEGs)

The FPKM data of tumor and normal tissues from TCGA-BLCA were analyzed using the “limma” R-package to obtain DEGs. The threshold for DEGs was as follows: adj P-value (FDR)<0.01 and —log2-fold change (FC)—≥2.

Univariate Cox regression analysis

A total of 109 chemotherapy samples were selected from 412 bladder tumor samples obtained from the TCGA-BLCA, and univariate Cox regression analysis was employed to explore the performance of DEGs in predicting overall survival (OS). Genes were determined as potential prognostic genes with p < 0.05. When HR ≤ 1 indicates that the gene may promote survival, HR ≥ -1 may have an adverse effect on survival. We identified genes with HR > 1 & p < 0.05 as risk prognostic genes.

Construction and evaluation of prognostic model

In this study, seven risk prognostic related genes with copy number amplification and up-regulation in tumors were selected. By integrating expression profile and clinical information, we selected 406 bladder cancer samples with complete survival data, and hub genes were included in the stepwise multivariate COX regression analysis, and the riskscore was constructed according to the standardized regression coefficient of the factors: Riskscore = Σi(Coefi × Expi). Expi represents gene expression of prognostic genes, Coefi represents the gene regression coefficient in multivariate Cox regression analysis.

According to the prognostic model, we calculated the riskscores of bladder tumor samples in the TCGA database, and divided the samples into high-risk group and low-risk group based on the median value of riskscores. t-SNE was drawn by R-package “Rtsne” and the high-dimensional data combined with samples survival information and riskscore was converted into low-dimensional data and visualized to distinguish between high- and low-risk groups. By Kaplan–Meier (K–M) survival curve we identified OS between high and low risk groups. In addition, time-dependent receiver operating characteristic (t-ROC) curve analyses of patients were plotted by R packet “pROC” to assess the sensitivity of the model in predicting half-year, one-year and two-year survival of patients (Robin et al., 2011).

Functional enrichment analysis

The Wilcoxon test was used to detect DEGs between the high- and low-risk groups. DEGs with —logFC—≥1 and FDR<0.05 were used for subsequent analysis. The R “clusterProfiler”(Yu et al., 2012) was used for the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment (Kanehisa, 2019; Kanehisa et al., 2021; Wixon & Kell, 2000) and Gene Ontology (GO) enrichment analysis to identify DEGs-enriched pathways.

Analysis of tumor immune microenvironment

R-package “GSVA” was used to calculate the immune score in the tumor microenvironment. We also explored the infiltration of 13 common immune cells including aDCs, B cells, CD8+ T cells, DCs, iDCs, macrophages, mast cells, neutrophils, NK cells, pDCs, T helper cells, Tfh, Th1 cells, Th2 cells, TIL, Treg. Besides, R-package “estimate” was used to calculate the immune microenvironment scores and tumor cell purity (Yoshihara et al., 2013).

Immunohistochemistry

Samples of bladder cancer diagnosed at the Department of Pathology of The Second Affiliated Hospital of Nanchang University between September 2021 and August 2022 (the project was approved by IBR EC opinion letter of The Second Affiliated Hospital of Nanchang University) were collected and made into tissue chips. Clinicopathological data and written informed consent of patients were obtained. Automatic slicing machine was used for slicing and XPO1 (Proteintech, Rosemont, IL, USA) was dyed by automatic dyeing machine (Roche, Basel, Switzerland).

Cell proliferation assays

The MB49, T24, and HT1376 cells were plated in a 24-well plate at 10,000 cells/well and then they were treated with Selinexor (BiochemPartner, Shanghai, China) or vehicle after cultured in an incubator for 24 h. Proliferation was assessed 24 h, 48 h, 72 h later using Cell Counting Kit-8 (CCK-8) (PointBio, Indianapolis, IN, USA). After the addition of CCK8, the cells were incubated in the incubator for 1–4 h, and then the absorbance of 450nm was measured with an enzyme marker. All conditions were performed in triplicate and repeated at least twice. Data is displayed as the mean ± standard deviation.

Animals

All animals experiments were conducted as approved by Nanchang Royo Biotech company Laboratory Animal Welfare Ethics Committee (RYE2021072501). Six-to-eight weeks old balb/c mice were obtained from Nanchang Royo Biotech company. Mice were raised in independent ventilation cage system, and in different cages according to male and female, with 3–5 mice per cage, and the lighting cycle was day and night every day. Given enough water and chow, mice could obtain it freely. The ambient temperature is maintained at 22−25 °C and the humidity is maintained at 40–60%. Executing animal method is Euthanized with CO2 inhalation. The Experimental end point was 21 days after drug intervention or mouse tumor volume >1500 mm³ or necrosis of tumor ulcer in mice or mice lost more than 20% of the body weight of normal animals (the amount of tumor should be taken into account).

Drug sensitivity test in vivo

Suspensions of MB49 cells (1.0*10 $\stackrel{ˆ}{}$6/200 uL) were subcutaneously injected to balb/c mice establish xenograft tumors. And when the MB49 xenograft tumors reached 50–100 mm³ size, the mice were randomly divided into two groups (n = 6), and then were treated with Selinexor 0 mg/kg/d (Control group) or 20 mg/kg (P.O; tiw/3; treatment). All the mice were treated for 2 weeks. The body weight and tumor volumes were measured every 3 days, and tumor volumes were calculated as (length × width2)/2. TGI was calculated as (1−T/C) × 100%, which T is the average tumor weight of the treatment group and C is the average tumor weight of the control group. The effect of the drug is considered as valid when TGI ≥40% and p < 0.05. Data were analyzed by two way ANOVA.

The prognostic genes with CNV and expression differences were screened

FPKM expression profiles in TCGA-BLCA were selected as the basis for differential expression analysis. Under the condition of FDR <0.01 and —log2FC— ≥2 we identified 2,996 DEGs between bladder tumor and normal samples, including 1,857 up-regulated genes and 1,139 down-regulated genes (Fig. 1A). Through univariate Cox regression analysis, we found 80 DEGs related to the risk prognosis of bladder tumor among 1,857 up-regulated DEGs. In addition, we used GISTIC 2.0 to identify genome amplification/deletion in bladder tumor samples (Figs. 1B and 1C). Venn diagram analysis showed that seven hub genes of the 80 up-regulated DEGs (Fig. 1D) were amplified with copy number (Fig. 2A): ACLY, CNP, NKIRAS2, P3H4, PDIA6, VPS25, XPO1.

A total of 400 bladder tumor samples with survival data were included in TCGA database. Based on these data, the software X-tile (version 3.6.1) (Camp, Dolled-Filhart & Rimm, 2004) determined the best cut-off value of these seven hub genes and divided the samples into high-expression and low-expression groups. K–M curves of the seven hub genes all confirmed that the survival of the high expression group was worse than that of the low expression group (P < 0.05) (Fig. 2B). In addition, in order to explore the relationship between these seven hub genes, we mapped the correlation network containing these seven hub genes and found that their expression was positively correlated with each other (Fig. 2C). Heatmap of seven hub genes amplified and upregulated DEGs in tumor tissue are shown in Fig. 2D.

Construction and evaluation of a prognostic model for bladder tumor

A total of seven hub genes were included in multivariate COX regression analysis (Fig. 3A), and a prognostic model of bladder tumor was established according to the regression coefficients: Riskscore = (0.014 *expression of ACLY) +(0.007 * expression of CNP) +(0.010 * expression of NKIRAS2) +(0.024 * expression of P3H4) +(0.002 * expression of PDIA6) +(−0.004* expression of VPS25) +(−0.010* expression of XPO1). After calculating the riskscore of all tumor samples, samples were divided into high- and low- risk groups with the median riskscore of 0.916. T-SNE analysis showed that high- and low-risk groups could be well distinguished (Fig. 3B). Riskplot shows that riskscore of the sample may be related to prognosis: as the patient’s riskscore increases, the patient’s survival time may decrease (Fig. 3C). The K–M curve further proved that the OS of the high-risk group was significantly worse than that of the low-risk group (P < 0.001) (Fig. 3D). The half-year of area under the t-ROC curve (AUC) was 0.663, 0.676 in1-year, 0.619 in 2-year (Fig. 3E).

Gene enrichment analysis

By applying the Wilcoxon test, 219 DEGs associated with high- and low- risk groups were identified (FDR<0.05). Including 116 up-regulation genes and 103 down-regulation genes. Subsequently, KEGG and GO enrichment analysis showed that most of the high risk-related pathways were PI3K-Akt signaling pathway, laminin binding and other pathways closely related to tumor occurrence, development and metastasis (Figs. S1A–S1B).

Analysis of the immune microenvironment

We calculated the immune score of tumor tissue using R-packet “GSVA”. After evaluating the infiltration of 13 subsets of immune cells in bladder cancer, we found that the mean proportion of T helper cells in all samples was the highest, followed by Treg (Fig. S2A). In the correlations between the expression of seven prognostic related genes and immune cells we found that that expression of XPO1 was negatively correlated to mast cells (cor = −0.31, p = 2.6e − 10) (Fig. S2B). We further explored the correlation between the immune microenvironment, tumor cell purity and seven genes using R-packet “estimate” (Yoshihara et al., 2013). And found that except for VPS25, the other six genes were correlated with the tumor immune microenvironment. In addition, ACLY, NKIRAS2, XPO1 and tumor purity were significantly positively correlated (p < 0.001), which further indicates that ACLY, NKIRAS2 and XPO1 may be potential targets for the treatment of bladder cancer.

Immunohistochemical validation of XPO1 as a hub gene

In order to verify the relationship between seven hub genes and the prognosis of clinical patients with bladder cancer, we selected XPO1 gene as the target gene to detect the expression in pathological samples of patients with bladder cancer. XPO1 is the target of Selinexor, an antitumor drug currently on the market. Immunohistochemical analysis showed that although a small amount of immune response could be observed in the cytoplasm, XPO1 was mainly expressed in the nucleus. The expression of XPO1 in bladder cancer tissues was significantly higher than that in normal tissues. Moreover, the expression level of XPO1 in advanced stages of bladder cancer was higher than that in lower stages (Fig. 4). Besides, the prognostic value of XPO1 was validated in 165 bladder cancer samples from the GSE13507 cohort of the GEO database. The best cut-off value was used to classify high and low XPO1 expression values, and the K–M curve showed that patients with high XPO1 expression had shorter survival in the external bladder cancer cohort.

Selinexor inhibited the proliferation of bladder cancer cell lines

To further verify the inhibitory effect of XPO1 on bladder cancer, we examined the inhibitory effect of Selinexor, an inhibitor of XPO1, on murine bladder cancer cell lines MB49 and human bladder cancer cell lines T24 and HT1376. Experimental results showed that positive controls (Cisplatin and Doxorubicin) were inhibited the bladder cancer lines (Figs. 5A–5B), and the IC50 of Selinexor against MB49 cells was 1.483 (in 24 h), 0.8886 (in 48 h), and 0.4408 (in 72 h) (Fig. 5C). The maximum inhibition rate of Selinexor on T24 and HT1376 cell lines was less than 50% in 24 h. IC50 of Selinexor on T24 cells was 0.8553 (in 48 h), and 0.2467 (in 72 h) (Fig. 5D); The IC50 of Selinexor on HT1376 cells was 0.2702 (in 48 h), and 0.2138 (in 72 h) (Fig. 5E). These results indicated that Selinexor could significantly inhibit the growth of MB49 and T24 bladder cancer cells.

Selinexor inhibited the growth of bladder cancer grafts in vivo

Next, we tested the in vivo effects of Selinexor on bladder cancer. We successfully constructed MB49 bladder cancer grafts. The mice of MB49 bladder cancer grafts were treated by Selinexor for 2 weeks. The results showed that the growth of bladder cancer grafts was significantly inhibited in the twelfth, fifteenth and eighteenth days of Selinexor treatment (Fig. 6). The TGI in Selinexor was 40.5% (much greater than 30%). That confirmed Selinexor significantly inhibited the bladder cancer grafts in vivo.