Identification of key Enz-resistant genes in PCa based on RNA-seq data

To identify potential target genes linked to resistance to Enz in prostate cancer, a search was carried out in the GEO database to identify genes resistant to Enz for further comprehensive analysis. A comparison of gene expression levels was conducted between LNCaP cells treated with Enz and the control group from the GSE44905 dataset, resulting in the identification of 954 upregulated DEGs and 728 downregulated DEGs that met the criteria of log2FC > 1.3 and p < 0.05 (Fig. 1A). The study conducted a comparative analysis of gene expression levels in Enz-resistant LNCaP cells and Control LNCaP cells using data from the GSE104935 dataset, identifying 1121 upregulated and 859 downregulated DEGs meeting the criteria of log2 FC greater than 1.5 and p-value less than 0.05 (Fig. 1B). Furthermore, a comparative analysis was conducted on VCaP cells treated with Enz and the control group from the GSE51872 dataset, identifying a total of 962 upregulated and 771 downregulated DEGs that met the specified criteria of a log2 FC greater than 4 and a p-value less than 0.01 (Fig. 1C). The heatmaps visually represented the clustering of these DEGs (Fig. 1D–F).

Fig. 1

Enz resistance-related genes identification. A Volcano plots of the DEGs in GSE44905. Red dots indicate upregulated genes; blue dots indicate downregulated genes. B Volcano plots of the DEGs in GSE104935. Red dots indicate upregulated genes; blue dots indicate downregulated genes. C Volcano plots of the DEGs in GSE51872. Red dots indicate upregulated genes; blue dots indicate downregulated genes. D Heatmap showing DEGs in different samples in GSE44905. E Heatmap showing DEGs in different samples in GSE104935. F Heatmap showing DEGs in different samples in GSE51872. G Circular enrichment of GO pathways among GSE44905 DEGs. H Circular enrichment of GO pathways among GSE104935 DEGs. I Circular enrichment of GO pathways among GSE51872 DEGs. J Circular enrichment of KEGG pathways among GSE44905 DEGs. K Circular enrichment of KEGG pathways among GSE104935 DEGs. L Circular enrichment of KEGG pathways among GSE51872 DEGs

To clarify the biological characteristics linked to the emergence of Enz resistance, distinct analyses were conducted on the overexpressed DEGs in the cohorts treated with Enz and those displaying resistance to Enz, employing GO and KEGG pathway analysis. The predominant GO terms related to the reaction to Enz are depicted in Fig. 1G–I. DEGs in LNCaP and VCaP cells subjected to Enz treatment (GSE44905, GSE51872), as well as in Enz-resistant LNCaP cells (GSE104925), exhibited enrichment in pathways linked to cellular homeostasis, mitotic DNA replication, lipid oxidation, and regulation of apoptotic signaling pathways, suggesting a potential role in regulating drug resistance in PCa. Additionally, the primary KEGG terms associated with Enz response are depicted in Fig. 1J–L. The DEGs were mainly associated with key signaling pathways like PI3K–Akt, NF-kappa B, and P53, consistent with GO enrichment analysis. This emphasizes the importance of studying molecular mechanisms to improve understanding and treatment of Enz-resistant PCa.

Enz-resistant gene co-expression networks determined by WGCNA

DEGs were primarily focused on the most significantly regulated genes, potentially neglecting others. This study utilized WGCNA to construct gene co-expression modules to investigate key gene modules and potential mechanisms associated with Enz resistance. Hierarchical clustering was employed to create color-coded modules for the cluster dendrogram, with heatmaps presented in Fig. 2A (GSE44905), Fig. 2B (GSE104935), and Fig. 2C (GSE51872). The blue module exhibited the strongest positive correlation with the Enz response trait in GSE44905 (corresponding correlation, CC = 0.92, p = 4e−04), as depicted in Fig. 2D. Additionally, the magenta module demonstrated a significant positive correlation with the Enz resistance trait in GSE104935 (CC = 0.87, p = 0.001), while another green module exhibited the strongest positive correlation with the Enz response trait in GSE51872 (CC = 0.98, p = 4e−06), as depicted in Fig. 2E, F.

Fig. 2

