Data from bulk RNA-seq, clinical records, and SNP mutation profiles of TCGA-PRAD were obtained from the TCGA database (https://portal.gdc.cancer.gov/). Clinicopathological information for TCGA-PRAD was acquired from the GDC Data Portal (https://portal.gdc.cancer.gov/projects/TCGA-PRAD) and the TCIA clinical repository (https://www.cancerimagingarchive.net/collection/tcga-prad/). The dataset included 52 normal prostate tissue samples and 492 PRAD tumor samples, all of which contain detailed survival and clinicopathological annotations. The scRNA-seq datasets PRAD_GSE137829, PRAD_GSE141445, and PRAD_GSE172301 were also incorporated into this study. In total, scRNA-seq data from 32 patients with PRAD were analyzed, with all datasets retrieved from the TISCH database (http://tisch.comp-genomics.org/gallery/).

Detailed patient information and sequencing statistics are available in Table 1. The integration of samples was performed using the anchors method, a technique implemented in the R “Seurat” (10). Core cells were subsequently created through a filtering step applied to the scRNA-seq data. Unqualified cells refer to genes that exhibit detectability exclusively in three or fewer cells, as well as in less than 200 cells of low quality. Genes identified under these criteria were omitted from subsequent analysis. The analysis of gene expression in core cells involved utilizing a standardized linear regression model, followed by the application of ANOVA to identify the top 2000 genes exhibiting highly variable fingerprints. Principal component analysis (PCA) was conducted on individual cell samples, identifying the first 20 principal components (PCs). These components were then chosen for further study. The UMAP algorithm was employed to comprehensively analyze dimensionality reduction on the initial 20 principal component pairs of the samples (11). R “singleR” was employed to utilize reference datasets such as HumanPrimaryCellAtlasData, BlueprintEncodeData, and ImmuneCellExpressionData for auxiliary annotation (12). Additionally, the CellMarker database and prior research were consulted to identify marker genes for the manual annotation of distinct clusters (13). The visualization process was facilitated by utilizing the R “SCP” (https://github.com/zhanghao-njmu/SCP).

The Seurat package’s FindAllMarkers function was utilized to locate marker genes for each cluster, with parameters set as min.pct=0.2 and only (14). The Wilcoxon rank sum test was performed to identify differentially expressed genes (DEGs) and screen marker genes when pos=TRUE. Using “clusterProfler” in R software (15), core cells of marker genes underwent enrichment analysis for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functions, respectively. ‘Monocle 2’ was utilized to infer the developmental trajectory of cells to elucidate the underlying molecular causes of PRAD progression (16). The “cellchat” algorithm was employed to explore potential cell-to-cell interactions, while “nichenet” was utilized to derive the regulatory links between ligands and immune cell and stromal cell target genes (17, 18).

A total of 52 control and 492 PRAD patient data were subjected to a difference analysis utilizing the limma package. The criteria of P < 0.05 and Log2FC>1 were utilized to identify genes as DEG. The visualization of the heat maps and volcano plots of DEGs was carried out utilizing the ggplot2 and pheatmap packages, respectively. The DEGs were subsequently analyzed by KEGG and GO using the online instrument Metascope (19).

The “goodSamplesGenes” function, a component of the R “WGCNA”, was employed to assess the necessity of gene filtering within a sample and to determine an optimal soft threshold (20). The analysis involved the establishment of a co-expression network to examine the link between the genes associated with the module (ME) and PRAD.

The identification of candidate genes involved the identification of pairwise intersecting genes from three sources: marker genes of core cells, DEGs associated with PRAD, and module genes. These candidate genes were subsequently used in the least absolute shrinkage and selection operator (LASSO) regression analysis, which was executed using the R “glmnet”. The aim of the LASSO regression was to minimize the number of genes in the final risk model, thereby screening for characteristic genes associated with prognosis (21). In order to decrease the number of genes in the final risk model, variables with P-values 0.05 were subjected to the LASSO regression analysis, which was executed utilizing the “glmnet” R package. The prognosis-predictive model was developed employing the subsequent formula: risk score = gene exp11+gene exp22+…+gene expression (where “gene expression” represents the gene expression value and represents the corresponding LASSO regression coefficient). The R programs “survminer” and “ggrisk” were employed to visualize the survival curves and risk maps of patients. Additionally, the “survROC” program was utilized to plot ROC curves, assessing the performance of risk scores in predicting overall survival (OS) at 1, 3, and 5 years among individuals with PRAD.

Univariate analysis was performed to examine the predictive significance of risk models and clinical parameters (age, T stage, M stage, N stage, riskScore). In contrast, multivariable Cox analysis of OS was utilized to determine independent risk factors and predict new events in PRAD. Utilizing the “cph” function in R, a nomogram model was plotted to visualize this predictive model and predict potential new events for patients at 1, 3, and 5 years. The validity of the bar chart was validated through regression analysis (22).

To investigate functional variations and related pathways between the high-risk and low-risk groups, a GSEA was performed utilizing the clusterProfler package on all genes in the TCGA samples from both groups (23). A collection of 50 human cancer marker pathway genes was obtained from the Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) (24). Subsequently, GSVA was performed on all genes from samples in both groups, and the variations in GSVA scores between both risk samples were assessed utilizing the “limma” (25, 26).

The ssGSEA, implemented through the R “gsva”, was performed for all samples within the TCGA-PRAD dataset to derive profiles of 28 immune cell infiltrations (25). This was done to facilitate the comparison of immune cell characteristics between both risk groups. Pearson’s correlation coefficient was employed to examine the link between risk score and immune infiltrating cells. The R “maftools” was utilized to analyze the mutational variations between both risk groups (27).

The half-maximal inhibitory concentration (IC50) of chemotherapeutic agents was estimated using the pRRophetic R package, which is widely applied for drug sensitivity prediction in cancer research. The algorithm constructs a ridge regression model trained on drug response profiles from the Genomics of Drug Sensitivity in Cancer (GDSC) database and predicts the IC50 values of each drug for individual TCGA-PRAD samples based on their bulk RNA-seq expression profiles. All parameters were set to default, and batch correction was performed as recommended in the original pipeline to ensure prediction robustness.

Human tissue samples were collected at Shenzhen Hospital, Southern Medical University under approval from the Institutional Ethics Committee (Approval No.: NYSZYYEC20210011). Written informed consent was obtained from all participants. All samples used in this study were fresh-frozen tissues, which were snap-frozen immediately after collection and stored at –80 °C until further processing.

Total RNA was extracted using RNAiso Plus (TaKaRa), and reverse transcription was performed with an RNA reverse transcription kit (TaKaRa). Actin served as the internal control, and gene expression levels were quantified using the ΔΔCt method. The results were expressed as fold changes compared to the PBS-treated control group. The detailed primer sequences have been included in the Supplementary Materials and can be found in the corresponding table (Supplementary Table 1).

The DU145 human prostate cancer cell line was obtained from the American Type Culture Collection (ATCC) in February 2025. Cell line authentication was performed prior to experimentation using short tandem repeat (STR) profiling, confirming a match with the reference STR profile and excluding cross-contamination. Mycoplasma contamination was tested using the MycoAlert™ Mycoplasma Detection Kit (Lonza), and all results were negative.

Cells were seeded in 6-well plates at a density of 2.0–3.0 × 10⁵ cells/well and cultured overnight to reach 60–70% confluence. For plasmid transfection, cells were transfected with 1.5–2.0 μg DNA using Lipofectamine 3000 (Thermo Fisher) according to the manufacturer’s protocol. For siRNA experiments, cells were transfected with 10–30 nM siRNA using Lipofectamine RNAiMAX (Thermo Fisher). All transfections were performed in antibiotic-free medium, followed by replacement with complete medium 24 h later. Cells were collected 24–48 h after transfection for RNA and protein analyses.

Targeted silencing of GDPD3 in DU145 cells was achieved using siGDPD3, while a scrambled siRNA (siNC) was used as a negative control. Cell viability and proliferation were assessed using the Cell Counting Kit-8 (CCK8, Beyotime Biotechnology). Cells were seeded in 96-well plates at a density of 1 × 10^4 cells per well and allowed to adhere for 24 hours at 37°C in a 5% CO2 incubator. After treatment, 10 µL of CCK8 reagent was added to each well, followed by an additional incubation for one to four hours at 37°C (28). The absorbance at 450 nm was measured using a microplate reader (BioTek Instruments) to determine cell viability. To monitor cell growth over time, the assay was conducted at multiple time points (1, 2, 3, and 4 days), and relative cell viability was calculated by normalizing the absorbance values to those of the control group.

The wound healing assay was employed to evaluate cell migration. DU145 cells were transfected with shRNA targeting GDPD3 (shGDPD3) to suppress its expression; shNC served as the control. Cells were seeded onto 6-well plates at a density of 1 × 10^6 cells per well and cultured until they reached 90% confluence. Cells were pretreated with mitomycin-C (10 µg/mL, 2 h) to inhibit proliferation. A linear wound was created in the cell monolayer using a sterile 200 µL pipette tip, and any detached cells were removed by washing the wells twice with PBS. Fresh media, either with or without treatment, was then added, and the cells were incubated at 37°C in a 5% CO2 incubator. Images of the wound area were captured at 0 and 24 hours using a light microscope (×40 magnification). The wound closure rate was calculated by measuring the wound width multiple times using ImageJ software.

In DU145 cells, GDPD3 expression was knocked down using shGDPD3, with a non-targeting shRNA (shNC) used as the negative control. The Transwell migration assay was conducted using a Transwell chamber (Corning, USA) with an 8 µm pore size to measure cell migration ability (29). Cells were seeded in the top chamber at a density of 5 × 10^5 cells per well in 200 µL of serum-free media. The lower chamber contained 600 µL of complete medium with 10% FBS as a chemoattractant. After a 24-hour incubation, non-migratory cells on the top surface of the membrane were removed with a cotton swab. The migrated cells on the bottom surface were fixed with 4% paraformaldehyde for 15 minutes and stained with crystal violet for 30 minutes. Images of the stained cells were taken using a light microscope (×40 magnification), and the number of migratory cells in five random fields was counted.

LPA levels were quantified using a commercial ELISA kit (Invitrogen, EEL166) following the manufacturer’s protocol. Briefly, cell culture supernatants were collected, clarified by centrifugation, and incubated in LPA-specific antibody-coated wells. Absorbance at 450 nm was measured using a microplate reader, and LPA concentrations were calculated from a standard curve.

Cells were fixed with 4% paraformaldehyde for 15 min, permeabilized using 0.1% Triton X-100 for 10 min, and blocked with 5% BSA for 1 h at room temperature. Cells were then incubated overnight at 4°C with the anti-GDPD3 primary antibody (Abcam, ab214375), followed by incubation with fluorophore-conjugated secondary antibodies for 1 h at room temperature in the dark. Nuclei were counterstained with DAPI, and fluorescence images were acquired using a Leica fluorescence microscope.

Proteins were extracted from cells or tissues using lysis buffer containing PMSF, and the concentration was determined using a BCA kit. Protein samples (20μg per lane) were mixed with 5× loading buffer and boiled for 10 min. Proteins were separated on 10% SDS–polyacrylamide gels, followed by transfer onto PVDF membranes (Millipore). The membrane was blocked with 5% skim milk in TBST for 2 hours, then incubated with primary antibody (1:1000) overnight and secondary antibody (1:5000) for 1 hour, with three washes in TBST between incubations. Finally, the proteins were detected using chemiluminescent substrate (ECL) and imaged with a chemiluminescence imaging system. Primary antibodies, including E-cadherin (Proteintech, 20874-1-AP), Vimentin (Proteintech, 10366-1-AP), α-SMA (Proteintech, 14395-1-AP), LPAR1 (Proteintech, 86655-2-RR), AKT (Proteintech, 10176-2-AP), phospho-AKT (p-AKT; Proteintech, 66444-1-Ig), and β-Actin (Proteintech, 66009-1-Ig) as a loading control.

R software (version 4.2.1) was used for all further data processing and statistical analyses. The methodologies for the bioinformatics assessments are detailed in the relevant subsections. Statistical significance was denoted as follows: * p < 0.05, ** p < 0.01, and *** p < 0.001. A p-value below 0.05 indicates significance.

The fundamental framework of the current research is depicted in Supplementary Figure S1. Initially, a collection of scRNA-seq samples was obtained from a database consisting of 32 patients diagnosed with PRAD. Subsequently, a rigorous screening process was conducted to eliminate cells that did not meet the predetermined criteria. As a result, a total of 162,897 cells were successfully retrieved and deemed suitable for subsequent investigations. The annotation of cell clusters was executed using the CellMarker database (30), and relevant references, generating seven cell clusters (Figure 1A). The cell types present include fibroblasts, endothelial cells, epithelial cells, monocyte/macrophages, B cells, T cells, and mast cells. Bubble charts were utilized to visually represent the expression of critical marker genes for each cell type (Figure 1B). Subsequently, the feature genes associated with cell types were subjected to GO and KEGG enrichment analyses. The majority of these pathways were enriched in tumor- and immune-related pathways (Supplementary Figures S2A–G). Then, we performed PPI network analysis on the feature genes of these cell types to identify key genes associated with PARD (Figure 1C).

Furthermore, Pseudo-time trajectory analysis using Monocle 2 revealed three distinct differentiation pathways in PRAD, demonstrating a gradual progression from progenitor states. Notably, immune and non-immune cells exhibited temporally divergent differentiation patterns: non-immune lineages predominantly underwent early-stage specification, while immune cell maturation was biased toward later trajectory phases (Supplementary Figure S3A). The usage of ligand-receptors is expected to facilitate intercellular communication. The data showed that fibroblasts, endothelial cells, and epithelial cells communicated more often (Supplementary Figure S3B). Moreover, it was observed that the MIF cell signaling pathway exhibited significant activity in facilitating intercellular communication. The major ligand pairs involved in the MIF signaling pathway are CD74-CXCR4 and CD74-CD44 (Supplementary Figure S3C). Non-immune cells are primarily responsible for generating MIF signals, while immune cells consist of signal receptors and influences (Supplementary Figures S3D, E). Furthermore, our analysis of the regulatory network for ligands and target genes indicated that 28 ligands were selectively activated in immune cells (Supplementary Figure S3F). A heat map was generated to visualize the distribution of 28 ligands on immune cells and their impact on the regulation of target genes in non-immune cells. (Supplementary Figure S3G).

In summary, we identified the characteristic genes of different cell types and highlighted their functions.

In the TCGA-PRAD database, we analyzed the differentially expressed genes (DEGs) between the tumor and normal groups. The corresponding clinical information used in this analysis has been provided in the Supplementary Materials (Supplementary Table 2). A total of 1,219 significant DEGs were identified. The distribution of DEGs was visually represented using volcano plots and heatmaps (Figures 2A, B). Then, we conducted KEGG and GO enrichment analyses of DEGs. The KEGG analysis revealed that these DEGs were primarily enriched in the neuroactive ligand-receptor interaction pathway (Figure 2C). The GO analysis showed that the DEGs were mainly enriched in pathways associated with various tumor-related biological processes (Figures 2C–F).

WGCNA was utilized to identify participants in the progression and development of PRAD. While constructing the co-expression network, the power was set to 5 when the fit index of the soft-threshold scale-free topology attained 0.80 (Figures 3A, B). In order to integrate a similar module tree algorithm through dynamic shear analysis, MEDissTres was set to 0.2. Following the merging process, ten modules were obtained (Figure 3C), and a correlation heatmap depicting the relationships between these modules was generated (Figure 3D). As per the correlation coefficient and P-value, MEbrown and MEblue (consisting of 2545 genes) were chosen as the key modules (Figure 3E). A scatter plot was generated to illustrate the clinical relevance of the brown and blue modules (Figure 3F).

Using Venn diagrams, we illustrated the intersection of feature genes from single-cell sequencing, WGCNA module genes, and DEGs of cell Bulk RNA-seq. A total of 145 overlapping genes were identified (Figure 4A). Then, we performed Lasso regression to select the key genes (Figures 4B, C). Subsequently, a riskScore model was constructed based on the expression of 21 risk-associated genes (Table 2). The Kaplan-Meier indicated a statistically significant difference in OS between individuals with high-risk and low-risk scores (Figures 4D, E). To enhance the credibility of the risk model, an assessment was conducted by calculating the Receiver Operating Characteristic (ROC) curve for OS. The Area Under the Curve (AUC) values at 1, 3, and 5 years were found to be larger than 0.80, suggesting a higher level of effectiveness of the risk model (Figure 4F). In conclusion, our prognostic model demonstrated outstanding predictive accuracy for PRAD. To evaluate the robustness of our risk model, we validated it using an external dataset from the ICGC database. The Kaplan–Meier survival analysis showed that patients in the high-risk score group had a significantly higher mortality rate, which is consistent with the results of our risk model (Supplementary Figure S7A). Furthermore, we analyzed the prognostic relevance of the 21 risk-related genes in PRAD and found that the majority of these genes were negatively correlated with patient prognosis (Supplementary Figure S7B). These findings support the predictive accuracy of the risk score model based on these genes in PRAD.

The study employed univariate and multivariate Cox analyses to determine independent predictive factors to assess clinical features and risk scores. RiskScore demonstrated superior predictive capability compared to the T, M, and N stages. (Figures 5A, B). The nomogram model (Figure 5C) incorporated the independent prognostic variables. Furthermore, the calibration curve demonstrated the substantial predictive efficacy of the model (Figure 5D). The accuracy of 1-, 3-, and 5-year predictions for was assessed using the Net benefit strategy (Figures 5E–G). Consequently, the findings of the study revealed that risk score served as an independent prognosis-predictive factor and that the nomogram exhibited a high predictive capability for DFS in individuals with PRAD.

To explore the influence of both risk groups on cancer advancement, GSEA was conducted to identify the most significant enrichment pathways distinguishing between the two. The findings indicated that the high-risk group exhibited significant enrichment in immunological processes, including TNF-α and Interferon-α pathways, whereas the low-risk group exhibited enrichment primarily in androgen and acid metabolism pathways (Supplementary Figures S4A–C). Furthermore, GSVA was performed on all genes in both risk groups. The findings indicated that the high-risk group exhibited significant expression in processes related to cell communication, cyclase regulation, and cyclic nucleotide phosphodiesterase. Moreover, the low-risk group exhibited activation in the areas of tyrosine kinase, biosynthesis, and metabolism (Supplementary Figure S4D). It was discovered that missense mutations and SNPs dominated the PRAD patients (Supplementary Figure S4E). The mutation analysis outcomes between both risk groups revealed that the majority of mutation types in both groups were missense mutations. The high-risk group exhibited a larger percentage of mutations in comparison to the low-risk group (Supplementary Figures S4F, G). Overall, the findings of this study highlighted that predicted genes, pathways, and site mutations are relevant in PRAD.

The ssGSEA methodology was utilized to calculate infiltration scores for 28 distinct immune cell populations across various risk groups. The findings unveiled statistically significant differences in infiltration levels among 17 immune cell species (Supplementary Figure S5A). The heatmap illustrates the correlation between prognostic genes and immune cells (Supplementary Figure S5B). The CIBERSORT algorithm was used to assess the correlation between riskScore and tumor-infiltrating immune cells (Supplementary Figure S5C), while the ESTIMATE algorithm was employed to evaluate the immune and stromal components of the tumor microenvironment (Supplementary Figure S5D).

We assessed the drug sensitivity of 6 potential drugs in PRAD patients by comparing their IC50 values, and our findings revealed that the high-risk cohort demonstrated a greater responsiveness to MTX, 5-FU, ADM and DDP in comparison to the low-risk cohort. Interestingly, patients with low-riskScore displayed an elevated sensitivity to CTX (Supplementary Figures S6A–F). Our study indicated that riskScore critically involved in the regulation of immune checkpoints, including CD274, PDCD1LG2, and IGSF8, were found to be expressed at elevated levels within the high-risk cohort, while CTLA4, HAVCR2, LAG3, PDCD1, TIGIT, and ITPRIPL1 showed higher expression in the low-risk group. These findings may provide valuable insights for immunotherapy (Supplementary Figure S6G).

In order to validate the accuracy of our analysis, we collected 12 tumor tissues samples from PRAD patients. 12 normal samples were used as controls. First, we detected the mRNA expression of 21 risk genes and found that GDPD3, ITGAX, and FAM167B were significantly overexpressed in tumor tissues (Figure 6A). We further analyzed the expression patterns of these three genes at the single-cell level. GDPD3 and FAM167B were predominantly expressed in progenitor cells, whereas ITGAX was mainly expressed in monocyte–macrophage populations (Supplementary Figures S8A–I). Subsequently, we examined the protein expression of all three genes and found that GDPD3 was highly expressed in tumor tissues, with levels increasing progressively across clinical stages. In contrast, ITGAX and FAM167B showed no significant differences in expression. (Figure 6B). These findings were further confirmed through immunofluorescence analysis (Figure 6C). Subsequently, in vitro experiments were conducted using DU145 cells transfected with shRNAs targeting GDPD3, ITGAX, and FAM167B to assess their effects on cell proliferation, migration, and invasion. Knockdown of GDPD3 significantly reduced all three cellular functions, whereas silencing of ITGAX and FAM167B showed no significant effects (Figures 6D, E). To further validate the functional relevance of the candidate genes, we employed the LNCaP prostate cancer cell line as an additional model system. Among the genes examined, only GDPD3 knockdown produced a substantial phenotypic effect in LNCaP cells. Silencing GDPD3 significantly reduced cell proliferation and markedly impaired migration and invasion, as evidenced by the pronounced decrease in transwell‐traversing cells. In contrast, knockdown of the other genes did not result in observable changes in any of these phenotypes. These results confirm that GDPD3 is the only gene among those tested that is required to maintain the proliferative, migratory, and invasive capacities of LNCaP prostate cancer cells (Supplementary Figures S9A–E).These results identify GDPD3 as a key regulator of prostate cancer cell proliferation, migration, and invasion, whereas ITGAX and FAM167B appear to have limited functional impact.

Given that GDPD3 encodes a glycerophosphodiesterase involved in lysophospholipid metabolism (31), we investigated whether its upregulation influences LPA production. First, we analyzed the correlation between GDPD3 and LPA in the TCGA-PRAD dataset and found a strong positive association (Figure 7A). To further investigate, we overexpressed GDPD3 and observed an increase in LPA levels (Figure 7B). We then explored whether LPA could induce epithelial–mesenchymal transition. DU145 cells were treated with LPA for 24, 48, and 72 hours, and EMT induction was detected as early as 24 hours (Figure 7C). Next, we used the LPAR1 inhibitors Ki6425 and AM095 to block LPA signaling. Inhibition of LPAR1 led to decreased expression of EMT markers Vimentin and α-SMA (Figures 7D, E). LPA stimulation selectively increased the mRNA expression of Twist, thereby enhancing the EMT program, while the levels of Slug and Snail remained unchanged (Supplementary Figures S10A, B). Conversely, genetic deletion of LPAR1 significantly reduced Twist expression without affecting Slug or Snail, indicating that LPA promotes EMT primarily through an LPAR1–Twist axis(Supplementary Figures S10C, D). These findings suggest that GDPD3 promotes epithelial–mesenchymal transition in prostate cancer cells by elevating LPA levels, which activate downstream signaling through LPAR1.

To further validate the role of LPAR1 in LPA-mediated epithelial–mesenchymal transition (EMT), we performed siRNA-mediated knockdown of LPAR1 (Figures 8A, B). Silencing of LPAR1 markedly suppressed EMT, as evidenced by increased E-cadherin and decreased Vimentin expression (Figures 8C, D). Mechanistically, LPA treatment induced phosphorylation of AKT, whereas LPAR1 knockdown reduced this phosphorylation, indicating that LPA activates AKT signaling via LPAR1(Figures 8E, F). Moreover, treatment with an AKT pathway inhibitor attenuated LPA-induced EMT, further supporting the involvement of the LPA–LPAR1–AKT axis in regulating EMT(Figures 8G, H). These findings reinforce that GDPD3-driven LPA production promotes EMT in prostate cancer cells through LPAR1-dependent activation of AKT signaling.