Prediction of Prostate Adenocarcinoma Recurrence Prognosis and Immune Status Through 6-Acetoxy-Anopterine Resistance-Associated Programmed Cell Death Genes

Introduction

Prostate adenocarcinoma (PRAD) is a common cancer among men, with its incidence and mortality rates steadily increasing worldwide.¹ Despite significant advancements in the early diagnosis and treatment of PRAD, recurrence and metastasis remain critical clinical challenges.² Especially for patients with advanced or metastatic prostate cancer, traditional treatment methods such as androgen deprivation therapy and chemotherapy often fail due to tumor cell resistance, leading to poor prognosis.³ Thus, revealing the detail molecular mechanism and key gene associated with drug resistance are critical for the clinical intervention of PRAD.

In recent years, with the development of molecular biology techniques, many researchers focused on programmed cell death (PCD) in tumor cells and its role in anticancer drug resistance, since recent studies have indicated that chemoresistance is facilitated by disrupting or circumventing apoptotic cell death mechanisms.⁴ PCD is an important mechanism for cellular self-regulation, playing a crucial role in maintaining tissue homeostasis and resisting tumor development.⁵ Recent studies discovered that YAP signaling plays an important role in the regulation of PCD in cancer, and it is reported that YAP signaling regulates PCD and drug resistance in cancer.⁶ Also, the abnormal regulation of PCD-related genes (DRGs) or pathways are closely associated with occurrence, development and molecular drug resistance.⁷ It has been proved that PCD contributes to the progression of PRAD via participating in the necroptosis and pyroptosis process.⁸ Li et al have highlighted that PLAG1 interacts with GPX4 to conquer vulnerability to sorafenib induced ferroptosis through a PVT1/miR-195-5p axis-dependent manner in hepatocellular carcinoma.⁹ Wang et al have revealed that inhibition of TPI1 sensitizes cisplatin-resistant oral cancer to ferroptosis.¹⁰ Guo et al developed an immune-based prognostic model utilizing four identified DRGs, which were subsequently chosen as predictors for the robust prognostic model of PRAD.¹¹ Actually, during the biological function of PCD, the DRGs are commonly associated with inflammation-related pathway, which further participating in the immune therapy and drug resistance of cancer.¹² For instance, Liu et al have found that CPLX1 is a novel prognostic biomarker in colorectal cancer correlating with immunotherapy resistance and ferroptosis.¹³ NOXA is a BCL2 family protein, and it is reported that NOXA as a pivotal regulator of resistance to CAR T-cell therapy by impairing apoptosis of tumor cells, and low NOXA expression in tumor samples is correlated with worse survival in a tandem CD19/20 CAR T clinical trial in relapsed/refractory B-cell lymphoma.¹⁴ As a novel anticancer drug, 6-acetoxy-anopterine (6-AA) is a natural product belonging to the organic compound family. Extracted from a plant called anopterus macleayanus, 6-AA shows promising potential as an anti-cancer compound. Studies have revealed its efficacy in inhibiting various cancer cells, suggesting its role as a potent anticancer agent in PRAD.¹⁵ However, its detailed pharmacology and clinical applications are still under research. Although 6-AA is not widely used in clinical currently, the development of drug resistance might greatly limit its clinical efficacy. Therefore, exploring the role of DRGs in PRAD, especially those related to drug resistance of 6-AA, has important clinical significance.

In this study, we screened differentially expressed programmed cell death genes (DEPGs) associated with 6-AA resistance based on the TCGA-PRAD database. Then, a prognostic risk prediction model was established to evaluate the potential role of these genes in the survival outcomes of PRAD. Moreover, the relationship between the risk model and the immune microenvironment of PRAD samples was revealed, followed by the characteristics and mechanisms investigation of immune cell infiltration in different risk groups. Finally, the application prospects of the risk model in predicting drug response sensitivity were explored. We aim to reveal the key roles of 6-AA resistance-related PCD genes in the prognosis and immune status of disease recurrence, thereby providing new insights and targets for clinical treatment strategies for patients with PRAD.

Materials and Methods

Microarray Data and Pre-Processing