WGCNA reveals Enz-resistant gene co-expression networks. A WGCNA analysis of Enz response samples in GSE44905. The dendrogram represented the clusters of differentially expressed genes based on different metrics. Each branch represented one gene, and each color below the branches represented one co-expression module. B WGCNA analysis of Enz resistance samples in GSE104935. The dendrogram represented the clusters of differentially expressed genes based on different metrics. Each branch represented one gene, and each color below the branches represented one co-expression module. C WGCNA analysis of Enz resistance samples in GSE51872. The dendrogram represented the clusters of differentially expressed genes based on different metrics. Each branch represented one gene, and each color below the branches represented one co-expression module. D The heatmap showed the correlation between gene modules and Enz response. The correlation coefficient in each cube represented the correlation between gene modules and traits, which decreased from red to blue. E The heatmap showed the correlation between gene modules and Enz resistance. The correlation coefficient in each cube represented the correlation between gene modules and traits, which decreased from red to blue. F The heatmap showed the correlation between gene modules and Enz resistance. The correlation coefficient in each cube represented the correlation between gene modules and traits, which decreased from red to blue. G The common hub genes shared between DEGs and WGCNA derived from GSE44905 were visualized in a Venn diagram. H The common hub genes shared between DEGs and WGCNA derived from GSE104935 were visualized in a Venn diagram. I The common hub genes shared between DEGs and WGCNA derived from GSE51872 were visualized in a Venn diagram. J The top enriched GO pathways among common hub genes from the GSE44905. The horizontal axis represented the p-value of GO terms on Metascape. K The top enriched GO pathways among common hub genes from the GSE104935. The horizontal axis represented the p-value of GO terms on Metascape. L The top enriched GO pathways among common hub genes from the GSE51872. The horizontal axis represented the p-value of GO terms on Metascape

