This study enrolled two independent PCa cohorts. Firstly, Gene expression (RNA-seq) transcriptome (raw counts and transcripts per million reads (TPM)) and sample clinical profiles were downloaded from The Cancer Genome Atlas (TCGA) database (2022.5.1) (https://portal.gdc.cancer.gov/). Certain clinical characteristics, such as Gleason score, which were not accessible in the Genomic Data Commons (GDC) data portal were curated from UCSC XENA (https://xenabrowser.net/datapages/). In total, 491 patients were incorporated into this study. GSE70770 cohort, including 203 PCa patients with microarray expression profile and matched clinicopathological information were curated from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/). The baseline characteristics of patients in the two cohorts are shown in Supplementary Table S1. Chromatin regulators (CRs) were summarized from the previous literature (Lu et al., 2018).

We performed DESeq2 algorithm using DESeq2 R package (R v1.36.0) to explore differentially expressed genes between the normal prostate and tumor part. LfcShrink and “ashr” method was used to generate more accurate estimates (Love et al., 2014). Differential expressed chromatin regulator genes (DE-CRGs) were a

cquired by comparing normal and tumor part in TCGA and filter criteria was |Log2Foldchange| greater than 1.2 and FDR lower than 0.05. “ClusterProfiler” R package (R v4.4.4) was applied to perform Gene Oncology (GO) functional enrichment analysis (Yu et al., 2012). Furthermore, we visualized co-expression potential using the R package corrplot (R v4.1.2), to easily present co-expression and anti-correlation between genes.

Besides, chi-square test was also utilized to assess clinicopathological characteristics distribution between two CRs-related subtypes. Next, we obtained ten oncogenic pathways which are somatically altered in varying cancers. And ssGSEA algorithm was performed to explore the distinct of the ten oncogenic pathways between two CRs-related subtypes (Sanchez-Vega et al., 2018). Moreover, the somatic mutation transcriptome of PCa patients was downloaded from the TCGA database. We analyzed the somatic mutation data by using the “maftools” R package (Mayakonda et al., 2018). We further applied a metagene approach utilized previously for immune cell subpopulations for PRAD tumor microenvironment evaluation (Bindea et al., 2013). The gene set variation analysis (GSVA) method was used to estimate the relative infiltration score of immune cells (Hänzelmann et al., 2013). Furthermore, immune profile differences between two subgroups were assessed by Wilcoxon test.

Firstly, the differentially expressed genes (DEGs) between two subtypes were obtained by using the R package “limma” (R v3.52.2) (Foldchange> 1.5 and FDR <0.05) (Ritchie et al., 2015). Then we intersected these DEGs with differentially expressed genes between the normal prostate and tumor part in TCGA database for subsequent analysis. We then used univariate Cox regression analysis to identify genes correlated with RFS in the training cohort (p < 0.05). The least absolute shrinkage and selection operator (LASSO) algorithm was applied to remove the overfitting between the prognosis-associated genes and reduce the scope of the prognosis-associated genes with penalty parameter tuning conducted via 10−fold cross−validation according to the R package “glmnet”. Next, the genes curated from LASSO regression analysis were incorporated in the multivariate Cox regression analysis. The signature risk score was calculated according to the average expression of each gene and matched regression coefficients generated from multivariate Cox regression analysis. The risk score formula was calculated as follows:

Risk score = betagene1× exprgene1+ betagene2×exprgene2+ betagenen × exprgene.

Then the patients were clustered into high- and low-risk groups based on the median risk score. In addition to the Kaplan–Meier survival curves, we also applied the time-dependent receiver operational feature curves (ROC) generated by the R packages “suvminer” and “survival ROC” to evaluate the performance of the signature risk score in predicting the RFS of PCa patient.

Independent prognostic factors and related clinical parameters were obtained from Cox stepwise regression analysis for constructing a prognostic nomogram to predict 1-, 2-, and 3-year RFS for PCa patients. Then we used Calibration curves of 1-, 2-, and 3-year to evaluate the reliability of this nomogram. Moreover, we also utilized the decision curve analysis (DCA) to compare the performance of clinical parameters and the nomogram model.

We exerted chi-square tests on the related clinical characteristics in different risk groups including T stage, Gleason score, and clinical outcome. CIBERSORT is an online tool for quantifying the infiltration abundance of 22 types of immune cells according to the basement of linear support vector regression (Chen et al., 2018). Aiming to estimate the immune cell function between the different risk groups, we used CIBERSORT to calculate the abundance of Tumor-infiltrating immune cells (TIICs) in PRAD. We further drew a comparison in the abundance of tumor infiltrating immune cells between the high- and low-risk groups by Wilcoxon test. Besides, the ESTIMATE algorithm was utilized to assess immune infiltration in PRAD patients (Li et al., 2017). Moreover, a Spearman correlation analysis between immune infiltrating cells and risk score or core genes constituting the risk score was performed. Additionally, we also performed chi-square test analysis to explore the TMB differences in the HR and LR subtype.

To explore the differences of molecular characteristics between the CRGs subgroups, we performed the gene set enrichment analysis (GSEA) based on the “ClusterProfiler” R package. GSVA enrichment analysis was also carried out to compare the transcriptomic remodeling using GSVA R package (R v1.44.2). “Limma” package was also utilized to nominate the distinct molecular pathways between the two subgroups (p < 0.05 was considered significant). Besides, we used “pRRophetic” R package to explore the therapeutic sensitivity and the concentration inducing 50% reduction growth (IC50) of targeted inhibitors. Then we performed Wilcoxon test to compare the IC50 difference between the two subgroups.

The human Pca cell lines PC3 and C4-2 were purchased from the American Type Culture Collection (Manassas, United States) and cultured in RPMI-1640 medium (Corning, Inc., Corning, NY, United States) containing 10% fetal bovine serum (GIBCO) and 1% penicillin/streptomycin (GIBCO). All cells were grown at 37 °C in a 5% CO2 humidified incubator.

Short hairpin RNA (shRNA) expression sequences are documented in Additional file: Table S3. Then these sequences were cloned into the pLKO.1 vector. And plasmids were transfected into HEK293FT cells utilizing PEI 25K (23966–1; Polysciences, Warrington, PA, United States) based on the manufacturer’s instructions. PC3 and C4-2 cells were transduced with lentivirus, and stable transformants were selected with puromycin (5 μg/ml) for 7 days.

RNA was isolated using TRIzol reagent (Invitrogen) in accordance with the manufacturer’s instructions. Then 1 ug of total RNA was reverse transcribed into cDNA using a PrimeScript™ first Strand cDNA Synthesis Kit (6110A; TaKaRa). qRT-PCR was conducted using TB Green Premix ExTaq (Tli RNaseH Plus) (RR420; TaKaRa) on the ABI7500 System (Applied Bio Systems, Foster City, CA, United States). Then we used 2 − ΔΔCt method to calculate the relative expression levels of genes. GAPDH was considered as an internal control for RT-PCR. Primers used to amplify genes of interest were listed in Additional file Table S3.

Cell growth was detected by Cell Counting Kit-8 (CK04; Dojindo, Kumamoto, Japan) at indicated time points according to the manufacturer’s instructions. The cells were cultured in 96-well plates (1000 cells per well) for 6 days. Ten μL CCK8 reagent was added to 100 μl complete medium in each well and then cultured at 37 °C for 3 h. Then the absorbance values were assessed using a microplate reader (Tecan, Mechelen, Belgium) at 450 nm (A450). Colony formation assay was conducted by seeding 1000 cells in complete medium for 10–12 days depending on colony size. Then the cells were fixed using methanol for 10 min and stained using 0.5% crystal violet for 1 h. Images were captured following wash of PBS.

Statistical analysis in this study was performed using R software v4.0.4 and Prism software, version 8 (GraphPad Software, San Diego, CA, United States). p values less than 0.05 were thought to be statistically significant.

The main analysis workflow is presented in Figure 1. Firstly, we identified differentially expressed genes (DEGs), including 2001 upregulated genes and 1928 downregulated genes by performing differential analysis between prostate tumor and normal part acquired from The Cancer Genome Atlas (TCGA) dataset (Figure 2A). Subsequently, we obtained 36 differentially expressed CR-related genes based on the intersection of CR-related genes and DEGs in TCGA (Figure 2B). Gene Ontology enrichment analysis revealed that the above differentially expressed CR-related genes were mainly enriched in “histone modification” in the biological process category, “condensed chromosome” in the cellular component category and “hydrolase activity, acting on carbon nitrogen (but not peptide) bonds” in the molecular function (Figure 2C). Univariate Cox analysis suggested that 18 CR-related genes are significantly associated with recurrence-free survival (RFS) (Figure 2D). Furthermore, correlation analysis suggested that most genes were significantly associated with each other (Figure 2E).

For exploration of CRs-related genes heterogeneity in PCa, we applied consensus clustering analysis to construct CRs-related molecular clusters of Pca using a TCGA cohort. The results of the consensus package in R were used to identify two CRs-related subtypes (Figure 3A) (Wilkerson and Hayes, 2010). We conducted principal component analysis (PCA) using “prcomp” function in R (v4.0.4) to verify the distinct classification between Cluster 1 and Cluster2 (Supplementary Figure S1). The two subtypes showed distinct clinical outcomes. Kaplan–Meier (K-M) plots revealed that patients in cluster 2 exhibited inferior RFS (Figure 3B). Similarly, the GSE70770 cohort could also be divided into two subgroups based on the CRs-related genes expression (Supplementary Figure S2A), and K-M plots also showed a similar difference in RFS between the two subtypes (Supplementary Figure S2B). Additionally, comparison of the two subtypes based on clinical parameters revealed that cluster 2 contained a greater proportion of patients with a higher Gleason Score (GS), advanced T stage and recurrent status present more proportions in cluster 2 (Figure 3C). Next, we performed single-sample gene set enrichment analysis (ssGESA) method to compare ten distinct oncogenic pathways between the two subtypes. Notably, the score of oncogenic pathways involved in the deterioration of prostate cancer, including Cell Cycle, MYC and PI3K-AKT pathways, were significantly elevated in cluster2 (Figure 3D). Collectively, these results implied that the clinical and molecular characteristics of PCa were more aggressive among patients in cluster 2.

We obtained simple nucleotide variation data from TCGA to explore differences in genomic mutations between CRs-related subtypes. The top 10 genes with the highest mutation frequencies are presented in Figure 4A; TP53, SPOP, and FOXA1 were more frequently mutated in cluster 2. Additionally, we evaluated the infiltration abundances of immune cells using the ssGSEA method to determine the association of each CR-related subtype with immune status. Patients in cluster 2 showed higher infiltration abundances of activated CD4 T-cell, CD56^dim natural killer cells, gamma-delta T-cell, and type 2 helper cells (Figure 4B). Furthermore, the ESTIMAT score further revealed that cluster2 had lower immune and stromal scores compared to the cluster1 (Figure 4C). Comprehensively, the results implied that CRs-related subgroups of prostate cancer can accurately indicate immunity status.

For further investigation of mechanisms that contribute to the heterogeneity of CRs-related subtypes and application of these subtypes to clinical prediction and treatment, we firstly performed differential analysis of the two subtypes. We recognized 963 DEGs that were correlated with the two CRs-related subtypes (fold-change >1.5 and false discovery rate <0.05) (Figure 5A). The intersection of these DEGs with DEGs identified through comparisons of normal prostate and tumor tissue in TCGA yielded 483 DEGs (Figure 5B). Subsequently, univariate Cox analysis identified 261 genes that were correlated with RFS (supplementary Table S2). Next, we used lasso regression analysis to optimize the number of genes; we selected the six genes listed in Figure 5C. Finally, 4 genes including MXD3, SSTR1, AMH, and PPFIA2 were obtained based on the stepwise cox multivariate regression analysis. Then we constructed a 4 genes prognostic signature to predicate the RFS of PCa patients. Then, patients in the TCGA and Gene Expression Omnibus (GEO) cohorts were clustered into a high-risk group (HRG) and low-risk group (LRG) based on the following risk score (RS) formula:

Risk score = 0.796*MXD3+0.138*AMH+0.137*SSTR1+0.15*PPFIA2.

We used the median RS value as the cut-off for classifying patients into the HRG and LRG. K-M survival analysis indicated that the HRG exhibited a worse prognosis than the LRG in both the training and validation cohorts (Figures 6A, C). For additional assessment of risk model accuracy, we conducted receiver operating characteristic (ROC) analysis in both the training and validation cohorts. In the training cohort, the areas under the ROC curve of 1-year, 2-year, and 3-year RFS were 0.79, 0.76, and 0.76, respectively. In the validation cohort, the areas under the ROC curve of 1-year, 2-year, and 3-year RFS were 0.71, 0.65, and 0.64, respectively (Figures 6B, D).

Univariate and multivariate Cox analyses were conducted to assess the relationships of RFS with potential variables. Pathological T stage, RS and GS were identified as independent risk factors (Supplementary Figure S3A, B). Consequently, a nomogram with an integrated prognostic risk score model, pathological T stage and GS was constructed for RFS prediction in PRAD samples from prostate adenocarcinoma patients (Figure 7A). The calibration curves at 1, 2, and 3 years showed good linearity and suggested that the nomogram could accurately predict the RFS of PRAD patients (Figure 7B). Furthermore, the decision curve analysis (DCA) suggested that, compared with clinical parameters such as GS or pathological T stage, the nomogram showed superior net clinical benefit (Figure 7C).

To gain insights into the correlation between clinical features and risk model, we performed Chi-square test to study the distribution of patients in the HR group and LR group based on clinical characteristic including GS, T stage and clinical outcome. Results from the analysis provided the evidence that in the TCGA cohort HR group was significantly associated with higher GS, aggressive T stages and poor prognosis. Furthermore, patients in cluster 2 occupied a greater proportion of the HR group (Figures 8A, B). Additionally, to determine whether the signature score was associated with tumor immunity, we compared the numbers of distinct tumor-infiltrating immune cells between the two subtypes using the CIBERSORT algorithm. The results suggested that, compared with the LRG, the HRG contained more regulatory T-cell (p < 0.01) and more M2 macrophages (p < 0.01). In contrast, low-risk group had more plasma cells (p < 0.01) and more resting mast cells (p < 0.05) (Figure 9A). Notably, the signature score was positively associated with the enrichment scores of regulatory T-cell and M2 macrophages (Figures 9C, D). Conversely, the RS was negatively correlated with plasma cells and mast cells (Figures 9E, F). Previous studies demonstrated that high levels of infiltrating M2 macrophages and regulatory T-cell were correlated with biochemical recurrence (Andersen et al., 2021). This partly explain the poor prognosis of PCa patients in the high-risk group. Importantly, correlation analysis between risk-related genes and tumor-infiltrating immune cells revealed that MXD3 had the strongest positive correlation with the number of regulatory T-cell (Figure 9B).

To investigate correlation between the CRs-related signature and the mutational landscape, we calculated the tumor mutational burden (TMB). Compared with patients in the LRG, patients in the HRG exhibited higher TMBs (Figure 10A). Strikingly, K-M plots suggested that patients with higher TMBs were more likely to exhibit progression, compared with patients who had lower TMBs (Figure 10B). Furthermore, use of combined risk models showed that patients in the high-risk + high-TMB group exhibited the worst prognosis, according to K-M survival analysis (p < 0.01) (Figure 10C). These findings indicated that both the risk score and TMB can predict poor prognosis in Pca patients.

Gene set enrichment analysis (GSEA) suggested that patients in the HR and LR groups have different transcriptomic alterations. Gene Ontology terms enriched in the HRG were myc targets V2, myc targets V1, oxidative phosphorylation, DNA repair, E2F targets, and G2-M checkpoint (Figure 10D). Similarly, gene set variation analysis (GSVA) revealed that Gene Ontology terms enriched in the HRG were cell cycle, base excision repair, mismatch repair, DNA replication, and homologous recombination (Figure 10E). These findings indicated that the risk model was closely associated with cell cycle-related pathways or DNA repair-related pathways, which require the participation of multiple CRs. Because a high RS is associated with poor prognosis and multiple oncogenic signaling pathways contribute to progression in PCa patients, we used the pRRophetic package in R to explore the relationship between RS and potential targeted inhibitors. As shown in Figure 9F, high-risk score samples were more sensitive to Erlotinib, Sunitinib, VX-680, TAE684 and Crizotinib. These drugs may be used as alternative treatment for PCa progression in high-risk patients.

Because MXD3 is a hub gene with a central role in the PCa signature, we focused on MXD3 during in silico and in vitro analyses. Analysis of TCGA PCa data revealed that MXD3 was strongly upregulated in PCa tissue (Figure 11A). Kaplan-Meier analysis of TCGA data indicated that the level of MXD3 expression was significantly associated with RFS of Pca patients (Figure 11B). Moreover, K-M analysis of mCRPC patient data from the West Coast Prostrate Cancer Dream Team (WCDT) cohort (Quigley et al., 2018) showed that high expression of MXD3 was strongly associated with overall survival (Figure 11C). Additionally, we analyzed the differential expression level of the MXD3 in various pathological stages and Gleason score of Pca patients using TCGA data; the findings indicated that MXD3 is significantly upregulated in higher T stage and higher Gleason score groups (Figure 11D). For additional exploration of the biological function of MXD3 in PCa, we performed in vitro experiments to validate the oncogenic role of MXD3 in the PC3 PCa cell line. We silenced MXD3 in PC3 cells and used reverse transcription polymerase chain reaction to confirm MXD3 knockdown (Figure 11E). Cell proliferation was evaluated using the CCK-8 method, and the results suggested that MXD3 knockdown significantly reduced the growth of PCa cell (Figure 11F). Colony formation assays indicated that MXD3 inhibition considerably reduced the numbers of PC3 cell colonies (Figures 11G, H). Considering the role of AR in the progression of prostate cancer, we selected the C4-2 PCa cell for further validation. Similarly, MXD3 knockdown significantly reduced the growth of C4-2 (Supplementary Figure S4). Collectively, these results demonstrated that MXD3 is essential for growth of PCa cell.