RNA-seq data (log2(fpkm+1)) from The Cancer Genome Atlas (TCGA)-PRAD, along with related clinical information, were retrieved from the TCGA database. The samples were included in this study with the following exclusion criteria: (a) Samples with missing or zero recurrence-free survival (RFS) times; (b) Deletion of missing values and unexpressed genes exceeding 50% of the total sequenced number; (c) Samples with over 50% missing values or unexpressed genes; (d) All expression values are log2 logarithmic processed. Finally, a total of 52 normal samples and 490 tumor samples were enrolled in current TCGA dataset. Meanwhile, the microarray datasets GSE54460 (including 106 PRAD samples with complete RFS prognostic information)¹⁶ in Gene Expression Omnibus (GEO) database were enrolled in current study. The platform for GSE54460 was GPL11154 Illumina HiSeq 2000 (Homo sapiens). In addition, GSE81277¹⁵ (platform: GPL16604 Agilent-032034 VPC Human 180K v3) including 3 6-AA-sensitive samples and 4 6-AA-resistant samples were obtained from GEO database.

DEPGs Investigation

The classic Bayesian approach from the limma package (version 3.10.3) in R¹⁷ was employed to identify the DEGs between tumor and normal samples in the TCGA dataset. in R was used to reveal the DEGs in tumor samples vs normal samples in TCGA dataset. The Benjamini & Hochber (BH) adjusted P-value < 0.05 and |log₂FC| > 0.5 were selected as the thresholds for DEGs investigation. The results were visualized by volcano plot. Then, DEGs between 6-AA-resistant samples and 6-AA-sensitive samples in GSE81277, identified using the Wilcoxon test with P < 0.05, were considered as 6-AA-specific genes. Finally, based on a review of previous literature,^18,19 the genes related to PCD were investigated. Based on the online tool Venny, we intersected the PCD genes with DEGs and 6-AA-specific genes to identify DEPGs for further analysis.

Enrichment Analysis on DEPGs

GO function and KEGG pathway analysis were performed on DEPGs using clusterProfiler package (version: 4.0.5) in R.²⁰ The GO functions including biological process (BP), cellular components (CC), and molecular function (MF). The Benjamini & Hochberg method was used for multiple test corrections, and the corrected P-value (adj. P value) was obtained. The adj. P < 0.05 and Count ≥ 2 were used as the cut-off value. The top 10 result of GO and KEGG enrichment was visualized using bar chart.

Prognostic Signature Investigation

The univariate Cox regression investigation was applied on DEPGs genes to reveal prognostic genes using a cut-off value of P < 0.05. Then, we utilized Least absolute shrinkage and selection operator (LASSO) Cox regression within R (version: 3.6.1)²¹ to explore the most advantageous gene collection among prognostic genes within TCGA samples with a 10-fold cross-validation. Then, the following formula was utilized to establish Riskscore (RS) model:

In the formula, β denoted the regression coefficient, h₀(t) represented the baseline hazard rate, and h (t, X) was the hazard rate at time t associated with X (the covariates). Then, the RS for each sample in TCGA and GEO dataset were computed using the risk score calculation formula, and samples were divided into two risk groups (high and low) based on the RS median value for survival analysis comparison. Subsequently, the kaplan-Meier (KM) curve in survival package of R was used to assess the association between risk grouping and the actual survival prognosis information.

The Prognostic Evaluation for Signature

By integrating clinical data of prostate cancer, the association among the RS and clinical information including T stage and Gleason score were investigated. To further explore the prognostic independence of clinical factors and RS, the clinical factors and RS from the PRAD data were analyzed by using both univariate and multivariate Cox regression. Independent prognostic factors were identified using a threshold of P < 0.05, followed by visualized by forest plot. Then, a nomogram was constructed by the value nomoScore of signatures by using rms package (version: 6.3-0)²² in R. We assessed group differences and used ROC analysis to compute the survival AUC.

Immune Infiltration and Correlation Analysis

The CIBERSORT methodology²³ was applied to estimate the levels of infiltration of 22 types of immune cells within the tumor samples with the cutoff value of P < 0.05, and the parameters were set as “perm = 50, QN = TRUE”. To reveal disparities in immune cell compositions between two different risk groups, ESTIMATE analysis was utilized.²⁴ Additionally, we examined the correlation between genes and immune cells, as well as among different cell types, using the Wilcoxon test.

Gene Set Enrichment Analysis (GSEA) Between Subtypes

