Expression landscape of VEs-associated genes in PCa

We obtained 27 genes associated with VEs from the CTD database. (Supplement material 1 ). The levels of differences between tumor and normal tissues were further explored. The expression of B3GAT1, CHEK2, DNMT3B, MAP2K1, and RNF43 was markedly upregulated in tumor, whereas the expression of BCL2, CDKN1B, ITGB8, MBNL1, MECON, PAX9, PKNOX2, POU2F2, PPFIBP2, RNASEL, TCF4, TCF7L2, and TTC28 was markedly downregulated in tumor tissues (Fig. 1A). Furthermore, we investigated copy number changes (CNVs) in these genes and reported that CNVs were prevalent in 29 VEs-associated genes. Among them, KDM2A, PAX9, and MECOM presented extensive CNV gains, whereas RNLS, CHD3, CDKN1B, TSHZ1, BCL2, CDH1, TCF4, TTC28, PKNOX2, POU2F2, and CHEK2 presented CNV loss (Fig. 1B). The gene networks describe a comprehensive picture of VEs-associated gene interactions and their impact on the progression-free interval (PFI) in PCa (Fig. 1C). We subsequently performed an enrichment analysis of the effects of VEs on prostate cancer via the CTD database. Among the VEs-enriched diseases, the main correlation was prostate-related diseases, especially prostate cancer (Fig. 1D). In Fig. 1D, it was mainly shown that VEs mainly affect biological processes in GO enrichment analysis. In addition, VEs affected pathways related to the PI3K/AKT signaling pathway, BH3 proteins/BCL-2 members, and EGFR tyrosine kinase inhibitor resistance, as shown by pathway enrichment analysis (Fig. 1D).

Fig. 1

The VEs-associated gene landscape in prostate cancer. (A) Comparison of VEs-associated genes between the tumor and normal groups. (B) Frequencies of CNV gain and loss among VEs-associated genes. (C) Expression and interaction of VEs-associated genes in prostate cancer patients. The size of each cell represents the impact of each gene on the patient’s PFI. (D) Disease, GO, and pathway enrichment analyses of genes associated with VEs

Screening of genes for constructing the VEs score model

Univariate Cox regression analysis was used to further investigate the associations between VEs-associated genes and the PFI. (Supplemental material 2 and Fig. 2A). Subsequently, LASSO analysis was used to construct the PFI model (Fig. 2B). Finally, DNMT3B, KDM2A, PAX9, and TSHZ1 were selected for PFI model construction (Supplemental material 3 and Fig. 2C). The sensitivity of the VEs score was assessed by comparing the concordance indices (C-index) values. Obviously, the C-index of the VEs scoring models were greater than those of the individual genes. (Fig. 2D). Kaplan-Meier analysis manifested that patients with high DNMT3B expression had relatively poor PFIs, whereas the opposite was observed for PAX9. However, the expression of KDM2A and TSHZ1 was not significantly associated with PFI ( Supplement material 4A ). GSEA indicated that high expression of DNMT3B was connected to E2F transcription factors, the G2/M checkpoint and mitotic spindle assembly (Fig. 2E). Moreover, high expression of PAX9 was connected to early response factors, the interferon gamma response and TNFα signaling (Fig. 2E). We subsequently performed protein‒protein interaction analysis for DNMT3B and PAX9 via GeneMANIA, aiming to explore the potential interactions between proteins related to these three proteins (Supplemental material and Fig. 2F). The crystal structure of the DNMT3B-DNMT3L complex was obtained from the Protein Data Bank.

Fig. 2

Screening VEs score-related model genes. (A) Univariate Cox regression analysis of VEs-associated genes. (B) LASSO regression analysis of prognosis-related genes. (C) Multivariate Cox regression analysis of VEs-associated genes. (D) C-index curve of the PFI. (E) Kaplan‒Meier analysis of the PFI of prostate cancer patients with low and high expression of DNMT3B and PAX9 was performed with data from the TCGA database. (F) GSEA of DNMT3B and PAX9 in different expression groups. (G) Determination of genes that interact with DNMT3B and PAX9 in tumors via physical interactions. Crystal structure of the DNMT3B-DNMT3L complex

The VEs score model as a predictor of the PFI

