Transcriptomic data and corresponding clinical features [including age, tumor (T) and node (N) stage, Gleason grade and relapse-free survival (RFS)] for The Cancer Genome Atlas (TCGA)-prostate adenocarcinoma (PRAD) were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/), with a total of 551 samples. After screening, 532 samples were obtained, with PCa:Control=481:51. A total of 409 of the PCa samples had survival information and biochemical recurrence (BCR) information. GSE103512 (GPL13158), GSE21034 (GPL10264), GSE70770 (GPL10558) and GSE54460 (GPL11154) were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/), and each dataset was subjected to independent analysis. GSE103512 comprised 60 PCa samples along with seven control samples. Similarly, GSE21034 and GSE70770 contained PCa:Control=150:29 and PCa:Control=219:74, respectively, and GSE54460 contained 106 PCa samples, of which 96 samples without a recurrence time 0 were retained. The sample information of the five datasets are presented in Table SI. In addition, 982 DRGs were identified from the literature (26) and 78 EHRGs were obtained from the Gene Ontology (GO) database (GO:0097009; R package org.hs.eg.db; http://bioconductor.org/packages/org.Hs.eg.db/).

To identify DEGs between PCa and control samples in the TCGA-PRAD dataset, the DESeq2 (v 1.36.0) package (https://bioconductor.org/packages/DESeq2/) was employed (adj.P<0.05 and |Log₂ fold change|>0.5). Subsequently, the ggplot2 (v3.3.6; http://CRAN.R-project.org/package=ggplot2) and pheatmap (v1.0.12) packages (https://CRAN.R-project.org/package=pheatmap) were utilized to plot the volcano map and heatmap, respectively (27). Finally, the clusterProfiler (v4.7.1.001) package (https://bioconductor.org/packages/clusterProfiler/) analyzed the biological functions for DEGs for GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment (P<0.05). The results were then visualized using the GOplot (v1.0.2; http://CRAN.R-project.org/package=GOplot) and enrichplot (v1.16.2) (https://bioconductor.org/packages/enrichplot/) packages.

Candidate biomarkers were identified from the intersection between DEGs, DRGs and EHRGs using the ggvenn (v0.1.9) package (https://CRAN.R-project.org/package=ggvenn). Additionally, the Wilcoxon rank-sum test was used to assess the expression level of candidate biomarkers between PCa and control groups in the TCGA-PRAD, GSE103512, GSE21034 and GSE70770 datasets (P<0.05). The mRNA expression and protein contents of candidate biomarkers were determined using western blotting and reverse transcription-quantitative (RT-q) PCR, respectively. Finally, genes with consistent expression and significant difference in both four datasets and experiments (western blotting and RT-qPCR) were identified as biomarkers.

According to the expression levels of biomarkers among PCa samples in the TCGA-PRAD dataset, they were divided into high- and low-expression groups, utilizing the optimal cut-off value for expression as the boundary using the surv_cutpoint function in the survminer (v0.4.9) package (https://CRAN.R-project.org/package=survminer). Kaplan-Meier (KM) survival curves for the expression groups were then plotted using the survminer (v0.4.9) package to determine the difference in overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI) progression-free interval (PFI) and RFS between the two expression groups. Meanwhile, the survival difference between two expression groups was compared using the log-rank test. Furthermore, the difference in biomarkers between different subtypes of clinical characteristics were compared in the TCGA-PRAD dataset and the results was visualized using a violin plot using the ggplot2 (v3.3.6) package (P<0.05). Finally, the χ² test was used to assess the distribution of several clinical features between two expression groups.

To elucidate the biological processes associated with the identified biomarkers which were involved in PCa, the c2.cp.kegg.v2023.2.Hs.symbols.gmt gene set was used as reference gene set from the Molecular Signatures Database (MSigDB; http://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp) using the clusterProfiler (v4.7.1.001) package. Spearman's correlation analysis of the biomarkers with all genes in the TCGA-PRAD dataset was performed and they were then ranked according to the correlation coefficient. The ranked genes were then analyzed for GSEA enrichment and the results of the top 5 pathways from largest to smallest were visualized based on their |normalized enrichment score (NES)| values (|NES|>1; P<0.05). To further assess pathway changes between two expression groups in the TCGA-PRAD dataset, and based on the c2.cp.kegg.v2023.2.Hs.symbols.gmt gene set from the MSigDB, the GSVA (v1.44.5; http://bioconductor.org/packages/GSVA/) and limma (v 3.52.4) packages (https://bioconductor.org/packages/limma/) were used to calculate GSVA scores for each pathway and to compare the GSVA score differences of all KEGG pathways between the two expression groups (|t|>2; P<0.05).

To evaluate the changes in the immune microenvironment of PCa, the content and relative infiltration abundance of 22 immune infiltrating cells in the TCGA-PRAD dataset were calculated utilizing the CIBERSORT (v1.03) algorithm (https://cibersort.stanford.edu; P<0.05) and visualized using the ggh4× (v4.2.2) package (https://CRAN.R-project.org/package=ggh4×). Additionally, the Wilcoxon rank-sum test (P<0.05) was used to analyze the difference in immune cell infiltration between the two expression groups, and the ggplot2 (v 3.4.1) package was used to depict the data in a box plot.

The maftools (v 2.12.0) package (https://bioconductor.org/packages/maftools/) was used to analyze two patient cohorts with mutation data to assess the genetic differences between the two expression groups in the TCGA-PRAD dataset.

To predict immunotherapy response in patients with PCa, the tumor immune dysfunction and exclusion (TIDE) score was used to evaluate the response to immunotherapy in the Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard.edu/). Differences between expression groups was assessed utilizing the Wilcoxon rank-sum test (P<0.05). A higher TIDE score indicates a greater immune escape of tumor cells. Subsequently, Spearman's correlation analysis was utilized to determine the correlation between biomarkers and TIDE score.

Furthermore, using the TCGA-PRAD dataset, the estimate (v1.0.13) package (https://R-Forge.R-project.org/projects/estimate/) was utilized to calculate the StromalScore, ImmuneScore and EstimateScore of biomarkers in the expression groups, and the differences in these groups were compared using the Wilcoxon rank-sum test (P<0.05). Lastly, the immune subtypes of disease samples were divided according to the gene expression of immune cells in tumors using the ImmuneSubtypeClassifier (v0.1.0) package (https://github.com/CRI-iAtlas/ImmuneSubtypeClassifier), and the differences among immune subtypes were then compared using the Wilcoxon rank-sum test (P<0.05).

Subsequently, using the oncoPredict (v0.2) package (28), the 50% inhibitory concentration (IC₅₀) values of 198 chemotherapeutic drugs were determined according to the Genomics of Drug Sensitivity in Cancer (GDSC; http://www.cancerrxgene.org/), and the IC₅₀ values of each chemotherapeutic drug were compared between the expression groups using the Wilcoxon rank-sum test (P<0.05).

PCa and adjacent non-cancerous tissues were obtained from the First Hospital of Shanxi Medical University (Taiyuan, China). A total of 6 cancer tissues and 6 adjacent non-cancerous tissue were collected from the pathological specimens of 6 patients after radical prostatectomy. All the participants were male, and the age range of patients was 60–75 years, with an average age of 68±5 years. Ethical approval was granted by the hospital Ethics Committee (approval no. KYLL-2024-100). All participants provided written informed consent to participate in the present study.

Total RNA was extracted using the TRIzol^® method (Invitrogen; Thermo Fisher Scientific, Inc), and cDNA was synthesized following the instructions of the PrimeScript^™ RT reagent Kit (cat. no. RR047A; Takara Bio, Inc.). For RT-qPCR, the reactions were performed using an Applied Biosystems 7300 Real-Time PCR machine (Thermo Fisher Scientific, Inc). The reaction mixture consisted of 10 µl 2X SYBR Green qPCR MasterMixII (Universal; Saiwen Innovation (Beijing) Biotechnology Co., Ltd), 2 µl cDNA, 0.4 µl each of upstream and downstream primers and ddH₂O to a final volume of 20 µl. The specific reaction conditions were the pre-denaturation in the first stage: 95°C, 3 min. The second stage of denaturation, annealing and extension: 95°C, 5 sec; 60°C, 30 sec; 60°C, 30 sec, respectively (40 cycles). Final extension of the third stage: 60°C, 5 min. A total of 12 samples were collected and each sample was tested three times. GAPDH served as the internal control, and the data were analyzed using the 2^−ΔΔCq method (29). The primers used are listed in Table SII.

Protein was extracted from PCa tissues using RIPA lysis buffer (cat. no. P0013B; Beyotime Biotechnology). The bicinchoninic acid protein assay was used to determine protein concentration. The percentage of gel was 10%, 30 µg of the sample was added to each well and was then subjected to SDS-PAGE and transferred to a PVDF membrane (cat. no. 0000279048; MilliporeSigma). To block nonspecific binding, the membranes were incubated with 5% skimmed milk (cat. no. GC310001-100g; Wuhan Servicebio Technology Co., Ltd.) and 0.1% Tween-20 (cat. no. GC204002-100 ml; Wuhan Servicebio Technology Co., Ltd.) in Tris buffer at room temperature for 30 min. RMI1 primary antibodies (cat. no. 14630-1-AP; Proteintech Group, Inc.) were diluted with antibody dilution buffer (1:1,000; cat. no. P0023A; Beyotime Biotechnology) and incubated overnight with the membrane at 4°C, followed by three washes with TBST. Subsequently, the appropriate goat anti-rabbit IgG (H+L) secondary antibodies (cat. no. 31460; Invitrogen; Thermo Fisher Scientific, Inc.) was diluted with 1X TBST (1:3,000) and incubated at room temperature for 1 h. Protein bands were visualized using an ultra-sensitive ECL chemiluminescence kit (cat. no. PK10003-100 ml; Proteintech Group, Inc.), and images were captured using an imaging system (Bio-Rad Laboratories, Inc.). The relative gray values of the protein bands were measured by ImageJ (v2.16.0; National Institutes of Health), using tubulin as the internal control. The molecular weights of RMI1, GAPDH and tubulin are 75, 36 and 55 kDa, respectively.

Human PCa cells PC-3 were purchased from Immocell (Xiamen Yimo Biotechnology Co., Ltd.; cat. no. IM-H075; http://www.immocell.com/?ryxbx/1680.html). The PC-3 cell line was verified using short tandem repeat profiling and passaged in the laboratory for <6 months.

Knockdown in the present study was performed using RNA silencing. The plasmids used in the experiment were constructed by our research group independently. The plasmid concentrations of the short hairpin negative control (shNC) and RMl1-shRNA-1, RMl1-shRNA-2, and RMl1-shRNA-3 were 524, 465, 369, and 350 ng/µl, respectively. Transfection was performed at 37°C for 24 h using the jetPRIME^® transfection reagent (cat. no. 101000025; Polyplus-transfection SA). The medium was replaced, followed by continuous culture for 48 h. Subsequent functional experiments were performed when the cells reached confluence. RT-qPCR and western blotting were used to assess the silencing efficiency of short hairpin (sh)RNA-transfected cells. The targeted sequences of shRNA and the shNC sequence were shown in Table SIII.

Transfected cell suspensions were added to 96-well plates at a volume of 100 µl per well and subsequently incubated with 10 µl CCK-8 solution (Beyotime Biotechnology) for 4 h at 37°C. Absorbance was then measured at 450 nm using a microplate reader (Feyond-A300; Hangzhou Allsheng Instruments Co., Ltd.).

In the scratch assay, cells in the logarithmic growth phase were digested and plated at a designated density into 6-well culture plates. The seeding process aimed to achieve 100% confluence following overnight incubation, using a final medium volume of 2 ml per well. After 24 h of incubation at 37°C with 5% CO₂, the cells were scored perpendicular to the bottom of the wells using a 200 µl pipette tip, with the lid or ruler of the 6-well plate serving as a reference. Subsequently, the cells were washed three times with PBS to eliminate debris, and serum-free medium was added. Images were captured under a ×4 magnification though an inverted fluorescence microscope (cat. no. MF52-N; Mshot) to ensure that the scratch was centered and perpendicular with a consistent background. Test samples of different concentration gradients were added. Samples were taken at 0 and 40 h and images were captured. Results were analyzed using ImageJ software (v1.4.3.67; National Institutes of Health). In the scratched area, 6–8 horizontal lines were randomly drawn, and the average distance between cells was measured to determine the cell migration rate.

The Transwell migration assay was performed using a Transwell chamber. After 24 h of cell culture, a serum-free suspension of 1×10⁵ cells (100 µl) was added to the upper Transwell chamber, while 600 µl complete medium supplemented with 30% FBS (Thermo Fisher Scientific, Inc.) was added to the lower chamber. Following 24 h of incubation at 37°C, a wet cotton swab was used to carefully remove the non-migratory cells from the upper membrane surface, whilst crystal violet was used to stain the migrating cells on the bottom surface at room temperature for 10 min. After another 24-h incubation at 37°C, the process was repeated to ensure thorough staining. The migrating cells were then observed and images were captured under an inverted fluorescence microscope (cat. no. MF52-N; Mshot). Cell counts were performed in three randomly selected fields, and the average number of migrating cells was calculated to assess cell migration capability.

All analyses were performed using R version 4.2.2 (The R Foundation) and GraphPad Prism 9 (Dotmatics). For comparisons between groups, the Wilcoxon rank-sum test, log-rank test and unpaired two-tailed Student's t-test were used. For group comparisons of three or more groups, the Wilcoxon rank-sum test for pairwise comparisons and then adjusted the significance levels using Bonferroni correction was used. Data values were presented as mean ± standard deviation (SD) in the experimental analyses. Statistical analysis of the experimental section employed unpaired two-tailed Student's t-tests and one-way analysis of variance (ANOVA), followed by the least significant difference test. For the comparison of the three sets of data involved in the experiment section, Tukey's test was used for statistical analysis. The experimental data were analyzed using SPSS software (version 26; IBM Corp.). P<0.05 was considered to indicate a statistically significant difference.

The gene expression levels of DEGs between PCa and control samples were assessed. This analysis identified a total of 4,395 DEGs, comprising 2,039 genes that were upregulated and 2,356 genes that were downregulated. The volcano plot and heatmap are presented in Fig. 1A and B. In addition, according to the results of GO analysis, 1,075 biological processes, 102 cellular components and 93 molecular functions were enriched by DEGs. They were mostly engaged in ‘positive regulation of DNA-binding transcription factor activity’, ‘ribosomal subunit’ and ‘regulation of stress-activated protein kinase signaling cascade’ (Fig. 1C). Furthermore, 39 KEGG pathways were enriched by DEGs, such as ‘regulation of actin cytoskeleton’, ‘ribosome’, ‘adrenergic signaling in cardiomyocytes’, ‘MAPK signaling pathway’ and ‘cAMP signaling pathway’ (Fig. 1D). These analyses suggested that the identified DEGs may be involved in a wide range of biological mechanisms underlying PCa development.

Through Venn diagram analysis, it was revealed that growth differentiation factor 15 (GDF15) and RMI1 were jointly screened out at the intersection of DEGs, DRGs and EHRGs (Fig. 2A). The research on the GDF15 gene in the field of PCa has been relatively in-depth, with a number of studies exploring its functions and mechanisms of action (30–34). Therefore, RMI1, whose biological functions in PCa have not yet been fully elucidated, was chosen as the key research target. Further analysis demonstrated that the expression trends of RMI1 were consistent and significant across four datasets, and were revealed to be highly expressed in the PCa group compared to the control group based on the Wilcoxon test results (Fig. 2B-E). Additionally, RMI1 expression demonstrated consistent and significant differences (P<0.05) in both RT-qPCR and western blot analyses, establishing it as a robust biomarker (Fig. 2F and G). The consistent and statistically significant expression patterns of RMI1 observed across multiple datasets and experimental validations suggest that RMI1 might serve a pivotal role in the progression of PCa.

To evaluate the effect of RMI1 expression on patient survival, patients were categorized into high and low expression groups according to their levels of RMI1 expression. Survival outcomes were then compared, including OS, DSS, DFI, PFI and RFS, between these groups. The analysis indicated notable differences in OS, DSS and PFI between patients with high and low levels of RMI1 expression in the TCGA-PRAD dataset (Fig. 3A-C). However, this dataset did not show any significant differences in DFI and RFS (Fig. S1A and B). In addition, although the difference in RFS did not reach statistical significance in the TCGA-PRAD dataset (despite there being a trend towards a lower recurrence-free rate noted in the high expression group with P<0.1), significant differences in RFS were demonstrated in other datasets, specifically GSE70770 and GSE54460 (Fig. 3D-E). This suggested that RMI1 may influence BCR in PCa. Furthermore, notable variations in RMI1 expression levels were observed across different subgroups based on Gleason grading, T-staging and N-staging (Fig. 3F). Furthermore, the χ² test results indicated a significant difference in T-staging among patients with varying levels of RMI1 expression (Fig. 3G).

Enrichment analysis was used to delve deeper into the functions of RMI1. GSEA demonstrated that RMI1 is involved in 71 KEGG pathways, primarily associated with the cell cycle, oocyte meiosis and cancer-related pathways (Fig. 4A). In addition, GSVA revealed significant differences in 143 pathways between the groups with high and low RMI1 expression. Specifically, pathways related to the metabolism of xenobiotics by cytochrome P450 were activated in the low expression group, whilst pathways such as aminoacyl-tRNA biosynthesis were activated in the high expression group (Fig. 4B). Additionally, the role of somatic cell mutations on PCa was assessed using the mutated landscapes. According to the waterfall plot, the low-expression group had the highest frequency of speckle type POZ protein (SPOP) gene mutations (13%), whilst the high-expression group had the highest frequency of TP53 mutations (15%). SPOP and TP53 were both missense mutations (Fig. 4C and D).

The results of immune cell infiltration are presented in Fig. 5A. A total of eight differential immunity cells were identified, such as naive B cells, plasma cells and CD8 T cells (P<0.05; Fig. 5B). At the same time, the TIDE score was markedly greater in the high-expression group (P<0.05; Fig. 5C). In addition, the StromalScore, ImmuneScore and EstimateScore of the high-expression group were all markedly higher than those of the low-expression group (P<0.05; Fig. 5D). This suggests that high levels of the RMI1 may be associated with immune escape and poor prognosis.

Furthermore, based on the 198 chemotherapy/targeted therapy drugs provided by GDSC, 164 drugs with significant differences in their IC₅₀ values were identified in different expression groups. In the boxplot demonstrating differences in IC₅₀ between the expression groups for three common chemotherapy drugs for PCa (5-fluorouracil, oxaliplatin and irinotecan), only irinotecan showed a significant difference between the two expression groups (P<0.001) and had a positive correction with RMI1 (cor=0.37) (Fig. 5E).

KM survival analysis indicated that higher RMI1 expression was associated with a worse prognosis in patients with PCa. As a result, PC-3 cells with an improved prognosis were constructed by transfecting sh-RMI1. The mRNA and protein expression levels of the transfected cells confirmed that the transfection was successful (Fig. 6A and B). Furthermore, the CCK-8 assay results indicated a significant reduction in cell viability for sh-RMI1 cells with decreased RMI1 expression (P<0.05; Fig. 6C), suggesting that lowering RMI1 levels can inhibit the proliferation of PCa cells. Furthermore, the findings from the scratch assay (Fig. 6D) and the Transwell migration assay (Fig. 6E) revealed that cells with reduced RMI1 expression demonstrated markedly decreased migration (P<0.05). These experiments suggested that RMI1 may serve a role in enhancing the proliferation and migration of PCa cells. A reduction in RMI1 expression led to decreased proliferation and migration in these cells, indicating that RMI1 could be crucial for the progression and metastasis of PCa.