Enrichment analysis comparing subtypes utilized the h.all.v7.4.symbols.gmt background from MSigDB v7.1.²⁵ Enrichment scores for each Hallmark gene set in periodontitis samples were computed and ranked using the GSVA algorithm in R,²⁶ applying a BH adjusted P (FDR) < 0.05 and Normalized Enrichment Score (NES) >1.

Drug Sensitivity Analysis of 6-AA

The GDSC database was utilized to evaluate the responsiveness to chemotherapy drugs. The difference of sensitivity (IC50 value difference) for 138 chemotherapy drugs were quantify by using the pRRophic algorithm.²⁷ Then, the Spearman correlation coefficient and corresponding P < 0.05 were used to reveal whether the difference in IC50 value of each drug associated with 6-AA. The result was visualized using box plot.

Immunotherapeutic Responses Prediction

TIDE analysis predicts immunotherapeutic responses by assessing two primary mechanisms of tumor immune escape.²⁸ In this study, the Wilcoxon test was employed to uncover differences in TIDE scores between the two risk groups. Moreover, the cytolytic activity (CYT) score is defined as the average expression value of Granzyme A (GZMA) and Perforin 1 (PRF1). The tertiary lymphoid structure (TLS) score is calculated using the GSVA algorithm from the GSVA package in R, based on the expression of TLS signature genes. Subsequently, the Wilcoxon test was used to compare the differences in CYT and TLS scores between different risk groups. Furthermore, the Wilcoxon test was used to revealed the difference of immune checkpoint genes between two risk groups.

Single-Cell Analysis

The GSE176031 single-cell dataset includes precise cell type annotations at the individual cell level, enabling detailed analysis of gene expression across various cell types. In this study, the differential gene expression across various cell types elucidates the heterogeneity of the TME in different PRAD patients, thereby contributing to our understanding of cancer heterogeneity.

Clustering Analysis for DEPGs

To investigate different clusters of PRAD associated with PCD, the ConsensusClusterPlus package²⁹ in R software (version: 4.3.2) was used for analysis by improving the clustering variable (k) from 2 to 9. KM survival analysis was then performed using the Survival package (version 2.41-1) to assess survival among different clusters. Following this, the t-test was employed to examine the relationship between clusters and clinical data, such as age and stage.

The qRT-PCR Analysis

To further investigate the expression of biomarkers (TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1, and LAPTM4B) revealed in the current study, a verification study using qRT-PCR was performed using cell lines and clinical samples, respectively. For cell experiment, the expression of all six genes in cultured LNCaP cells (ATCC^® CRL-1740™, 6-AA sensitive group) and DU145 cells (ATCC^® HTB-81™, 6-AA resistant group) were investigated. In addition, the expression of all six genes was revealed based on the clinical samples from 6-AA resistant (n = 20) and 6-AA sensitive (n = 20) PRAD patients. The clinical samples were obtained from our hospital. All procedures involving human participants were conducted in accordance with the ethical standards of our committee. This study protocol was complied with the Declaration of Helsinki and was approved by the Ethics Committee of the Affiliated Hospital of Shandong Second Medical University. Informed consent was obtained from all individual participants included in the study. Briefly, total RNAs were extracted using TRIZOL reagent (Invitrogen, USA) and reverse transcripted using RevertAidTM First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) in accordance with manufacturers’ instructions. The PCR was performed on ABI7500 (Applied Biosystems, U.S.A.). The detailed information for all primers used in current study were listed in Table 1. The PCR program included 95°C for 5 min, 35 cycles of 95°C for 30 s and 52°C for 30 s. The relative expression was calculated using the 2^−ΔΔCt method.³⁰

Table 1 The Detail Information for All Primers Used in Current Study

Results

DEPGs Investigation

A total of 2585 DEGs were explored between tumor samples and normal samples (Supplementary Figure 1A). Then, based on GSE81277, totally 2198 6-AA-specific genes were further revealed (Supplementary Figure 1B). Finally, the PCD genes with DEGs and 6-AA-specific genes were used to identify DEPGs. The VENN plot analysis showed that there were 57 DEPGs enrolled for subsequent analysis (Supplementary Figure 1C).

The Enrichment Analysis on DEPGs