Patients with high VEs scores had significantly worse PFIs in the TCGA database (Fig. 3A). In the two validation cohorts, we reached the same conclusion (Fig. 3B). The evaluation of clinical information for TCGA prostate cancer patients revealed that the VEs score was significantly different from the Gleason score, T stage, and N stage (Fig. 3C and D). We observed an increase in the VEs score with increasing Gleason score (Fig. 3C). GSVA revealed that in addition to E2F-related pathways and the G2/M checkpoint pathway, other signaling pathways were observably activated in the low VEs score group (Fig. 3E). However, according to the results of our analysis, patients with high VEs scores had poor PFIs. Moreover, the cell cycle signaling pathway was activated in the high VEs score group (Supplement material 4B). Androgen receptor (AR) pathway is the core factor in the pathogenesis of PCa. Nevertheless, androgen response was worse in the high VEs score group (Fig. 3E and Supplement material 4B). Subsequently, we explored the associations of AR signaling pathway-related genes with VEs scores, which generally showed a negative trend with AR (Fig. 3F). These results indicate that the high VEs score group has greater malignant potential. Furthermore, there was no difference in AR mutation among the different VEs score groups (Fig. 3I). We subsequently observed a highly positive correlation between TMB and VEs scores (Fig. 3G). The top 10 mutation genes in the low VEs score group included SPOP, TTN, FOXA1, KMT2D, MUC16, TP53, SYNE1, ATM, RP1, and RYR2. The top 10 mutation genes in the high VEs score group included TP53, TTN, SPOP, KMT2D, MUC16, FOXA1, KMT2C, LRP1B, PTEN, and SPTA1 (Fig. 3H).

Fig. 3

VEs score model and associations with survival and clinical features. (A-B) Kaplan‒Meier analysis of the VEs score low and high groups in the TCGA, GSE116918 and GSE94767 cohorts. (C) VEs scores of different Gleason score groups in TCGA. (D) Differences in clinical characteristics between the low- and high-VEs score groups. (E) GSVA revealed pathways enriched with different VEs scores. (F) Pearson’s correlation analysis between the VEs score and the expression of AR signaling pathway-related genes. (G) Relationship between VEs score and TMB. (H) Top 10 mutated genes in the two VEs score groups. (I) Mutation comparison of the AR gene in different VEs score groups

Predictive potential of the VEs score in the response to cancer immunotherapy

To reveal the predictive effect of the VEs score on cancer immunotherapy efficacy, we downloaded and evaluated the activity score of the cancer immune cycle from the Tracking Tumor Immunophenotype database. The VEs score was positively correlated with the release of cancer cell antigens, priming and activation, T-cell recruitment, eosinophil recruitment, basophil cell recruitment, B-cell recruitment, Th2 cell recruitment, and Treg cell recruitment. However, Th17 cell recruitment, MDSC recruitment and recognition of cancer cells by T cells were inversely associated with the VEs score (Fig. 4A and Supplementary material 6 ). In addition to immune cells, immune-related regulatory factors also play a key role in immunotherapy [19]. Except for CTLA4, there was no difference in the expression of other genes between the two groups. (Fig. 4B). Meanwhile, we investigated the association between the VEs score and major histocompatibility complex (MHC) and found that TNFSF14, TNFSF18, TNFSF25, and CD28 were increased in the high VEs score group. (Fig. 4C). In terms of the expression of CTLA4, the high VEs score group was more sensitive to immunotherapy (Fig. 4D). Moreover, the cytolytic activity (CYT) presented a positive correlation with VEs score (Fig. 4E). GSEA demonstrated that immune-associated pathways were upregulated in the high VEs score groups (Fig. 4F). Therefore, we believe that patients with high VEs scores have a more active immune status and are more suitable for combination immunotherapy.

Fig. 4

Estimation of anticancer immune activity via the VEs score. (A) Heatmap of the seven-step cancer immunity cycle and VEs score. (B-C) The expression of immune checkpoint-related genes and the MHC gene set in the two VEs score groups. (D) IPS between the two VEs score groups. (E) Relationships between VEs scores and CYT scores. (F) GSEA for the high VEs score group

Effect of the VEs-associated gene DNMTB3 on prostate cancer progression

To explore the function of DNMT3B in PCa, two small interfering RNAs (siRNAs) targeting DNMT3B were transfected into PCa cells. The mRNA level of DNMT3B was markedly knocked down (Fig. 5A). Subsequently, CCK8 and colony formation assays were used to verify the interference of DNMT3B knockdown on proliferation. The proliferation ability of tumor cells was significantly inhibited when DNMT3B was knocked down (Fig. 5B and C). A cell migration assay revealed that decreased DNMT3B expression resulted in decreased migration (Fig. 5D). Therefore, the above results suggest that DNMT3B has a critical role in prostate cancer progression. Moreover, we observed that the expression level of DNMT3B was increased after VEs treatment ( Supplement material 7).

Fig. 5

DNMT3B promotes the proliferation of prostate cancer cells in vitro. (A) qRT‒PCR analysis of DNMT3B expression levels in DNMT3B-knockdown and control cells. (B) Cell viability was assessed in DNMT3B-knockdown prostate cancer cells. (C) Colony formation assays were performed in DNMT3B‐knockdown prostate cancer cells. (D) Representative images of prostate cancer cells migration assays showing cell migration after downregulation of DNMT3B