To elucidate the crucial genes involved in the development of resistance to Enz, an effort was made to identify shared hub genes among DEGs and hub genes derived from WGCNA modules. This investigation resulted in the identification of 79 genes in GSE44905, 42 genes in GSE104935, and 91 genes in GSE51872. Subsequently, an enrichment analysis was conducted on the common hub genes identified from the three GEO datasets using Metascape (https://metascape.org/). The primary Metascape terms identified in GSE44905 included the regulation of the G protein-coupled receptor signaling pathway and positive regulation of the MAPK cascade. Furthermore, the principal GO enrichment outcomes of GSE104935 comprised the immune response-regulating signaling pathway and lysosomal transport. Additionally, the top terms of GSE51872 identified through GO analysis included response to steroid hormone, PID HIF1 TFPATHWAY, and regulation of chemokine production.

Expression and prognosis of hub genes in PCa

To inform clinical decision-making and stratify Enz-resistant conditions, it is essential that the information obtained from Enz-resistant samples remains consistent. Consequently, a comparative analysis was performed on hub genes identified in three GEO datasets, revealing that CBX2, DCAF6, FAM117B, and TMEM141 were consistently identified across all three datasets (Fig. 3A). The mRNA expression levels of these four candidate genes in TCGA dataset are depicted in Fig. 3B. Consequently, CBX2 demonstrates elevated expression levels in PCa samples and is notably associated with the prognostic outcomes of PCa patients as indicated by Kaplan–Meier survival analysis (Fig. 3C). Moreover, among the four genes identified through Venn diagram analysis, CBX2 displays statistically significant differences between the control and Enz-treated LNCaP groups (Fig. 3D–E), as well as between the control and Enz-resistant VCaP groups (Fig. 3F). Taken together, these findings suggest that CBX2 may serve as a promising biomarker for Enz resistance in PCa.

Fig. 3

Expression and prognosis of hub genes. A The common hub genes shared among three GSE datasets were visualized in a Venn diagram. B The expression of four selected genes in the TCGA database, PCa (red) and normal (gray). C Kaplan–Meier curves, respectively, for disease-free survival (DFS), and of PCa patients with high versus low expression of four selected genes in TCGA. D Estimation of hub gene expression in GSE44905. E Estimation of hub gene expression in GSE104935. F Estimation of hub gene expression in GSE51872. *p < 0.05

CBX2-related downstream signaling pathways identification

After establishing CBX2's involvement in Enz resistance in PCa, our objective was to investigate the molecular mechanisms underlying CBX2's pro-resistance role in PCa through gene expression analysis of scRNA-seq data from a variety of PCa samples. Subsequently, gene expression profiles were extracted from two Enz-treated samples and one control sample within the GSE215943 dataset. Following appropriate sequencing depth adjustment, gene count detection, and data normalization steps, a subset of 2000 genes showing significant variability was chosen for subsequent analysis. Subsequently, after the application of quality control filters, a total of 31,312 cells were analyzed, each containing a median of 24,371 genes per cell.

The dimensionality reduction procedure was executed utilizing the "RunPCA" function, leading to the discernment of 12 clusters at a resolution of 0.5 (Fig. 4A, B). Next, an evaluation of the distribution of these cellular clusters across the two sample types was conducted and depicted (Fig. 4C). The heatmap exhibited the top 5 genes displaying notable variability within each cluster (Fig. 4D). Following this, a comprehensive analysis of the visualization outputs generated by CellChat was undertaken for various analytical objectives. The hierarchy plot and circle plot employed distinct color schemes to depict distinct cell groups, whereas the bubble plot utilized color intensity to indicate the relative likelihood of communication (Fig. 4E, F). The results suggest that clusters 7 and 10 primarily functioned as sources, clusters 3, 5, and 6 as recipients, and clusters 0, 1, 2, 4, 8, 9, and 11 as both sources and recipients. Additionally, an investigation of the interactions among the nine clusters entailed examining interactions in which each cluster served as both a target and a source. The bidirectional interactions between outgoing and incoming signaling within these clusters involved pathways such as NOTCH, WNT, and IGF (Fig. 4G).

Fig. 4

ScRNA-seq data validated the selected hub genes. A The UMAP was utilized to partition cells into 12 distinct clusters, with each cluster being represented by a unique color that corresponded to its numbered phenotype. B The UMAP was used to partition cells according to different sample types. C The distribution of cell proportions across distinct groups is depicted on a graph, where the vertical axis denotes the proportion of each cell cluster. D The heatmap presents the expression of the top 5 DEGs (rows) in each cell cluster (column). E Cell–cell communication signaling network among the 9 clusters analyzed with CellChat. The right panel showed that cell clusters were located based on the count of their significant incoming (Y-axis) or outgoing (X-axis) signaling pattern. F Number of interactions in 12 cell clusters. The width of the lines indicates the number of pairs. Different colors represent different signal sources. G Heatmap of the CellChat signaling in each cluster. The left panel shows the outgoing signaling patterns (expression weight value of signaling molecules) and the right panel shows the incoming signaling patterns (expression weight value of signaling receptors). A gradient of white to dark green indicates a low to high expression weight value in the heatmap. H The inferred NOTCH2 signaling pathway network. I Feature plots showing the distribution of CBX2, DCAF6, FAM117B, and TMEM141 in 12 cell clusters. J Violin plots showing the distribution of indicated genes in various cell clusters

Prior studies have demonstrated a strong association between CBX2 and the Notch signaling pathway in kidney renal papillary cell carcinoma [32]. To investigate the role of CBX2 in the development of Enz resistance in PCa, our study examined the Notch signaling pathway in cluster communications. Our findings indicate that clusters 0 and 6 serve as primary sources, while the other clusters act as main receivers. Additionally, the distribution and expression patterns of hub genes CBX2, DCAF6, FAM117B, and TM3M141 are depicted in Fig. 4I, J.

Furthermore, to emphasize the predominant signaling pathway activated within each cluster, GSVA was conducted across the 12 clusters. The results indicated that in clusters 0 and 1, there is an upregulation of the HALLMARK_INTERFERON_ALPHA_RESPONSE. Besides, there is an upregulation of the HALLMARK_ANDROGEN_RESPONSE in cluster 2, and an upregulation of the HALLMARK_PANCREAS_BETA_CELLS in cluster 3, among others (Fig. 5A). Subsequently, the cells were categorized into CBX2+ and CBX2− groups for the purpose of analyzing gene expression disparities. GO enrichment analysis revealed enrichment in cell survival pathways such as the cell cycle process, DNA metabolic process, and mitotic cell cycle. KEGG enrichment analysis revealed the involvement of the P53 signaling pathway, Cell cycle, and Apoptosis (Fig. 5B, C). Additionally, GSEA demonstrated the essential role of CBX2 activity in promoting prostate cancer (PCa) cell survival by upregulating the cell cycle signaling in CBX2+ cells. Moreover, the inhibition of the P53 signaling pathway was associated with CBX2− cells, indicating that the upregulation of CBX2 suppresses the P53 signaling pathway, thereby contributing to enzalutamide resistance in PCa (Fig. 5D, E).

Fig. 5

In vitro experiments validated the role of CBX2 in Enz-response. A Heatmap showing different pathways enriched in the 12 cell clusters by GSVA. Each column represents different groups or subpopulations of cells, and each row represents a pathway. The redder the color, the higher the score, and the bluer the color, the lower the score. B The top enriched GO pathways among DEGs of CBX2+ and CBX2− cells were identified and presented in a graphical representation. C The top enriched KEGG pathways among DEGs of CBX2+ and CBX2− cells were identified and presented in a graphical representation. D The P53 signaling pathway was identified and presented in a graphical representation by GSEA. E The cell cycle pathway was identified and presented in a graphical representation by GSEA

In vitro experiments validated the role of CBX2 in Enz-response

To validate the hypothesis that Enz response can activate CBX2 and enhance prostate cancer cell survival through modulation of the P53 signaling pathway, resulting in Enz resistance, an in vitro model of PCa was established using LNCaP cells. The expression levels of CBX2 and P53 were assessed following treatment of LNCaP cells with Enz at a concentration of 10 μM for a specified duration. Western blot analysis confirmed that stimulation with Enz led to downregulation of the P53 signaling pathway in LNCaP cells (Supplement Fig. A, B).

Our functional analysis indicated that the acquisition of Enz resistance is linked to alterations in the P53 signaling pathways. Collectively, the findings from our colony formation assay suggest that resistance stemming from CBX2 activation could potentially be mitigated through a combination of CBX2 knockdown and Enz treatment (Supplement Fig. 1C, D). Consequently, we posit that the CBX2-related pathway may contribute to the development of Enz resistance. Our collective in vitro experimental results indicate that the inhibition of CBX2 may lead to the reactivation of the P53 pathway and suppression of resistance to Enz in prostate cancer. This underscores the therapeutic potential of targeting CBX2 in the treatment of advanced prostate cancer.