The GO enrichment analysis performed on DEPGs showed that these genes were predominantly assembled in functions like regulation of cysteine-type endopeptidase activity (BP, GO:2000116, Genes: BCL2L12, MYC, LTF, etc) (Figure 1A), ErbB-3 class receptor binding (MF, GO:0043125, Genes: PIK3R1 and ERBB3) (Figure 1B) and late endosome membrane (CC, GO:0031902, Genes: MCOLN2, LAPTM4B and TMEM106B, etc) (Figure 1C). Meanwhile, these genes were primarily enriched in p53 signaling pathway (hsa04115) (Figure 1D).

Figure 1 The enrichment analysis performed on the differentially expressed 6-acetoxy-anopterine (6-AA) resistant programmed cell death-related genes (DEPGs). (A) The Top 10 GO-BP function that significantly assembled by DEPGs. (B) The Top 10 GO-MF function that significantly assembled by DEPGs. (C) The Top 10 GO-CC function that significantly assembled by DEPGs. (D) The Top 10 GO-BP function that significantly assembled by DEPGs. The X-axis represented the number of genes enriched in certain function or pathway, while the Y-axis represented the name of function or pathway.

Prognostic Signatures Investigation and Validation

Totally 18 prognostic genes were revealed by using univariate Cox regression (Supplementary Figure 2). LASSO regression pinpointed 11 significant genes (Figure 2A and B). Then, multivariate COX regression identified six signature genes including TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B (Figure 2C). Moreover, based on the RS of each patient, we separated the sample into two risk groups. For TCGA training dataset, the expression of all six signature genes was noticeable divergence between these two risk groups (all P < 0.05) (Figure 2D). Meanwhile, there was a significant different in survival probability between two risk groups (Figure 2E). The ROC analysis for signatures showed that the AUC value of 1-year, 3-years and 5-years were larger than 0.717, which showed a well diagnostic value for these signature genes (Figure 2F). The trend of the results in the GEO validation dataset was basically consistent with that in the TCGA training dataset (Figure 2G–I).

Figure 2 The prognostic signature investigation and validation. (A) The least absolute shrinkage and selection operator (LASSO) Cox analysis revealed 11 significant genes. (B) The best penalty coefficient in LASSO model. (C) The forest plot showed six signature genes based on multivariate COX regression. (D) The distribution of genes based on RS, the difference of survival time between two outcomes of disease, and gene expression based on two risk groups in the TCGA training dataset. (E) The KM analysis showed the survival probability in two different risk groups in TCGA training dataset. (F) The ROC analysis showed the AUC value of 1-year, 3-years and 5-years survival for current signatures in TCGA training dataset. (G) The distribution of genes based on RS, the difference of survival time between two outcomes of disease, and gene expression based on two risk groups in the GEO validation dataset. (H) The KM analysis showed the survival probability in two different risk groups in GEO validation dataset. (I) The ROC analysis showed the AUC value of 1-year, 3-years and 5-years survival for current signatures in GEO validation dataset. *P < 0.05; **P < 0.01; ***P < 0.001.

Independent Prognostic Analysis

The univariate Cox regression on clinical characteristics showed that T_Stage, Gleason, Resection_Status and RS were independent prognostic factors for PRAD (Figure 3A). Meanwhile, Gleason and RS were independent prognostic factors for PRAD (Figure 3B) based on multivariate Cox regression. To further reveal the correlation between the two factors (Gleason and RS) and prognosis, these factors were enrolled for nomogram construction based on signature genes (Figure 3C). The numerical rating of each parameter was assigned along a continuum. Summing these individual ratings provided a method to predict survival outcomes, indicating the applicability of the nomogram for assessing PRAD risk. The survival prediction analysis showed that the 1-year, 2-years and 5-years survival rate were accordance with the actual survival rates (Figure 3D), which further indicated the prognostic value of current nomogram. Furthermore, the survival time was significantly different between two risk groups based on nomogram (Figure 3E). In addition, the ROC analysis for nomogram showed that the AUC value of 1-year, 3-years and 5-years were larger than 0.769, which showed a well diagnostic value for current nomogram (Figure 3F).

Figure 3 The interdependent analysis and nomogram investigation. (A) The univariate Cox regression revealed the correlation between the RS and clinical factors. (B) The Cox multivariate regression revealed the correlation between the RS and clinical factors. (C) A nomogram for predictive ability study. (D) The consistency line chart between 1-year, 2-years and 5-years survival rate predictions and actual survival rates: the X-axis represented the predicted survival rate, while the Y-axis represented the actual survival rate. (E) The KM analysis showed the survival probability in two different risk groups based on nomogram. (F) The ROC analysis showed the AUC value of 1-year, 3-years and 5-years survival for current nomogram.

Immune Infiltration and GSEA Analysis

Based on TCGA dataset, the immune score and matrix score were calculated using the estimate algorithm. The results showed significant differences in matrix score and immune score among different risk groups (Figure 4A). Then, the CIBERSORT was used to calculate the proportion of 22 immune cells in TCGA dataset. With P < 0.05, a total of eight immune cells’ expression were revealed dramatically different between two risk groups (Figure 4B). For example, the expression of CD4+ resting memory T cells in low-risk group was dramatically higher than that in high-risk group (P < 0.001). However, the expression of Regulatory T cells (Tregs) in low-risk group was dramatically lower than that in high-risk group (P < 0.0001). In addition, the ssGSEA algorithm revealed totally 10 immune cells that differentially expressed between high-risk group and low-risk group (Figure 4C). For example, the expression of CD4+ central memory T cells in low-risk group was dramatically higher than that in high-risk group (P < 0.05). Furthermore, with P < 0.05 and NES > 1, a total of 13 KEGG pathways that significantly expressed between groups were revealed by GSEA (Figure 4D). The heatmap of hallmark gene sets analysis showed that these pathways were separated by different risk groups (Figure 4E).

Figure 4 Immunological analysis and GSEA investigation. (A) Comparison of infiltration in different risk groups using violin plots for immune scores and stromal scores: X-axis represented different risk group, while the Y-axis represented different scoring. (B) The correlation between expression of immune cells and different risk group based on CIBERSORT algorithm. (C) The correlation between expression of immune cells and different risk group based on ssGSEA algorithm. (D) The significantly different KEGG pathways revealed by GSEA. (E) The heatmap showed the distribution of different KEGG pathways in two risk groups by using hallmark gene sets analysis. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Drug Sensitivity and Immunotherapy Response Analysis

Based on TCGA dataset, the IC50 values of drugs associated with PRAD expression were revealed. The result showed the association between IC50 values of 138 drugs and different risk groups. For instance, BL.D1870, SL.0101.1, and RO.3306 demonstrated notably reduced IC50 values in the high-risk cohort compared to the low-risk cohort, suggesting heightened drug sensitivity among high-risk PRAD patients (Figure 5A). Moreover, the variation in immunotherapy response between two risk groups was explored using TIDE score analysis. The results demonstrated that the TIDE score did not exhibit a marked disparity between the two risk groups (Figure 5B). In addition, except for BTLA, SIRPA, TNFRSF9, and VTCN1, the expression of the remaining checkpoints was dramatically different between two risk groups (Figure 5C) (all P < 0.05). Furthermore, except for TMB and CYT, the score of MSI and TLS between high-risk group and low-risk group was dramatically different (all P < 0.05) (Figure 5D).

Figure 5 The results of drug sensitive and immunotherapy response investigation in current study. (A) Comparison of differences in IC50 levels among three chemotherapy drugs between two risk groups. (B) The comparison of differences in TIDE scores among different risk groups. (C) The expression of immune checkpoint genes among different risk groups. (D) Comparison of differences in Tumor Mutational Burden (TMB), Cytolytic Activity Score (CYT) scores, Microsatellite Instability (MSI), and Tertiary Lymphoid Structures (TLS) scores among different risk groups. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; ns, no significant difference.

Single-Cell Analysis Revealed the Differential Expression of Signatures in the Immune Micro-Environment

Based on GSE176031 from the TISCH database, the expression of six signatures in the immune microenvironment. The result showed that a total of 23 clusters (Figure 6A) and 9 immune cell types (Figure 6B and C) were revealed in current analysis. Immune microenvironment analysis revealed that the six signatures identified in current study were expressed in almost all cell types (Figure 6D).

Figure 6 Single-cell analysis revealed the differential expression of signatures in the immune micro-environment. (A) A total of 23 clusters revealed in GSE176031 dataset: the different color represented different clusters. (B) A total of 9 immune cell types revealed in GSE176031 dataset: the different color represented different immune cell types. (C) The number and distribution of cell types. (D) The expression of six signatures including TRIM14, PIK3R1, PABPN1, BCL2L12, LAPTM4B and TOP2A by using immune micro-environment analysis.

Clustering Analysis of DRGs Associated with 6-AA in PRAD

Based on the calculation of Cumulative Distribution Function (CDF), totally three clusters of DRGs associated with 6-AA in PRAD (Figure 7A) were identified. The heatmap analysis confirmed the result of clustering (Figure 7B). The results of clusters distribution were confirmed by the Phase-Amplitude Coupling (PAC) (Figure 7C). The PCA analysis showed that all genes could be separated by different clusters (Figure 7D). The patients with PRAD in cluster 3 had a better clinical survival prognosis than in cluster 1 and cluster 2 (Figure 7E). The heatmap analysis showed that all six signatures could be separated by different groups according to different clinical characteristics (Figure 7F). The box plot analysis for six signatures showed that all these genes were dramatically differentially expressed between cluster 1 and cluster 2 (all P < 0.0001) (Figure 7G). Moreover, the Sankey diagram analysis showed that high-risk patients were more frequently associated with cluster 1/2 and had higher Gleason scores, indicating a more aggressive disease profile and higher drug resistance. Meanwhile, low-risk patients were more prevalent in cluster 3 which correlating with lower Gleason scores, indicating a better clinical survival prognosis and lower drug resistance than cluster 1 and cluster 2 (Figure 7H). Finally, the distribution differences of RS among different clusters showed that cluster 2 had a high RS than both cluster 1 and cluster 3 (all P < 0.0001) (Figure 7I).

Figure 7 The clustering analysis of PCD-related genes (DRGs) associated with 6-AA in PRAD. (A) The clustering analysis based on DRGs associated with 6-AA in PRAD by using consistency clustering empirical cumulative distribution function (CDF). (B) The heatmap analysis for two clusters. (C) The results of Phase-Amplitude Coupling (PAC) for clusters distributed. (D) The results of Principal Component Analysis (PCA) for clusters. (E) The KM survival analysis for clinical survival prognosis of clusters. (F) The clinical heatmap for signatures expression in current study: the red color represented up-regulation, while the blue color represented the down-regulation. (G) The box plot analysis showed that expression of six signatures in different clusters. (H) The Sankey diagram of different clusters, different risk groups, and clinical distribution. (I) The distribution differences of risk scores among different clusters. ***P < 0.001; ****P < 0.0001.

The qRT-PCR Analysis

The relative expression of TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B in cultured LNCaP cells and DU145 cells were investigated based on qPCR analysis (Figure 8A). The result showed that when compared with LNCaP cells, the expression of TOP2A, PABPN1, BCL2L12 and TRIM14 were all significantly increased in DU145 cells (all P < 0.001). Meanwhile, the expression of PIK3R1 and LAPTM4B in DU145 cells were significantly decreased than those in LNCaP cells (all P < 0.001). In addition, the result of qRT-PCR analysis based on clinical samples showed that when compared with 6-AA sensitive group, the expression of TOP2A, PABPN1, BCL2L12 and TRIM14 were all significantly increased in 6-AA resistant group (all P < 0.001). Meanwhile, the expression of PIK3R1 and LAPTM4B in 6-AA resistant were significantly decreased than those in 6-AA sensitive group (all P < 0.01) (Figure 8B). The expression of six signatures in the validation analysis aligned with the results from our ongoing bioinformatics investigation affirming the reliability of our results.

Figure 8 The results of qPCR analysis for the expression of signatures including TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B. (A) The expression of six signatures between cultured LNCaP cells (6-AA sensitive group) and DU145 cells (6-AA resistant group). (B) The expression of six signatures between 6-AA sensitive samples and 6-AA resistant samples. **P < 0.01; ***P < 0.001.

Discussion

Drug resistance is one of the key factors contributing to the high recurrence rate and poor prognosis of PRAD clinical treatments. Although recent molecular studies, including those on PCD, have partially elucidated the mechanisms of drug resistance in PRAD, the mechanisms of 6-AA resistance in PRAD and the role of PRGs in predicting recurrence and prognosis of PRAD remain unexplored. In this study, we investigated totally 57 DRGs that associated with 6-AA resistance. A total of 6 signatures including TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B were explored, followed by expression validation by qPCR. In addition, the enrichment analysis showed that these signatures such as CBL2L12 were mainly enriched in GO functions like regulation of cysteine-type endopeptidase activity. The nomogram and survival analysis proved the prognostic value of signatures. Immune infiltration analysis based on two algorithms revealed the dyregulation of memory CD4+ T cells between different risk groups.

The prognostic genes play crucial roles in cancer recurrence and immune processes, and their study contributes to revealing the pathogenic mechanisms of tumors and identifying new therapeutic targets.³¹ DNA topoisomerase II Alpha (TOP2A) is a gene encoding DNA topoisomerase IIα. A previous study indicates that there is a close association between TOP2A and prognosis of PRAD.³² Overexpression of TOP2A has been associated with cancer cells’ resistance to certain chemotherapy drugs.³³ Some studies suggest that TOP2A plays a role in regulating apoptosis and DNA damage repair processes, potentially linking it to PCD.³⁴ Phosphoinositide-3-Kinase Regulatory Subunit 1 (PIK3R1) encodes the regulatory subunit of phosphoinositide 3-kinase. This pathway is often aberrantly activated in cancer and is associated with drug resistance. The abnormal expression of PIK3R1 is associated with poor prognosis in PRAD.³⁵ It has been proved that abnormal expression of PIK3R1 affects the activity and infiltration of immune cells in the tumor microenvironment, thus further intervening in the progression of human cancers, including PRAD.³⁶ Lysosomal Protein Transmembrane 4 Beta (LAPTM4B) is a gene involved in tumor invasion and recurrence.³⁷ Its overexpression may affect tumor cell sensitivity to therapeutic agents through mechanisms such as modulation of apoptosis and drug transport.³⁸ A decade ago, LAPTM4B is already proved to be a useful marker for the poor prognosis of PRAD.³⁹ Tripartite Motif Containing 14 (TRIM14) is a gene involved in regulating cell proliferation and immune response.⁴⁰ A previous study shows that aberrant expression of TRIM14 may be associated with cancer cell resistance and responses to therapy-induced cell death.⁴¹ Lysosomal Protein Transmembrane 4 Beta (PABPN1) is a gene encoding polyadenylate-binding protein, playing a crucial role in RNA synthesis and mRNA stability.⁴² It has been found to be associated with PRAD recurrence, and its overexpression is correlated with poor prognosis.⁴³ While the exact role of PABPN1 in drug resistance remains to be fully elucidated, its involvement in mRNA stability and transcriptional regulation suggests potential impacts on cellular responses, including responses to therapeutic agents. BCL2 Like 12 (BCL2L12) encodes a protein that inhibits cell apoptosis, and its increased expression correlates with drug resistance and recurrence in cancer.⁴⁴ It has been confirmed that overexpression of BCL2L12 may suppress tumor cell apoptosis, thereby affecting the tumor’s immunogenicity.⁴⁵ In this study, a total of 6 key genes including TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B were explored to be involved in the regulation of 6-AA resistance in PRAD patients, which indicating that TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B might play a crucial role in mediating 6-AA resistance in PRAD. In addition, TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B were successfully used to construct the prognostic model for predicting the prognosis PRAD. The performance of this prognostic model was confirmed in both TCGA and GEO cohorts, with Kaplan-Meier analysis results showing significant differences in survival between the low-risk and high-risk groups. Also, the ROC curve for the prognostic model showed that the AUC value of 1-year, 3-years and 5-years were all above 0.712 in both TCGA and GEO cohorts, which showed a well diagnostic value for the prognostic model. Importantly, the qRT-PCR proved that the expression of six signatures in the validation analysis aligned with the results from our ongoing bioinformatics investigation, affirming the reliability of our results. Taken together, the above findings and evidence indicated that TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B were crucial in 6-AA resistance development in PRAD, targeting these genes will facilitate their usage in overcoming 6-AA resistance in PRAD.

A previous study shows that certain DEGs such as BIRC5 contribute to the drug resistance via participating in regulation of cysteine-type endopeptidase activity pathway in human cancers including chronic myelogenous leukemia, indicating an important role of this pathway in the progression of cancer and associated clinical drug therapy.⁴⁶ The regulation of cystaine-type endopeptidase activity pathway involves the regulation of the activity of a class of enzymes called cystaine-type endopeptidases. These enzymes include many important protease families, such as caspases, which are particularly important in apoptosis and other cell death processes.⁴⁷ By regulating the activity of cystinase, cancer cells can resist apoptosis induced by chemotherapy and radiotherapy, and thus obtain drug resistance. For example, overexpression of IAPs is common in a variety of cancers including PRAD, and this overexpression can inhibit cystinase activity and lead to anti-apoptosis.⁴⁸ BCL2L12, a member of the Bcl-2 protein family, functions as an anti-apoptotic factor by inhibiting cysteine-type endopeptidases. By stabilizing the mitochondrial membrane, BCL2L12 prevents cytochrome c release, a key event in intrinsic apoptotic signaling.⁴⁹ This stabilization diminishes the downstream activation of caspase-9 and the apoptosome complex, further contributing to drug resistance.⁵⁰ A previous study based on combination treatment of Cisplatin and Titanium Dioxide Nanoparticles indicates that BCL2L12 gene is a therapy target in the clinical treatment of PRAD.⁵¹ However, the detailed molecular mechanism of BCL2L12 via cystaine-type endopeptidase activity pathway during the development of PRAD is still unclear. In the current study, we elucidated the pivotal role of BCL2L12 in mediating resistance to 6-AA through the modulation of cysteine-type endopeptidase activity. Thus, BCL2L12 may regulate PRAD progression and induce 6-AA resistance by interacting with these pathways. By inhibiting key apoptotic pathways, BCL2L12 contributed to the survival and chemoresistance of PRAD cells. Targeting this pathway may provide a novel approach to sensitizing PRAD to apoptosis-inducing therapies, offering a promising direction for future research and clinical management.

Memory CD4+ T are essential contributors to the immune system, particularly in anti-tumor immune responses.⁵² Studies have shown that the presence and functional status of memory CD4+ T cells in the prostatic microenvironment significantly impact tumor immune surveillance and escape mechanisms.⁵³ However, in prostate cancer, the tumor microenvironment often suppresses the function of memory CD4+ T cells, leading to reduced anti-tumor efficacy.⁵⁴ In PRAD treatment, 6-AA may exert its immunomodulatory effects by enhancing the anti-tumor activity of memory CD4+ T cells. For example, 6-AA can increase the expression of tumor antigens, promoting the recognition and attack of tumor cells by memory T cells.⁵⁵ The intrinsic pathway primarily revolves around alterations in mitochondrial membrane potential and the liberation of cytochrome c, whereas the extrinsic pathway triggers apoptosis by engaging death receptors like Fas and TNF receptors.⁵⁶ Activation of these pathways not only directly leads to tumor cell death but also may release tumor antigens, lead to promoting the recruitment and activation of immune cells and further enhancing the anti-tumor immune response.⁵⁷ In this study, we found that both CIBERSORT and ssGSEA algorithms indicated differential expression of memory CD4+ T cells between two risk groups defined by six signatures. Thus, we speculated that memory CD4+ T cells might play take part in tumor microenvironment and influencing the efficacy of 6-AA in PRAD treatment. Although this study provides valuable insights, there are still some limitations. Firstly, the datasets used in this study were obtained from public databases, and as there were differences in data quality, the results were easily affected, resulting in bias. Secondly, the study was based on existing gene expression data and requires further experimental validation. Besides, the specific functions and mechanisms of the drug-resistant gene clusters (TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B) needed to be further investigated by in vivo and in vitro experiments. In addition, whether BCL2L12 takes part in the resistance to 6-AA in PRAD via the cysteine-type endopeptidase activity pathway should be deeply explored. Future research should include functional experiments and clinical trials to validate the specific roles and clinical application potential of these gene clusters in prostate cancer.

Conclusion

In conclusion, PCD-related genes, TOP2A, PABPN1, BCL2L12, TRIM14, PIK3R1 and LAPTM4B, were crucial in 6-AA resistance development in PRAD, targeting these genes will facilitate their usage in overcoming 6-AA resistance in PRAD. The prognostic risk prediction model established in this study could effectively predict patients’ prognosis and provide a reference for the clinical diagnosis and immunotherapy of PRAD.