- Department of Urology, Tangdu Hospital, Fourth Military Medical University, Xi’an, China
Background: Myeloid cell differentiation (MCD) has an important correlation with prostate cancer (PCa), but the mechanism of action of the former in the latter is still under investigation. This study designed to investigate the prognostic genes related to MCD in PCa and the associated mechanisms.
Methods: The related data were downloaded from public databases. Differentially expressed genes (DEGs) were intersected with MCD related genes (MCDRGs) to acquire candidate genes. Candidate prognostic genes with a causal relationship to PCa were further obtained through Mendelian randomization (MR). Prognostic genes were acquired by univariate Cox regression analysis and Least Absolute Shrinkage and Selection Operator (LASSO) analysis. Then, the risk model was built based on prognostic genes. Immune infiltration, nomogram model, and drug sensitivity were employed to investigate the roles of prognostic genes in PCa. The manifestation of prognostic genes in key cells was also investigated by single-cell sequencing (scRNA-seq) analysis. Finally, the manifestation of prognostic genes were authenticated by in vitro experiments.
Results: The 23 candidate prognostic genes had a causal relationship with PCa. The 5 prognostic genes (NR3C1, BMP2, RACGAP1, TLR3, FASN) were identified. The risk models suggested that high risk group (HRG)’s survival rate was inferior to that of low risk group (LRG). The nomogram indicated that prognostic genes could effectively predict the survival status of PCa patients. There were 18 immune cells that suggested notable differences between the HRG and the LRG. The HRG and LRG suggested notable differences in sensitivity to 86 drugs such as AZD8186. Epithelial cells were considered as key cells. Only FASN was consistently active during critical cell differentiation. The in vitro results were consistent with the results of bioinformatics analysis, indicating that the analysis results were reliable.
Conclusion: This study identified 5 prognostic genes and a risk model, suggesting a fresh thought on the subsequent development of PCa related drugs.
1 Introduction
Prostate cancer (PCa) represents one of the most prevalent malignancies in the male genitourinary system, globally ranking as the second most frequently diagnosed cancer among males (1, 2). Its pathogenesis is multifactorial, involving age, genetic predisposition, racial disparities, and dysregulated gene expression, with approximately 20% of cases exhibiting aberrations in DNA repair pathways (1, 3). Although localized PCa in early stages can be effectively managed through radical prostatectomy or radiotherapy, biochemical recurrence occurs in 20-40% of treated patients (1, 3). Notably, individuals progressing to metastatic castration-resistant PCa (mCRPC) demonstrate a 5-year survival rate below 30%, underscoring the critical need for advanced therapeutic strategies (2, 3).
Current clinical management strategies encompass androgen deprivation therapy (ADT), novel androgen receptor pathway inhibitors, chemotherapy, and immune checkpoint inhibitors (ICIs) (3–5). However, these therapeutic approaches demonstrate suboptimal efficacy against metastatic lesions and biochemical recurrence, coupled with a substantial toxicity burden that often limits their clinical utility (3–5).
In-depth investigation of the molecular mechanisms underlying PCa represents a critical avenue to overcome current therapeutic limitations. scRNA-seq technologies have unveiled substantial heterogeneity within cancer-associated fibroblasts (CAFs) in the tumor microenvironment, which promote oncogenic niche formation through secretion of specific cytokines, a biological feature positively correlated with tumor progression (6, 7). Genomic studies have identified key driver events including PTEN deletion and BRCA1/2 mutations that induce homologous recombination repair (HRR) pathway dysfunction (3, 8, 9). PARP inhibitors developed based on these molecular characteristics demonstrate significant clinical benefits in BRCA-mutant patients, extending median radiographic progression-free survival (rPFS) to 7.4 months (3). Prognostic models systematically integrating 10 CAFs core regulatory genes such as THBS1 and LDHA exhibit robust discriminative power for stratifying patient survival outcomes through risk scoring (3). These findings not only elucidate the evolutionary biology of PCa but also provide critical theoretical foundations for developing targeted therapies and advancing personalized medicine. Therefore, identification of novel prognostic biomarkers and establishment of precision prediction frameworks will be pivotal to resolving therapeutic challenges in advanced PCa.
Myeloid cells, as central orchestrators of the tumor immune regulatory network, critically determine immune evasion, tumor progression, and clinical outcomes through their differentiation states and functional plasticity (10, 11). Accumulating evidence demonstrates that myeloid-derived suppressor cells (MDSCs) in solid tumor microenvironments drive oncogenesis through dual mechanisms, direct suppression of CD8+ T cell antitumor activity via effector molecules including arginase-1 (ARG1), inducible nitric oxide synthase (iNOS), and reactive oxygen species (ROS), and facilitation of metastatic dissemination through pro-angiogenic mediators such as VEGF (10, 12). Notably, LOX-1 surface expression on MDSCs exhibits significant inverse correlations with circulating tumor DNA burden and overall survival rates (10, 12). Mechanistically, tumor-derived oxidized lipids potentiate myeloid immunosuppressive capacity by activating the STAT3 signaling axis and CSF1R pathway, thereby inducing metabolic reprogramming that sustains immune tolerance (11, 12).
Previous study has indicated that PCa cells elicit functional reprogramming of myeloid lineages (including THP-1 and HL-60) through stress protein secretion, particularly heat shock protein 27 (Hsp27), manifesting as surface marker polarization and VEGF secretion dysregulation, thereby modulating tumor-immune crosstalk (13). Notably, therapeutic interventions targeting myeloid surface receptors like the CD47-SIRPα axis enhance macrophage-mediated tumor phagocytosis compared to controls and improve anti-tumor immune responses in preclinical studies (11). However, the precise molecular targets governing myeloid cell differentiation (MCD) in PCa progression remain poorly characterized, particularly key signaling pathways such as JAK/STAT or NF-κB signaling pathways and immunoregulatory cytokine networks (14). Moreover, the bidirectional regulatory network between neoplastic cells and myeloid populations involving cytokine crosstalk and surface receptor interactions requires comprehensive investigation to delineate their co-evolution mechanisms.
Mendelian randomization (MR) leverages the random assortment of genetic variants during gametogenesis to emulate randomized controlled trials, effectively circumventing confounding biases and reverse causation inherent in conventional observational studies (15, 16). ScRNA-seq offers high-throughput resolution of whole transcriptomes at individual cell resolution, enabling precise mapping of immune cell heterogeneity and functionally distinct subsets within tumor microenvironments. Pioneering studies have employed scRNA-seq to construct pan-cancer endothelial cell atlases, uncovering tumor-specific endothelial subpopulations orchestrating pro-angiogenic and immunosuppressive programs (5, 17, 18). In PCa pathobiology investigations, Miao et al. integrated scRNA-seq with bulk transcriptomics to identify MXRA8-mediated tumor progression through dysregulating oxidative stress pathways in prostate tumor niches (16), while Ye et al. established MR-based causal relationships between genetically proxied CD25+ naïve B cell abundance and PCa risk (19). Nevertheless, the synergistic application of MR framework with scRNA-seq technologies to decipher MCD-PCa crosstalk mechanisms remains substantially underexplored, representing a critical knowledge gap in the field.
In this study, we employed integrated multi-omics analyses (incorporating bulk transcriptomic and scRNA-seq data) to identify prognostically significant genes causally linked to MCD through MR framework, ultimately constructing a clinical-grade prognostic signature. Bioinformatics interrogation systematically delineated the molecular regulatory circuitry underlying these candidate genes in PCa progression. Furthermore, single-cell resolution analysis uncovered their expression dynamics within disease-associated cell subpopulations, while reverse transcription-quantitative polymerase chain reaction (RT-qPCR) experiments provided preliminary validation of their potential roles in tumorigenesis. To further validate the functional roles of candidate genes in PCa progression, we performed complementary in vitro assays: Western blot analysis to examine protein expression, Ki67 immunofluorescence staining to assess cell proliferation, and scratch-wound assays to evaluate cell migratory capacity, aiming to verify their regulatory effects on PCa cell malignant behaviors. Our findings provide multi-dimensional evidence elucidating MCD-associated molecular networks in PCa, establishing both conceptual and experimental foundations for developing personalized therapeutic strategies and targeted drug discovery.
2 Method
2.1 Data collection
The transcriptome data (TCGA-PRAD) on gene expression matrix and relapse information of PCa were extracted from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). The TCGA-PRAD dataset (access time: November 20th, 2024) included 502 PCa tissue samples (397 samples with relapse information) and 52 paracancer (control) tissue samples. Among the 397 samples, 70% samples (278 samples) were employed as the training set (TCGA-PRAD-train), and 30% samples (119 samples) were employed as the validation set (TCGA-PRAD-validation). The single-cell RNA sequencing (scRNA-seq) data (GSE141445) of PCa was scoured from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE141445 (GPL24676, access time: November 20th, 2024) included 13 PCa tissue samples. Validation dataset GSE116918 contained 248 samples with complete recurrence information. The MR data (EBI-A-GCST90018905) of PCa and eQTL data were acquired from the Integrative Epidemiology Unit (IEU) Open Genome-wide Association Study (GWAS) database (https://gwas.mrcieu.ac.uk/). The “EBI-A-GCST90018905” included 24,119,306 SNPs from 211,227 Europeans (case: control=11,599: 199,628). The 423 MCD-related genes (MCDRGs) used in the study were obtained by merging and removing duplicates from three myeloid cell differentiation-related datasets: GOBP_NEGATIVE_REGULATION_OF_MYELOID_CELL_DIFFERENTIATION (95 genes), GOBP_POSITIVE_REGULATION_OF_MYELOID_CELL_DIFFERENTIATION (104 genes), and GOBP_MYELOID_CELL_DIFFERENTIATION (423 genes), which were downloaded from the MCD related genes (MCDRGs) were acquired from Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/human/search.jsp) and relevant references (20) (Supplementary Table 1).
2.2 Differential expression analysis
To acquire differentially expressed genes (DEGs) between PCa and control samples in TCGA-PRAD, “DESeq2” package (v 3.4.1) was carried out (PCa vs control) (|log2Fold Change (FC)| > 0.5, P < 0.05) (21). On the basis of the log2FC value, DEGs were visualized and the top 10 up/down-regulated gene names were labeled by the volcano plot utilizing “ggplot2” package (v 3.4.1) (22). Similarly, the expressions of the top 10 up/down-regulated genes between PCa and control groups were displayed by the heat plot utilizing “ComplexHeatmap” package (v 2.14.0) (23).
2.3 Identification and functions of candidate genes
DEGs were intersected with MCDRGs to acquire candidate genes via “ggvenn” package (v 1.7.3) (24). To probe the organic activities and signal pathways involved in candidate genes, based on background gene set of “org.Hs.eg.db” package (v 3.16.0) (25), Gene Ontology (GO) (P < 0.05) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (P < 0.05) enrichment analyses were carried out utilizing “clusterProfiler” package (v 4.7.1.3) (26). Subsequently, to determine the interactions of proteins encoded by candidate genes, candidate genes were uploaded to STRING database (https://cn.string-db.org/) (interaction score > 0.4). The results were imported into Cytoscape software (v 3.9.1) and Protein-Protein Interaction (PPI) network was constructed (27).
2.4 MR analysis
To probe the causal dependence between candidate genes and PCa and acquire candidate prognostic genes, MR analysis was performed with PCa as outcome event and candidate genes as exposure factors utilizing “TwoSampleMR” package (v 0.6.1) (28). The MR consisted of 3 assumptions: (1) instrumental variables (IVs) were linked with exposure factors; (2) IVs could only influence outcomes by exposure factors; (3) IVs were not linked with potential confounders. To obtain effective IVs, exposure factors and outcome event were read and IVs were screened via “extract_instruments” function. The screening criteria were as follows: (1) IVs strikingly linked with exposure factors (P < 5×10−6); (2) IVs that exhibited linkage disequilibrium were removed with R2 = 0.001, kb=10, and clump=TRUE; (3) IVs that were strikingly linked with outcome were removed with rsq=0.8 and proxies=TRUE; (4) IVs whose F statistic < 10 were removed (, R2 represented the cumulative explanatory variance of SNPs, N represented the number of the samples, and K represented the number of SNPs). Then, the effect alleles and effect sizes were unified via “harmonCIe_data” function. The 5 algorithms of the “mr” function were utilized to conduct MR analysis for each exposure factor and outcome, which included Inverse Variance Weighted (IVW) (29), Weighted Mode (30), MR Egger (31), Simple Mode (28), and Weighted Median (32). MR analysis mainly relied on IVW results. Exposure factors with SNP > 2 and P < 0.05 were considered as the exposure factors that had a causal relationship with PCa. The odds ratio (OR) > 1 suggested that the exposure factor was a contributing factor to the risk of PCa, while an OR < 1 indicated that it was a protective factor. Notably, the scatter plot was drawn to further identify the correlation between exposure factors and outcome in combination with SNP-exposure effects and SNP-outcome effects via “mr_scatter_plot” function. The forest plot was drawn to evaluate the diagnostic power of the estimated exposure factors of each SNP site on the outcome via “mr_forest_plot” function. The funnel plot was drawn to judge whether the analysis was random and adhered to Mendel’s second law random grouping in combination with the β and standard error (SE) of each IV via “mr_funnel_plot” function.
Moreover, to evaluate the reliability of MR analysis results, sensitivity analysis was applied. Heterogeneity test was performed via “mr_heterogeneity” function (P > 0.05, Cochran’s Q test). Horizontal pleiotropy test was performed to find out whether confounding factors existed via “mr_pleiotropy_test” and “MR-Egger” functions (P > 0.05). To observe whether the SNPs of each IV caused considerable changes in the outcome, Leave-one-out (LOO) analysis was performed via “mr_leaveoneout” function.
To verify that the results of the forward analysis were not interfered by reverse causality and determine the validity of the causal sequence between the outcome and the exposure factors, Steiger test was carried out via “steiger_filtering” function (Steiger-dir=TRUE, P < 0.05).
At last, the candidate genes examined by sensitivity analysis and Steiger test were specified as candidate prognostic genes.
2.5 Identification of prognostic genes and construction of prognostic models
To identify prognostic genes, in TCGA-PRAD-train, univariate Cox regression analysis was performed on the basis of candidate prognostic genes via “survival” package (v 3.5-3) (33)(P < 0.2) (34, 35) and the result was presented via “forestplot” package (v 2.0.1) (36). Genes with consistent hazard ratios (HRs) or ORs in both univariate Cox regression analysis and MR analysis were utilized for proportional hazards (PH) assumption test via “cox.zph” function (P > 0.05). After that, Least Absolute Shrinkage and Selection Operator (LASSO) was implemented utilizing “glmnet” package (v 4.1.4) (37)(10-fold cross-validation). Finally, genes that had passed through the above analyses in sequence were defined as prognostic genes to construct a risk model.
In TCGA-PRAD-train, on the basis of the expression of prognostic genes and the risk coefficients gained from LASSO regression, PCa patients were scored with the following formula.
Expr represented the expression level of each prognostic gene and coef signified the risk coefficient of each prognostic gene. Notably, according to the median of risk score, PCa patients were categorized into high risk group (HRG) and low risk group (LRG) and risk score distribution and survival status of patients were displayed. According to relapse time and relapse status of PCa patients, the “survminer” package (v 0.4.9) (38) was implemented to draw and compare survival curves of HRG and LRG (P < 0.05). The Area Under Curve (AUC) values of receiver operating characteristic (ROC) curves at 1, 2, and 3 years were employed to evaluate the accuracy of risk model utilizing “survivalROC” package (v 1.18.0) (39) (AUC > 0.6). The validation was performed in the GSE116918 dataset. Additionally, the heat plot was drawn to display the expression levels of prognostic genes in HRG and LRG. Using the same method, 2 risk models were constructed in TCGA-PRAD-validation and all samples with relapse information of TCGA-PRAD to validate the accuracy of the above-mentioned model.
2.6 Construction of nomogram model
To evaluate the ability of prognostic genes to predict PCa recurrence rates, the “regplot” package (v 1.1) (40) was employed to build the nomogram model for PCa in the TCGA-PRAD-train. Each prognostic gene was scored separately, and the scores were added together to obtain the total scores. The higher the total scores, the higher the recurrence rate of the patient. Calibration curves built via “rms” package (v 6.5.0) (41) and ROC curves built via “survivalROC” package (v 1.18.0) (AUC > 0.6) were applied to compute the effectiveness of the nomogram model in clinical prediction at 1, 2, and 3 years for PCa.
2.7 Immune infiltration analysis
To estimate the tumor purity and the status of stromal and immune cells in the malignant tumor tissues within the tumor microenvironment of PCa patients. The data were normalized using the R package estimate (v 1.0.13, https://R-Forge.R-project.org/projects/estimate/), and the stromal score, immune score, and ESTIMATE score were calculated using the ESTIMATE algorithm. Wilcoxon test was employed to contrast the scores mentioned above between HRG and LRG (P < 0.05). For further evaluation of the situation of immune cells in the development process of PCa, based on single sample Gene Set Enrichment Analysis (ssGSEA) algorithm, in the TCGA-PRAD-train, “GSVA” package (v 1.46.0) (42) was applied to compute the enrichment scores of 28 immune cells (43) in the HRG and LRG. Wilcoxon test was implemented to contrast differences in immune cell infiltration between HRG and LRG (P < 0.05). What’s more, to understand the link between prognostic genes and differential immune cells, in the TCGA-PRAD-train, “psych” package (v 2.2.9) (44) was performed (|correlation (cor)| > 0.3, P < 0.05).
2.8 Pathways and GeneMANIA analysis
To determine the biological pathways involved in the occurrence and development of PCa in the HRG and LRG, in the TCGA-PRAD-train, “DESeq2” package (v 3.4.1) was employed to perform differential expression analysis between HRG and LRG and log2FC values were calculated. The log2FC values were sorted in descending order. Then, Gene Set Enrichment Analysis (GSEA) was performed via “clusterProfiler” package (v 4.7.1.3) (|Normalized Enrichment Score (NES)| > 1, P < 0.05). The top 5 pathways with notable P-values were presented. Besides, GeneMANIA database (http://genemania.org) was applied to predict the genes related to the functions of prognostic genes and their involved functions.
2.9 Drug sensitivity analysis
To probe the drug sensitivity of the HRG and LRG, drugs related to tumors were obtained from Genomics of Drug Sensitivity in Cancer 2 (GDSC2) database (https://www.cancerrxgene.org/). In the TCGA-PRAD-train, “pRRophetic” package (v 0.5) (45) was employed to determine the half maximal inhibitory concentration (IC50) of each tumor sample. Wilcoxon test was employed to contrast the differences in drug sensitivity between the HRG and LRG (P < 0.05). The result was displayed via “ggplot2” package (v 3.4.1).
2.10 Construction of molecular regulatory network
To explore the upstream regulatory factors of prognostic genes and their interaction relationships, based on NetworkAnalyst online website (https://www.networkanalyst.ca/NetworkAnalyst/uploads/ListUploadView.xhtml), Transcription Factors (TFs) were predicted by TRRUST database (https://www.grnpedia.org/trrust/) and miRNAs were predicted by miRTarBase database (v 9.0) (https://mirtarbase.cuhk.edu.cn/) and TarBase (v 9.0) (https://dianalab.e-ce.uth.gr/tarbasev9). Finally, the TF-mRNA-miRNA network was constructed utilizing Cytoscape software (v 3.9.1).
2.11 ScRNA-seq analysis
All scRNA-seq analyses were performed via “Seurat” package (v 5.0.1) (46). To ensure the accuracy and reliability of scRNA-seq data, in the GSE141445, the data was filtered via “PercentageFeatureSet” function. The screening criteria were as follows: (1) the number of genes in the cells ranged from 200 to 3,000; (2) genes with an expression level ranging from 200 to 4,531 and covered by at least 3 cells were retained. Then, the data after filtering were normalized via “NormalizeData” function. The 2,000 highly variable genes (HVGs) were acquired and the top 10 most mutated genes were labeled via “FindVariableFeatures” function and the result was presented via “LablePoints” function. Moreover, the samples of GSE141445 were subjected to normalization processing via “Scale Data” function. The HVGs were subjected to principal component analysis (PCA) via “runPCA” function. The P-value when the principal components (PCs) ranged from 1 to 30 was calculated via “Jackstraw” function. The values of the sudden variance drops when different values were taken for the PCs were calculated and the result was displayed via “Elbowplot” function. When P < 0.05, the PCs at the inflection point in the variance elbow plot were used for subsequent analysis. Additionally, unsupervised clustering analysis was applied via “FindClusters” and “FindNeighbors” functions (resolution=0.1). Notably, cells were clustered and result was visualized via “RunTSNE” function. Furthermore, to determine the cell types of the cell clusters, cell clusters were commented on the basis of marker genes (47). The expression patterns of marker genes in different cell clusters and cells were displayed through bubble plots.
2.12 Identification of key cells, cell communication, and pseudo-time analysis
To gain key cells associated with PCa development, in the GSE141445, the manifestation of prognostic genes in cells were displayed through bubble plots. The cells with the highest gene expression were labeled as key cells. Besides, the distribution of prognostic genes in cells was also presented.
In the GSE141445, the “CellChat” package (v 1.6.1) (48) was employed to explore the communication networks among different cells. The pairing of ligands and receptors among cells was presented through bubble diagrams using “ggplot2” (v 3.4.1). Furthermore, to investigate the differentiation of key cells, key cells were subjected to secondary clustering using the same method as that in 2.11. Then, the developmental trajectories of key cells were analyzed via the “Monocle” package (v 2.22.0) (49). The “DDRTree” package (v 0.1.5) (50) was employed to draw cell trajectory diagram. Finally, the expression of prognostic genes during the differentiation process of key cells was also explored.
2.13 Prognostic genes expression and reverse transcription quantitative PCR
To clarify the manifestation of prognostic genes, in the PCa tissue and paracancer tissue samples of TCGA-PRAD, Wilcoxon test was employed to compare the differences in manifestation of prognostic genes between PCa samples and control samples (P < 0.05). Subsequently, to verify the accuracy of the above results, RT-qPCR was performed. A total of 5 pairs of tissue samples were obtained at Tangdu hospital, including 5 PCa and 5 paracancer (control). Informed consent was obtained from all patients for the use of PCa tissue samples in this study. The informed consent form needed to be signed and filled out by all participants, while the ethical approval agency was The Clinical Ethics Committee of Tangdu Hospital of Air Force Medical University (permission number: TDLL-KY-202405-18). Then, total RNA of 5 pairs of tissue samples was extracted by TRizol reagent (Vazyme, Nanjing, Jiangsu). The RNA concentrations were computed by NanoPhotometer N50. Subsequently, mRNA was converted to cDNA by Hifair® III 1st Strand cDNA Synthesis SuperMix for qPCR test kit (Yeasen, Shanghai). Finally, RT-qPCR was carried out. The primers of reaction reagents, reaction conditions and genes were arranged in Supplementary Table 2. The internal reference gene was GAPDH, which was employed to normalize the results. The expression levels of prognostic genes were calculated by 2-ΔΔCt. The results were calculated by GraphPad Prism software (v 5.0) (51).
2.14 Western blotting
Cells were lysed in RIPA buffer (P0013B, Beyotime) to extract proteins, and protein concentrations were determined using a BCA assay kit (P0012, Beyotime). Samples were separated by SDS-PAGE (12%, Invitrogen) and transferred to PVDF membranes, which were then blocked with BSA for 2.5 h. After washing, membranes were incubated overnight at 4°C with primary antibodies, including anti-BMP2 (ab284387, 1:1000, Abcam) and anti-FASN (ab128870, 1:1000, Abcam). The next day, membranes were washed three times and incubated with HRP-conjugated goat anti-rabbit IgG secondary antibody (G-21234, 1:5000, Invitrogen) at room temperature for 2 h. Chemiluminescent substrate (ECL) was applied evenly to the membranes, and signals were captured using a ChemiDoc XRS+ imaging system (Bio-Rad, USA). Band intensities were analyzed with ImageJ, and protein levels were quantified as the ratio of the target band intensity to that of β-actin (ab272085, 1:1000, Abcam).
2.15 Ki67 staining
Vector, OE-BMP2, sh-NC, and sh-FASN cells were seeded and allowed to adhere for 24 h. Cells were washed twice with phosphate-buffered saline (PBS) and fixed with 4% paraformaldehyde for 20 min. Permeabilization was performed with 0.3% Triton X-100 for 15 min, followed by blocking with 6% donkey serum for 30 min. Cells were incubated overnight at 4°C with primary antibody against Ki67 (1:200 in PBS). The next day, cells were incubated with a fluorescent secondary antibody (1:500 in PBS) for 1 h at 37°C in the dark. After three PBS washes, nuclei were counterstained with DAPI for 10 min. Ki67 expression was assessed using a fluorescence microscope.
2.16 Scratch-wound assay
Vector, OE-BMP2, sh-NC, and sh-FASN cells (1×10^6 cells/well) were seeded in 6-well plates and allowed to adhere for 24 h. Once confluence exceeded 90%, a straight scratch perpendicular to the horizontal axis was made using a sterile 20-µL pipette tip. To remove debris, cells were rinsed with serum-free DMEM. Plates were incubated, and scratch closure was monitored at 24 h and 48 h after wounding under a high-power microscope. Migrating cells within the scratch area were quantified using ImageJ.
2.17 Cell culture and viral preparation
Human prostate cancer PC3 cells (ATCC® CRL-1435™) were maintained in RPMI-1640 medium (Gibco) supplemented with 10% fetal bovine serum (FBS; Corning) and 1% penicillin/streptomycin (HyClone) at 37°C with 5% CO2. Lentiviral particles expressing shRNA targeting human fatty acid synthase (FASN; target sequence: shRNA1:5′-GCATGGAGCGTATCTGTGAGAACTCGAGTTCTCACAGATACGCTCCATGTTTTTT-3′ shRNA2: 5′-GCTACGACTACGGCCCTCATTCTCGAGAATGAGGGCCGTAGTCGTAGCTTTTTT-3′) or non-targeting scramble control (shRNA-NC: 5′-GCACCCAGTCCGCCCTGAGCAAATTCAAGAGATTTGCTCAGGGCGGACTGGGTGCTTTTT-3′) in pLKO.1 vector were packaged in HEK293T cells using psPAX2 and pMD2.G plasmids. Viral supernatants were concentrated by ultracentrifugation (26,000 × g, 2.5 h) and titers determined via qPCR.
The human BMP2 ORF (NM_001200) was cloned into pLV-CMV-3flag-zsgreen vector via NotI and NheI site. Lentivirus production and titration followed procedures (titer ≥5×107 TU/mL). PC3 cells were transduced at MOI=20 with polybrene supplementation.
Cells were seeded in 6-well plates (5×105 cells/well) and incubated for 24 h to reach 50–60% confluency. Transduction was performed by replacing medium with viral suspension (MOI=20) containing 8 μg/mL polybrene (Sigma).
2.18 Statistical analysis
R software (v 3.9.1) was implemented to apply bioinformatics analyses. Wilcoxon test and t test were utilized to compare the disparities between 2 groups (P < 0.05).
3 Results
3.1 Identification and functions of candidate genes in PCa
By differential expression analysis, 13,091 DEGs were acquired, comprising 7,935 DEGs with up-regulated expression and 5,156 DEGs with down-regulated expression. All DEGs and the top 10 up/down-regulated gene names were presented (Figure 1A). Besides, the expressions of the top 10 up/down-regulated genes in PCa and control groups were also displayed (Figure 1B). After that, 141 candidate genes were acquired (Figure 1C). GO analysis indicated that candidate genes were enriched in 1,758 functions, including 1,619 biological process (BP) items such as MCD, 37 cellular component (CC) items such as chromosome centromeric core domain, and 102 molecular function (MF) items such as cytokine receptor binding (Figure 1D, Supplementary Table 3). KEGG analysis suggested that candidate genes were enriched in 85 pathways such as osteoclast differentiation (Figure 1E, Supplementary Table 4). Subsequently, PPI network indicated that there were interactions among proteins encoded by 122 candidate genes such as MYC, H4C11, and MMP9 (Figure 1F).

Figure 1. Multi-dimensional characterization of differentially expressed genes (DEGs) and functional annotation of candidate genes in pca. (A) Volcano plot of DEGs in pca. Orange dots represent upregulated genes, glaucous dots represent downregulated genes. (B) Heatmap analysis of top differentially expressed genes in PCa versus control groups. (C) Venn map to obtain a total of 141 candidate genes shared by DEGs and MCDRGs. (D) Functional enrichment of 141 candidate genes across GO categories. (E) Functional enrichment of candidate genes in KEGG pathways. (F) PPI network constructed using the 122 candidate genes.
3.2 Candidate prognostic genes with a causal relationship with PCa
In order to investigate the causal relationship between DE-MCDRGs and PCa outcomes, MR analysis was conducted using 141 candidate genes as exposure factors and PCa as the outcome. Firstly, the 46 candidate genes had a substantially causal relationship with PCa based on the IVW method (P < 0.05), including 23 risk factors (OR > 1) and 23 protective factors (OR < 1) (Table 1). The positive slope of the scatter plot indicated that the gene was a risk factor, while a negative slope indicated that it was a protective factor. The intercept close to 0 suggested the absence of confounding factors (Supplementary Figure 1). For protective factors, the effect size of each SNP was less than 0. In contrast, for risk factors, the effect size of each SNP was greater than 0 (Supplementary Figure 2). The symmetrical arrangement of SNPs around each exposure indicated that MR affirmed Mendel’s second law of randomization (Supplementary Figure 3). So 46 candidate genes were used for subsequent analysis. Notably, among the 46 candidate genes, heterogeneity was calculated for 35 of them. The P-values of 34 candidate genes were greater than 0.05, indicating the absence of heterogeneity. Since the P-value of SLC4A1 was less than 0.05, SLC4A1 was analyzed with random effects IVW (Supplementary Table 5). Meanwhile, among the 46 candidate genes, the P-values of 33 candidate genes were greater than 0.05, suggesting the lack of horizontal pleiotropy and the reliability of the results. However, the P-values of the remaining 13 candidate genes were less than 0.05, suggesting the presence of horizontal pleiotropy. To ensure the reliability of the results, these 13 candidate genes were excluded (Supplementary Table 6). LOO analysis demonstrated that the removal of any SNP had minimal impact on the results (Supplementary Figure 4). Steiger test indicated that only the P-values of 33 candidate genes were all less than 0.05 and the Steiger-dir values were all TRUE (Supplementary Table 7). Finally, taking all of the above screening results into account, only 23 candidate genes could be labeled as candidate prognostic genes.
3.3 The value of prognostic genes in risk models
In order to further screen for genes related to prognosis among the 23 candidate prognostic genes, univariate Cox regression analysis was conducted. The 9 candidate prognostic genes associated with PCa recurrence were gained utilizing univariate Cox regression analysis (P < 0.2), including 4 risk factors (HR > 1) and 5 protective factors (HR < 1) (Figure 2A). Notably, compared with the results of the MR analysis, only the OR/HR of 6 candidate prognostic genes was consistent. Besides, these 6 candidate prognostic genes were also identified by PH assumption test (P > 0.05) (Figure 2B). Among them, 5 candidate prognostic genes were identified by LASSO analysis (lambda.min=0.005517655) (Figure 2C). So NR3C1, BMP2, RACGAP1, TLR3, and FASN were considered as prognostic genes to build the risk model in TCGA-PRAD-train. On the basis of the median of the risk scores (1.111923), PCa patients were divided into a HRG (139 PCa patients) and a LRG (139 PCa patients). Patients’ risk score distribution (Figure 2D) and survival status (Figure 2E) were displayed. The survival curves indicated that the survival rate of the HRG was substantially lower than that of the LRG (P=0.00077) (Figure 2F). The AUC values of the ROC curves were all greater than 0.7, indicating that the risk model had good performance (Figure 2G). In addition, BMP2, RACGAP1, and FASN were prominently expressed in the HRG, while TLR3 and NR3C1 were prominently expressed in the LRG (Figure 2H).

Figure 2. Development and validation of a prognostic gene signature for PCa recurrence. (A) Univariate Cox regression analysis identified 9 candidate genes. (B) MR analysis and PH assumption testing confirmed consistency in 6 candidate genes. (C) LASSO regression selected 5 prognostic genes (NR3C1, BMP2, RACGAP1, TLR3, FASN) for model construction. (D) Risk score distribution (median=1.11) in TCGA-PRAD cohort (n=278) divided into high- (HRG) and low-risk (LRG) groups. (E) Survival status distribution between risk groups. (F) Significant survival difference by Kaplan-Meier analysis. (G) ROC validation showing predictive accuracy (AUC>0.7). (H) Differential expression patterns of signature genes across risk groups, red indicates HRG, blue indicates LRG.
3.4 Validation of risk models
In order to assess the accuracy of risk prediction, the ROC curve was used to calculate the AUC value based on the TCGA-PRAD training set and the GSE116918 dataset to evaluate the model’s effectiveness. In the TCGA-PRAD-validation, on the basis of the median of the risk scores (1.056659), PCa patients were separated into a HRG (59 PCa patients) and a LRG (60 PCa patients). Patients’ risk score distribution (Figure 3A) and survival status (Figure 3B) were also displayed. The survival curves indicated that HRG had a lower survival rate than LRG (P=0.033) (Figure 3C). As the AUC values were all above 0.7, it demonstrated that the risk model had excellent performance (Figure 3D). BMP2, RACGAP1, and FASN exhibited high expression in the HRG, with TLR3 and NR3C1 showing high expression in the LRG (Figure 3E). Similarly, in all samples with relapse information of TCGA-PRAD, on the basis of the median of the risk scores (1.085715), PCa patients were categorized into a HRG (198 PCa patients) and a LRG (199 PCa patients). The risk score distribution (Figure 3F) among patients and their survival status (Figure 3G) were likewise shown. According to the survival curves, the survival rate in the HRG was markedly lower than that in the LRG (P < 0.0001) (Figure 3H). Given that the AUC values all exceeded 0.7, it was evident that the risk model performed well (Figure 3I). High expression of BMP2, RACGAP1, and FASN was observed in the HRG, while high expression of TLR3 and NR3C1 was found in the LRG (Figure 3J). Based on the results of 3.4, the risk model was accurate. In the GSE116918 dataset, there was a significant survival difference between the high and low-risk groups (p=0.014), but the ROC curves for recurrence at 1, 2, and 3 years showed that the AUC values were all less than 0.6. This suggests that although the risk score was significantly associated with survival differences, its accuracy in predicting recurrence remains limited. Further validation and optimization may be needed, possibly by integrating other clinical indicators or more comprehensive models (Supplementary Figure 5).

Figure 3. Validation of PCa recurrence risk model in TCGA-PRAD cohorts. (A) Risk score distribution stratified by median cutoff (1.057) in validation cohort (HRG: 59 patients, red; LRG: 60 patients, blue). (B) Survival status mapping across ordered risk scores. (C) Significant survival divergence between groups by Kaplan-Meier analysis (P=0.033). (D) ROC validation demonstrating robust 1-/3-/5-year recurrence prediction (AUC>0.7). (E) Heatmap contrasting elevated expression of BMP2/RACGAP1/FASN in HRG (red) versus TLR3/NR3C1 in LRG (blue). (F-H) Relapse cohort stratification (median=1.086) with risk distribution (198 HRG vs 199 LRG) and survival status. (I) Superior prognostic accuracy in relapse-specific ROC analysis (AUC>0.7). (J) Conserved expression dichotomy of signature genes across relapse subgroups.
3.5 The predictive ability of prognostic genes for PCa patients
The nomogram indicated that prognostic genes could predict the return rate of PCa patients quite well, and the recurrence rate increased as the total score of the nomogram rose (Figure 4A). The slope of all calibration curves was close to 1, indicating the accuracy of the nomogram model (Figure 4B). The AUC values of ROC curves all exceeded 0.7, which further verified the accuracy of the nomogram model (Figure 4C). In conclusion, prognostic genes could be utilized to predict the recurrence rate of PCa patients.

Figure 4. Development and validation of a gene-based nomogram for PCa recurrence prediction. (A) Prognostic nomogram integrating FASN, TLR3, NR3C1, BMP2, and RACGAP1 expression scores to estimate biochemical recurrence (BCR) probabilities at 1-, 3-, and 5-year intervals. Cumulative risk scores correlate with escalating recurrence rates (Pr[BCR]). (B) Calibration curves demonstrate high concordance between predicted and observed outcomes (slope ≈1) for 1- and 3-year BCR predictions, with narrow 95% confidence intervals. (C) ROC curves validate robust discriminative capacity, with AUC values exceeding 0.7 across all time points.
3.6 Tumor microenvironment and immune cells in PCa
In order to assess tumor purity, immune infiltration, and stromal infiltration in the malignant tumor tissue within the patient’s tumor microenvironment, immune infiltration analysis was conducted. Within the microenvironment specific to PCa tumors, immune scores, ESTIMATE scores, and the stromal scores in the LRG were all strikingly higher than those in the HRG (P < 0.05) (Figure 5A). Subsequently, the infiltration levels of 28 immune cells in the HRG and LRG were analyzed (Figure 5B). There were 18 immune cells that showed notable distinctions between the HRG and the LRG, and they were defined as differentially expressed immune cells (P < 0.05) (Figure 5C, Supplementary Table 8). Except for the activated CD4 T cells, the extents of infiltration of the remaining cells in the LRG were all remarkably higher than those in the HRG (P < 0.05). TLR3 had the largest notable positive link with natural killer cells (cor=0.66, P < 0.0001) and FASN had the largest notable negative correlation with natural killer (NK) T cells (cor=-0.38, P < 0.0001) (Figure 5D, Supplementary Table 9). In summary, the occurrence and development of PCa might be related to changes in the tumor microenvironment or some immune cells.

Figure 5. Tumor microenvironment and immune landscape stratification in PCa risk groups. (A) Violin plots demonstrating elevated stromal, immune, and ESTIMATE scores in the LRG compared to the HRG. (B) Heatmap of 28 immune cell infiltration profiles (blue: low; red: high) revealing distinct immune microenvironments between risk groups. (C) Box plots identifying 18 differentially infiltrated immune cells (P < 0.05), with LRG showing higher infiltration in all cell types except activated CD4 T cells. (D) Correlation matrix highlighting TLR3-NK cell synergy (cor=0.66, P < 0.0001) and FASN-NKT cell antagonism (cor=−0.38, P < 0.0001).
3.7 Enrichment pathways and drug sensitivity in HRG and LRG
In order to determine the biological pathways involved in the development of PCa between the high and low-risk groups, GSEA enrichment analysis was conducted. By GSEA analysis, the HRG and LRG were strikingly enriched 41 pathways such as hcm, porphyrin and chlorophyll metabolism, pentose and glucuronate interconversions, and cell cycle (Figure 6A, Supplementary Table 10). These pathways suggested that the risk of PCa might be related to some biological role. In order to obtain molecular targeted drugs corresponding to the genes in the high and low-risk groups, drug sensitivity analysis was conducted. Among 198 drugs, the HRG and LRG suggested notable disparities in sensitivity to 86 drugs (P < 0.05) (Figure 6B, Supplementary Table 11). The IC50 value of AZD8186 in the HRG was strikingly higher than that in the LRG, indicating that the LRG was more sensitive to this drug. On the contrary, the IC50 value of ML323 in the HRG was remarkably lower than that in the LRG, indicating that the HRG was more responsive to this drug.

Figure 6. Pathway enrichment and therapeutic vulnerability profiling in PCa risk groups. (A) GSEA enrichment plot showing HRG-specific activation of 41 pathways (FDR <0.25), including hypertrophic cardiomyopathy (HCM), porphyrin metabolism, and cell cycle regulation. Colored traces represent pathway-specific enrichment scores, with dashed lines indicating significance thresholds. (B) Box plots comparing IC50 values of 86 clinically actionable drugs (Wilcoxon test, P <0.05). Yellow indicates HRG; blue indicates LRG. LRG exhibits heightened AZD8186 sensitivity (HRG median IC50: 6.3 μM vs. LRG: 3.1 μM, P=0.002), while HRG shows preferential ML323 response (HRG: 0.8 μM vs. LRG: 2.4 μM, P=0.0001).
3.8 The functions and related regulatory factors of prognostic genes
In order to identify other genes related to the function of prognostic genes, the GeneMANIA database was used to predict genes associated with the function of prognostic genes and the functions they are involved in. By GeneMANIA analysis, 20 genes related to the functions of prognostic genes were acquired such as BMPR1A and these genes were related to 7 functions such as response to BMP (Figure 7A). Furthermore, in order to explore the upstream regulatory factors and their interactions for the prognostic genes, a TF-mRNA-miRNA regulatory network was constructed. Two TFs and 32 miRNAs were predicted for FASN; 5 TFs and 32 miRNAs were predicted for NR3C1; 2 TFs and 32 miRNAs were predicted for RACGAP1; 4 TFs and 10 miRNAs were predicted for BMP2; 6 TFs and 10 miRNAs were predicted for TLR3 (Figure 7B). Then, the network was constructed to demonstrate the complex relationships between prognostic genes and regulatory molecules.

Figure 7. Functional interactome and regulatory networks of PCa prognostic genes. (A) GeneMANIA-derived functional interaction network of 20 co-expressed genes (e.g., BMPR1A) enriched in 7 pathways including BMP response (red edges) and cell adhesion (blue edges). Node size reflects interaction degree, with physical interactions (gray edges) and genetic linkages (gold edges) highlighting modular biology. (B) Multi-layer regulatory network mapping transcription factors (TFs, pink nodes) and miRNAs (green nodes) targeting core prognostic genes. Circular layout emphasizes combinatorial regulation, with edge thickness proportional to prediction confidence.
3.9 Key cells in PCa
In order to ensure the accuracy, reliability, and interpretability of the single-cell data, quality control was performed on all scRNA-seq data. There were 24,391 genes and 36,424 cells before quality control (Figure 8A), and 24,391 genes and 34,571 cells after quality control (Figure 8B). Then, the 2,000 HVGs and the top 10 most varied gene names were displayed (Figure 8C). After dimensionality reduction, the first 30 PCs were used for clustering analysis (Figure 8D). Finally, the 14 cell clusters were acquired (Figure 8E). After annotating cell clusters with markers genes, there were 7 cells acquired, which included mast cells (TPSB2, MS4A2, TPSAB1), stroma cells (COL1A2, TAGLN, ACTA2), endothelial cells (PECAM1, VWF, ACKR1), T cells (GZMA, CD3E, CD3D), epithelial cells (KRT18, EPCAM, KRT19), myeloid cells (CD14, FCGR3A, CD163), and B cells (MS4A1, CD79A, DERL3) (Figures 8F-H). Among all the prognostic genes, FASN was prominently expressed in epithelial cells, so epithelial cells were considered as key cells (Figure 8I). NR3C1 was more abundantly distributed in T cells, and FASN was more abundantly distributed in epithelial cells (Figure 8J).

Figure 8. Single-cell transcriptomic profiling identifies epithelial cells as key mediators in PCa progression. (A, B) Quality control metrics showing 36,424 cells (24,391 genes) pre-filtering (A) and 34,571 cells post-filtering (B) with mitochondrial/ribosomal thresholds. (C) Top 10 highly variable genes (HVGs) among 2,000 identified, ranked by dispersion. (D) Principal component analysis using first 30 principal components (PCs) for dimensionality reduction. (E) t-SNE visualization of 14 unsupervised cell clusters. (F-H) Cell type annotation using lineage-specific markers: mast (TPSB2), stroma (COL1A2), endothelial (PECAM1), T cells (CD3E), epithelial (KRT18), myeloid (CD14), B cells (MS4A1). (I) Violin plots revealing FASN overexpression (red gradient) in epithelial clusters versus other cell types. (J) Spatial expression mapping showing NR3C1 enrichment (blue) in T cells and FASN dominance (red) in epithelial compartments.
3.10 Communication networks, differentiation of key cells, and expression of prognostic genes in key cells
In order to identify the cell types of the samples in the dataset GSE141445 and describe the cellular states of the clustering results, annotations were made for the seven different cell clustering results. Among the annotated cells, the engagements between endothelial cells and epithelial cells were the most frequent (Figure 9A, Supplementary Figure 6A). Epithelial cells and myeloid cells had the strongest interactions with other cells (Figure 9B). Epithelial cells had the strongest interaction with endothelial cells, T cells, myeloid cells, and B cells (Supplementary Figure 6B). Interaction between epithelial cells and B cells was carried out by MIF−(CD74+CXCR4) (Figure 9C). Through reduction and clustering, epithelial cells were eventually categorized into 12 clusters (Figures 9D, E). Further investigation into the developmental trajectory of the key cell, epithelial cell, revealed that epithelial cells varied from dark blue to light blue during differentiation and were categorized into 9 stages and 12 clusters (Figure 9F). During the differentiation of epithelial cells, only the expression level of FASN continuously increased, and it remained relatively active throughout the entire cell development stage (Figure 9G). In conclusion, the development of PCa might be related to epithelial cells and FASN. Furthermore, the expression patterns of prognostic genes during the pseudotime process of myeloid cells differentiation were demonstrated. As shown in the figure, the FASN and TLR3 genes exhibited higher expression levels at the early stages of myeloid cell differentiation, while RACGAP1 and NR3C1 genes had higher expression at the late stages, and BMP2 gene showed higher expression at the mid-stage of differentiation. These results suggested that these five prognostic genes had elevated expression at specific stages of myeloid cell differentiation, which might have led to their overall expression being less prominent (Supplementary Figure 7).

Figure 9. Epithelial cell communication dynamics and FASN-driven progression in PCa. (A) Cell-cell interaction network showing dominant endothelial-epithelial crosstalk (edge thickness proportional to interaction frequency). (B) Force-directed graph quantifying epithelial cells as central interactors with myeloid, T, and B cells. (C) Ligand-receptor validation of epithelial-B cell communication via MIF-(CD74+CXCR4) axis (red: ligand expression; blue: receptor activity). (D-E) UMAP visualization of 12 epithelial subclusters (D) and pseudo-temporal ordering (E) from quiescent (dark blue) to advanced malignant states (light blue). (F) Pseudotime trajectory map resolving 9 differentiation stages across 12 subclusters. (G) Lineage-specific expression kinetics identifying FASN as the sole prognostic gene with monotonic upregulation, contrasting static patterns of NR3C1/RACGAP1/TLR3/BMP2.
3.11 Validation of prognostic gene expression
In order to determine the expression patterns of prognostic genes in PCa samples and control samples, expression analysis was conducted using the dataset and RT-qPCR experiments. The expression levels of BMP2, NR3C1 and TLR3 were all remarkably lower in PCa samples than in the control samples (P < 0.0001), while the expression levels of RACGAP1 and FASN were remarkably higher in PCa samples than in the control samples (P < 0.0001) (Figure 10A). In RT-qPCR, the expression levels of NR3C1 (P < 0.0001), BMP2 (P < 0.01) and TLR3 (P < 0.01) in PCa were significantly lower than those in control group, while the expression levels of RACGAP1 (P < 0.05) and FASN (P < 0.01) in PCa were significantly higher than those in control group (Figure 10B). The expression level and trend of prognostic genes in vitro samples were consistent with that of bioinformatics analysis, indicating the reliability of bioinformatics analysis results. Previous studies have highlighted the significance of BMP2 and FASN in PCa progression. Horvath et al. demonstrated that reduced BMP2 expression is associated with PCa progression, linking its loss to more aggressive phenotypes (52). Similarly, Tae et al. reported that decreased BMP2 expression correlates with a higher incidence of biochemical recurrence (BCR) and elevated Gleason scores (GS) (53). However, the precise mechanisms underlying BMP2’ s role as a tumor outcome determinant remain elusive. In contrast, the biological role of FASN as a key regulator of lipid metabolism is well-established. By driving lipid synthesis, FASN provides energy to fuel tumor proliferation and progression. Despite this, its specific regulatory effects within the PCa microenvironment are not fully understood. Chianese et al. observed FASN overexpression in PCa and proposed that FASN inhibition disrupts the metabolic axis, leading to lipid accumulation and subsequent lipotoxicity (54, 55). This metabolic dysregulation impairs replication mechanisms and arrests cells in the G0/G1 phase, thereby inhibiting proliferation (56, 57). These findings underscore the multifaceted roles of BMP2 and FASN in PCa biology and warrant further investigation into their underlying mechanisms and therapeutic potential.

Figure 10. Experimental validation of prognostic gene expression patterns in PCa. (A) Box plots of RNA-seq expression levels in control (teal) and PCa (orange) groups. Whiskers extend to 1.5×IQR; dots indicate outliers. (B) RT-qPCR confirmation using 15 paired PCa/control specimens. Bar graphs quantify concordant expression trends.
3.12 Overexpression of BMP2 or silencing of FASN suppresses malignant behaviors of PCa cells
Western blot analysis revealed that compared with the Vector group, the protein level of BMP2 was significantly upregulated in the OE-BMP2 group, while the FASN protein level showed no obvious change; conversely, in the sh-FASN group, the FASN protein level was notably downregulated compared with the sh-NC group, with no significant difference in BMP2 level (Figures 11A-D). Immunofluorescence staining demonstrated that the proportion of Ki67-positive cells was significantly higher in the OE-BMP2 group than in the Vector group (Figures 11E, G), and lower in the sh-FASN group than in the sh-NC group (Figures 11F, H), as quantified by statistical analysis. In the cell migration assay, the number of migrated cells in the OE-BMP2 group was significantly greater than that in the Vector group at both 24 h and 48 h time points, whereas the sh-FASN group showed a significantly reduced number of migrated cells compared with the sh-NC group at the same time points, as shown by the quantitative results (Figures 11I-L).

Figure 11. Overexpression of BMP2 or silencing of FASN suppresses malignant behaviors of PC3 cells. (A–D) Western blot results confirm that BMP2 overexpression effectively increases BMP2 protein levels in PC3 cells, whereas FASN knockdown effectively reduces FASN protein levels. (E–H) Ki67 immunostaining shows that BMP2 overexpression or FASN knockdown both suppress PC3 cell proliferative activity. (I–L) Scratch-wound assay indicate that BMP2 overexpression or FASN knockdown significantly reduces the migratory capacity of PC3 cells. Data are presented as mean ± SEM, n=3. **p < 0.01.
4 Discussion
PCa remains a leading contributor to global cancer-related morbidity and mortality among males, with rising incidence rates documented in recent epidemiological studies (58). The clinical management of PCa faces significant challenges due to the lack of discernible clinical manifestations in early-stage disease, resulting in delayed diagnosis for the majority of patients until intermediate/advanced phases or metastatic progression, which is the critical factor driving elevated mortality (59). These clinical realities underscore the urgent need to develop novel therapeutic strategies and personalized treatment paradigms to improve patient outcomes. In this study, we integrated transcriptomic profiling and scRNA-seq data to systematically identify MCDRGs with causal associations in PCa pathogenesis. A predictive risk model was developed and rigorously validated to assess the clinical utility of these biomarkers. Through comprehensive bioinformatic interrogation, we elucidated the mechanistic contributions of candidate genes to PCa progression, complemented by single-cell resolution analysis of their expression dynamics across tumor-associated cellular subpopulations. The computational findings were further substantiated through in vitro functional validation, confirming the biological relevance of identified molecular networks.
NR3C1 (Nuclear Receptor Subfamily 3 Group C Member 1), located at chromosome 5q31.3-q32, encodes the glucocorticoid receptor (GR), the sole mediator of glucocorticoid signaling through ligand-activated transcriptional regulation within the nuclear receptor superfamily (60). This receptor-ligand complex translocates to the nucleus, binding specific DNA response elements to orchestrate diverse physiological processes including glucose/lipid metabolism, inflammatory responses, and cellular differentiation (60, 61). Emerging evidence positions NR3C1 as a pivotal oncogenic regulator across malignancies: it drives progression in triple-negative breast cancer, ovarian carcinoma, urothelial cancer, and clear cell renal carcinoma (61–65), while mediating platinum and targeted therapy resistance in lung and ovarian cancers (66, 67). In PCa, a dynamic AR-NR3C1 axis governs therapeutic resistance. Androgen receptor (AR) suppresses NR3C1 expression in treatment-naïve states, whereas androgen deprivation therapy induces compensatory GR upregulation—a critical mechanism enabling treatment evasion through AR-GR crosstalk (68, 69). Mechanistically, Qian et al. delineated the role of NR3C1 in PCa lineage plasticity, demonstrating that the ONECUT2 (OC2) transcription factor activates NR3C1 and neuroendocrine splicing factor SRRM4 to drive adenocarcinoma progression and therapy-resistant stem-like/neuroendocrine variants (70).
BMP2 (Bone Morphogenetic Protein 2), a key TGF-β superfamily member located at chromosome 20p12, encodes a multifunctional regulator of cellular processes including proliferation, differentiation, migration, and apoptosis, though it remains best characterized for its osteoinductive role in skeletal development (71). During embryogenesis, BMP2 drives osteogenic differentiation of mesenchymal stem cells by stimulating extracellular matrix production (collagen, osteocalcin) and subsequent bone mineralization (72, 73). Beyond developmental biology, BMP2 exhibits remarkable therapeutic potential in fracture repair through site-specific upregulation that recruits progenitor cells and accelerates osteogenesis (74). The pleiotropic nature of BMP2 signaling manifests through context-dependent tumor modulation. While suppressing metastasis in breast cancer (75, 76) and inhibiting proliferation/biochemical recurrence in PCa (53, 77), its aberrant activation paradoxically enhances hepatocellular carcinoma progression via proliferative and invasive mechanisms (78, 79). Moreover, recent evidence extends the functional repertoire of BMP2 to fibroblast biology, demonstrating anti-inflammatory properties that mitigate atrial fibrosis (80).
RACGAP1 (Rac GTPase-Activating Protein 1), a critical regulator of cellular dynamics, encodes a member of the GTPase-activating protein (GAP) family that modulates Rac GTPase activity to control cytoskeletal reorganization and mitotic fidelity (81). During cell division, RACGAP1 ensures genomic stability through its indispensable role in spindle assembly and chromosome segregation (82, 83). Emerging oncogenic roles of RACGAP1 span multiple malignancies. Its dysregulated expression correlates with aggressive phenotypes in lung adenocarcinoma and bladder cancer, where it drives tumor progression via enhanced proliferation, invasion, and metastatic dissemination. In hepatocellular carcinoma, RACGAP1 emerges as both an independent prognostic biomarker and a modulator of tumor immune microenvironment (84–86). Notably, RACGAP1 intersects with therapeutic resistance pathways in PCa. Mechanistic studies reveal its capacity to activate downstream effectors of the PI3K/AKT axis, a compensatory signaling network implicated in ADT resistance and neuroendocrine differentiation (87). Clinical validation through qPCR analysis confirms RACGAP1 overexpression in castration-resistant PCa (CRPC), underscoring its functional relevance in treatment-refractory disease (88).
TLR3 (Toll-like Receptor 3), located at chromosome 4q35.1, encodes a pattern recognition receptor predominantly expressed on immune cells including dendritic cells, NK cells, and macrophages (89). This receptor initiates antiviral immunity through specific recognition of viral double-stranded RNA (dsRNA), triggering signal transduction cascades that activate NF-κB and induce interferon/cytokine production (90). Emerging evidence reveals context-dependent prognostic implications of TLR3 dysregulation across malignancies, low TLR3 expression correlates with favorable outcomes in gastric, prostate, and breast cancers, yet paradoxically associates with poor prognosis in clear cell renal carcinoma and hepatocellular carcinoma (91–93). Muresan et al. analyzed hormone-naïve and hormone-resistant PCa specimens, demonstrating TLR3 upregulation in therapy-refractory tumors alongside its functional role in promoting migratory and invasive capacities (94). These findings align with prior clinical observations documenting TLR3’s association with aggressive PCa behavior (95–97). Intriguingly, comparative studies reveal tissue-specific expression patterns, While González-Reyes et al. reported elevated TLR3 levels in PCa versus benign prostate tissue (97), our data paradoxically demonstrate significant TLR3 downregulation in tumor versus adjacent paracancerous tissues, which correlates with changes in the tumor microenvironment and immune evasion mechanisms, suggests that TLR3 suppression may play a critical role in the malignant transformation of prostate epithelium. As a critical bridge between innate and adaptive immunity, TLR3 plays a pivotal role in anti-tumor immunity (98). Extensive research has demonstrated that TLR3 can directly activate tumor-specific NK cells or mediate the release of interferon to enhance cytotoxic lymphocyte (CTL) infiltration and response, establish type 1 T helper cells (Th1) immunity, and upregulate genes involved in the recruitment and functionality of immune cells within the tumor microenvironment (99, 100). Notably, TLR3 executes its functions through diverse immune pathways. However, its expression at the tissue or organismal level may exhibit ectopic distribution, such as variations in subcellular localization (e.g., cytoplasmic versus membrane-bound) or differential expression between cell types (101, 102). Importantly, the functional impact of TLR3 cannot be comprehensively inferred from its overall expression levels alone, as shifts in cellular composition or differential expression across specific cell types may obscure its true immunological significance under various conditions. These complexities underscore the need for a nuanced exploration of the role of TLR3 in tumor immunity.
FASN (Fatty Acid Synthase), the rate-limiting enzyme catalyzing the final step of de novo fatty acid synthesis, is ubiquitously expressed with elevated activity in lipid-metabolizing organs including liver, adipose tissue, and mammary glands (103, 104). This multienzyme complex mediates the NADPH-dependent condensation of acetyl-CoA and malonyl-CoA to generate palmitate, the primary substrate for membrane phospholipid synthesis, energy storage, and bioactive lipid precursors (104). Under physiological conditions, FASN activity is tightly regulated, with endogenous synthesis suppressed under nutrient-replete conditions via insulin-mediated regulation to prioritize dietary fatty acid utilization. In oncogenic contexts, FASN undergoes pathological upregulation to fuel tumorigenic demands. Elevated FASN expression drives de novo lipogenesis, fulfilling the biosynthetic requirements of rapidly proliferating cancer cells for membrane remodeling and signaling lipid generation (105, 106). Beyond its role in lipid synthesis, FASN promotes lipid accumulation within tumor cells, which may disrupt antigen presentation or alter surface molecule profiles, thereby impairing the immune system’ s ability to recognize these cells (107, 108). Moreover, FASN-mediated metabolic reprogramming reshapes the tumor microenvironment through lipid-driven immunosuppression, where increased FASN activity inversely correlates with antitumor immune cell infiltration across multiple malignancies (109, 110). These findings position FASN as a compelling therapeutic target, with pharmacological inhibition strategies showing promise for disrupting cancer-specific lipogenic dependencies.
Our investigation revealed 18 differentially abundant immune cell populations between LRG and HRG in the PCa microenvironment, with NK cells demonstrating the most significant correlations with prognostic gene signatures. Unlike T cells requiring antigen-specific MHC recognition, NK cells execute innate immunosurveillance through non-antigen-directed cytotoxicity against malignant cells. Notably, NK cells constitute 2-9% of tumor-infiltrating lymphocytes in prostate carcinoma, underscoring their microenvironmental relevance (111). Gannon et al. reported reduced biochemical recurrence rates in treatment-naïve patients exhibiting elevated intraprostatic NK cell infiltration (112). Complementary studies by Pasero et al. linked high surface NKG2D expression on tumor-associated NK cells with attenuated disease progression (113). Mechanistically, Lundholm et al. identified tumor-derived exosomal NKG2D ligands as immunosuppressive agents that downregulate NK cell activation via receptor internalization—a plausible immune evasion mechanism (114). Multidimensional profiling by Zorko et al. further delineated NK cell-mediated clinical benefits, enhanced NK infiltration inversely correlated with driver mutations (e.g., AR-V7 variants) in primary tumors while positively associating with immunosuppressive checkpoints, suggesting dual roles in tumor editing and microenvironment modulation (111). These findings collectively nominate NK cell potentiation strategies—including endogenous activity enhancement and adoptive cell therapies—as promising therapeutic frontiers for CRPC. Additionally, in tumor immunotherapy, mitochondrial function in NK cells directly impacts their activity, survival, and antitumor capacity (115, 116). Thus, optimizing mitochondrial health in NK cells may emerge as a potential strategy to enhance therapeutic efficacy.
Our analysis identified clinically divergent drug sensitivities between HRG and LRG, including AZD8186 and JAK1. AZD8186, a potent selective PI3Kβ inhibitor, targets the oncogenic PI3K signaling axis implicated in tumor cell proliferation, metabolic adaptation, and angiogenesis (117–119). The first-in-human trial (NCT01884285) established its manageable safety profile and preliminary efficacy (120). Preclinically, Ruiz et al. demonstrated synergistic antitumor activity of AZD8186 combined with selumetinib in docetaxel-resistant murine models without additive toxicity, highlighting its therapeutic potential for taxane-refractory prostate cancer (121). Mechanistic insights into chemoresistance emerged from single-cell profiling of docetaxel-resistant tumors by Cheng et al. revealing IL-11 overexpression that activates the JAK1/STAT4 axis. This cascade facilitates STAT4-CBP complex formation, driving c-MYC transcription—a well-characterized oncogene promoting tumorigenesis and therapy resistance (122–124). IL-11 further orchestrates a chemoresistant niche via autocrine tumor cell signaling and paracrine stromal interactions involving extracellular matrix remodeling (125, 126). Therapeutic opportunitieslie in disrupting this axis through JAK1 inhibition, IL-11 neutralization, or STAT4-CBP interface targeting. Future work should employ advanced humanized models recapitulating tumor-immune-stroma crosstalk to validate these strategies and delineate microenvironmental influences on drug response.
GeneMANIA analysis identified 20 functionally interconnected genes associated with prognostic signatures, including BMPR1A—a pivotal receptor mediating BMP signaling through ligand binding and pathway activation (127). Yang et al. mechanistically demonstrated that GALNT12 enhances BMPR1A O-glycosylation to suppress metastatic PCa cell proliferation, migration, and invasion, nominating GALNT12 as a therapeutic target (128). Notably, several network components like GREM2 and FSTL1 remain underexplored in PCa contexts, warranting further investigation of their therapeutic potential. The reconstructed regulatory network revealed transcription factors with established oncogenic roles. Yin Yang 1 (YY1), a C2H2 zinc finger transcriptional regulator implicated in tumor-associated immune suppression, was shown to drive IL-6 production in M2-polarized macrophages while inhibiting anti-tumor T-cell activity—a mechanism suggesting YY1-targeted immunomodulation strategies (14, 129, 130). Furthermore, CUX1, a homeodomain transcription factor governing development and cell cycle progression (131), emerged as a network hub differentially regulated during ADT. Sharma et al. identified CUX1 as a key transcriptional regulator through whole-transcriptome profiling of pre-/post-ADT specimens (132), while Dorris et al. reported paradoxical effects: CUX1 knockdown enhanced migration in androgen-sensitive cells but increased invasion in castration-resistant lineages (133). These context-dependent phenotypes underscore the need for mechanistic elucidation of CUX1’s dual roles in PCa progression.
Epithelial cells, integral components of innate immune surveillance and tissue barrier defense (134, 135), emerge as pivotal mediators of PCa progression through multifaceted mechanisms. In this study, the results form scRNA-seq analyze revealed three interconnected oncogenic roles: (1) Prognostic gene FASN exhibits sustained overexpression in malignant epithelia, directly driving tumorigenic behaviors; (2) Epithelial cells orchestrate a pro-tumorigenic niche via crosstalk with endothelial and immune compartments; (3) Persistent FASN activation during epithelial differentiation suggests dynamic involvement in malignant transformation. These findings align with growing interest in PCa immunotherapy, though its unique immunosuppressive microenvironment poses translational challenges (136–138). Zhu et al. pioneered an epithelial cell marker gene prognostic signature (ECMGPS) derived from scRNA-seq data (GSE176031), demonstrating robust predictive accuracy for immunotherapy responses despite lacking clinical cohort validation (139). Mechanistic insights from Jiang et al. established its capacity to induce epithelial-mesenchymal transition (EMT) via transcriptional repression of E-cadherin and N-cadherin activation, facilitating peritoneal metastasis in ovarian cancer (140). To date, no studies have mechanistically delineated FASN-epithelial cell interactions in PCa. Our study addresses this critical knowledge gap by providing the first functional evidence linking FASN activity to epithelial cell-mediated oncogenic progression in PCa pathogenesis.
Experimental validation using RT-qPCR confirmed the consistency between bioinformatic predictions and experimental findings, with dysregulated expression of BMP2, NR3C1, TLR3, RACGAP1, and FASN closely associated with PCa progression. Specifically, BMP2 was found to suppress tumor cell proliferation and migration, as evidenced by downregulation of Ki67 expression and reduced migratory activity in PC3 cells. Conversely, FASN promoted both proliferation and migration, underscoring its role as a driver of tumor aggressiveness. These findings highlight the dual regulatory roles of BMP2 and FASN in PCa pathobiology, suggesting their potential as both biomarkers and therapeutic targets. The observed concordance between bioinformatic analysis and experimental validation strengthens the reliability of the identified prognostic genes. Furthermore, the mechanistic insights into BMP2 and FASN’s opposing effects on proliferation and migration provide a foundation for the development of targeted therapeutic strategies. Future studies should prioritize mechanistic investigation of these biomarkers, including development of targeted agents, analysis of stage-specific expression dynamics, and optimization of combinatorial therapeutic strategies. Such efforts will advance precision oncology frameworks to address unmet clinical needs in PCa management.
5 Conclusions
This study innovatively integrated MR with transcriptomic and scRNA-seq analyses to identify five prognostic genes associated with PCa progression, subsequently constructing a risk stratification model and evaluating its clinical utility. Mechanistic exploration revealed functional roles of these genes in tumorigenic pathways, complemented by scRNA-seq analysis of their expression patterns across PCa-associated cellular subpopulations. While these findings advance our understanding of PCa biology, several limitations warrant consideration. First, the prognostic model relied solely on a single dataset without external validation cohorts, potentially limiting its generalizability. Second, experimental validation focused primarily on preliminary functional assays, necessitating further exploration of the molecular mechanisms underlying BMP2 and FASN-mediated effects. Despite these constraints, our study establishes a robust theoretical foundation for the development of personalized therapeutic strategies in PCa management. By identifying and characterizing key prognostic genes, we provide a comprehensive framework for advancing precision oncology and addressing unmet clinical needs in PCa prognosis and treatment.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Ethics statement
The studies involving humans were approved by The Clinical Ethics Committee of Tangdu Hospital of Air Force Medical University. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
JC: Validation, Data curation, Visualization, Conceptualization, Investigation, Writing – review & editing, Software, Formal analysis, Writing – original draft. JQ: Writing – review & editing, Data curation, Conceptualization. WZ: Visualization, Conceptualization, Writing – review & editing, Project administration. ZN: Writing – original draft, Data curation, Methodology. XG: Writing – original draft, Validation. GX: Supervision, Writing – original draft. LK: Investigation, Writing – original draft. ZZ: Writing – review & editing, Resources, Validation, Software, Writing – original draft, Supervision, Project administration.
Funding
The author(s) declare that no financial support was received for the research, and/or publication of this article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1619194/full#supplementary-material
References
1. Gandaglia G, Leni R, Bray F, Fleshner N, Freedland SJ, Kibel A, et al. Epidemiology and prevention of prostate cancer. Eur Urol Oncol. (2021) 4:877–92. doi: 10.1016/j.euo.2021.09.006
2. Bergengren O, Pekala KR, Matsoukas K, Fainberg J, Mungovan SF, Bratt O, et al. 2022 Update on prostate cancer epidemiology and risk factors-A systematic review. Eur Urol. (2023) 84:191–206. doi: 10.1016/j.eururo.2023.04.021
3. Kulasegaran T and Oliveira N. Metastatic castration-resistant prostate cancer: advances in treatment and symptom management. Curr Treat options Oncol. (2024) 25:914–31. doi: 10.1007/s11864-024-01215-2
4. Cha HR, Lee JH, and Ponnazhagan S. Revisiting immunotherapy: A focus on prostate cancer. Cancer Res. (2020) 80:1615–23. doi: 10.1158/0008-5472.Can-19-2948
5. Wen XY, Wang RY, Yu B, Yang Y, Yang J, and Zhang HC. Integrating single-cell and bulk RNA sequencing to predict prognosis and immunotherapy response in prostate cancer. Sci Rep. (2023) 13:15597. doi: 10.1038/s41598-023-42858-9
6. Wang W, Li T, Xie Z, Zhao J, Zhang Y, Ruan Y, et al. Integrating single-cell and bulk RNA sequencing data unveils antigen presentation and process-related CAFS and establishes a predictive signature in prostate cancer. J Trans Med. (2024) 22:57. doi: 10.1186/s12967-023-04807-y
7. Zhang Y, Fan A, Li Y, Liu Z, Yu L, Guo J, et al. Single-cell RNA sequencing reveals that HSD17B2 in cancer-associated fibroblasts promotes the development and progression of castration-resistant prostate cancer. Cancer Lett. (2023) 566:216244. doi: 10.1016/j.canlet.2023.216244
8. Ciccarese C, Massari F, Iacovelli R, Fiorentino M, Montironi R, Di Nunno V, et al. Prostate cancer heterogeneity: Discovering novel molecular targets for therapy. Cancer Treat Rev. (2017) 54:68–73. doi: 10.1016/j.ctrv.2017.02.001
9. McCabe N, Hanna C, Walker SM, Gonda D, Li J, Wikstrom K, et al. Mechanistic rationale to target PTEN-deficient tumor cells with inhibitors of the DNA damage response kinase ATM. Cancer Res. (2015) 75:2159–65. doi: 10.1158/0008-5472.Can-14-3502
10. De Cicco P, Ercolano G, and Ianaro A. The new era of cancer immunotherapy: targeting myeloid-derived suppressor cells to overcome immune evasion. Front Immunol. (2020) 11:1680. doi: 10.3389/fimmu.2020.01680
11. Goswami S, Anandhan S, Raychaudhuri D, and Sharma P. Myeloid cell-targeted therapies for solid tumours. Nat Rev Immunol. (2023) 23:106–20. doi: 10.1038/s41577-022-00737-w
12. Hicks KC, Tyurina YY, Kagan VE, and Gabrilovich DI. Myeloid cell-derived oxidized lipids and regulation of the tumor microenvironment. Cancer Res. (2022) 82:187–94. doi: 10.1158/0008-5472.Can-21-3054
13. Singer D, Ressel V, Stope MB, and Bekeschus S. Heat shock protein 27 affects myeloid cell activation and interaction with prostate cancer cells. Biomedicines. (2022) 10:2192. doi: 10.3390/biomedicines10092192
14. Chen S, Lu K, Hou Y, You Z, Shu C, Wei X, et al. YY1 complex in M2 macrophage promotes prostate cancer progression by upregulating IL-6. J immunother Cancer. (2023) 11:e006020. doi: 10.1136/jitc-2022-006020
15. Cui H, Zhang W, Zhang L, Qu Y, Xu Z, Tan Z, et al. Risk factors for prostate cancer: An umbrella review of prospective observational studies and mendelian randomization analyses. PloS Med. (2024) 21:e1004362. doi: 10.1371/journal.pmed.1004362
16. Miao M, Song Y, Jin M, Du Y, Xin P, Jiang Y, et al. Single-cell RNA combined with bulk RNA analysis to explore oxidative stress and energy metabolism factors and found a new prostate cancer oncogene MXRA8. Aging. (2024) 16:4469–502. doi: 10.18632/aging.205599
17. Xin S, Liu X, Li Z, Sun X, Wang R, Zhang Z, et al. ScRNA-seq revealed an immunosuppression state and tumor microenvironment heterogeneity related to lymph node metastasis in prostate cancer. Exp Hematol Oncol. (2023) 12:49. doi: 10.1186/s40164-023-00407-0
18. Zhang N, Liu X, Qin J, Sun Y, Xiong H, Lin B, et al. LIGHT/TNFSF14 promotes CAR-T cell trafficking and cytotoxicity through reversing immunosuppressive tumor microenvironment. Mol therapy: J Am Soc Gene Ther. (2023) 31:2575–90. doi: 10.1016/j.ymthe.2023.06.015
19. Ye Z, Deng X, Zhang J, Shao R, Song C, Zhao J, et al. Causal relationship between immune cells and prostate cancer: a Mendelian randomization study. Front Cell Dev Biol. (2024) 12:1381920. doi: 10.3389/fcell.2024.1381920
20. Wu D, Liu Y, Liu J, Ma L, and Tong X. Myeloid cell differentiation-related gene signature for predicting clinical outcome, immune microenvironment, and treatment response in lung adenocarcinoma. Sci Rep. (2024) 14:17460. doi: 10.1038/s41598-024-68111-5
21. Love MI, Huber W, and Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
22. Gustavsson EK, Zhang D, Reynolds RH, Garcia-Ruiz S, and Ryten M. ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinf (Oxford England). (2022) 38:3844–6. doi: 10.1093/bioinformatics/btac409
23. Gu Z, Eils R, and Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinf (Oxford England). (2016) 32:2847–9. doi: 10.1093/bioinformatics/btw313
24. Mao W, Ding J, Li Y, Huang R, and Wang B. Inhibition of cell survival and invasion by Tanshinone IIA via FTH1: A key therapeutic target and biomarker in head and neck squamous cell carcinoma. Exp Ther Med. (2022) 24:521. doi: 10.3892/etm.2022.11449
25. Wang L, Wang D, Yang L, Zeng X, Zhang Q, Liu G, et al. Cuproptosis related genes associated with Jab1 shapes tumor microenvironment and pharmacological profile in nasopharyngeal carcinoma. Front Immunol. (2022) 13:989286. doi: 10.3389/fimmu.2022.989286
26. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Cambridge (Mass.)). (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
27. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
28. Hemani G, Zheng J, Elsworth B, Wade KH, Haberland V, Baird D, et al. The MR-Base platform supports systematic causal inference across the human phenome. eLife. (2018) 7:e34408. doi: 10.7554/eLife.34408
29. Burgess S, Scott RA, Timpson NJ, Davey Smith G, and Thompson SG. Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors. Eur J Epidemiol. (2015) 30:543–52. doi: 10.1007/s10654-015-0011-z
30. Hartwig FP, Davey Smith G, and Bowden J. Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption. Int J Epidemiol. (2017) 46:1985–98. doi: 10.1093/ije/dyx102
31. Bowden J, Davey Smith G, and Burgess S. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. Int J Epidemiol. (2015) 44:512–25. doi: 10.1093/ije/dyv080
32. Bowden J, Davey Smith G, Haycock PC, and Burgess S. Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol. (2016) 40:304–14. doi: 10.1002/gepi.21965
33. Lei J, Qu T, Cha L, Tian L, Qiu F, Guo W, et al. Clinicopathological characteristics of pheochromocytoma/paraganglioma and screening of prognostic markers. J Surg Oncol. (2023) 128:510–8. doi: 10.1002/jso.27358
34. Ye Y, Dai Q, and Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov. (2021) 7:71. doi: 10.1038/s41420-021-00451-x
35. Kang SJ, Cho YR, Park GM, Ahn JM, Han SB, Lee JY, et al. Predictors for functionally significant in-stent restenosis: an integrated analysis using coronary angiography, IVUS, and myocardial perfusion imaging. JACC Cardiovasc Imaging. (2013) 6:1183–90. doi: 10.1016/j.jcmg.2013.09.006
36. Jiang H, Li Y, and Wang T. Comparison of Billroth I, Billroth II, and Roux-en-Y reconstructions following distal gastrectomy: A systematic review and network meta-analysis. Cirugia espanola. (2021) 99:412–20. doi: 10.1016/j.cireng.2020.09.018
37. Engebretsen S and Bohlin J. Statistical predictions with glmnet. Clin Epigenet. (2019) 11:123. doi: 10.1186/s13148-019-0730-1
38. Liu TT, Li R, Huo C, Li JP, Yao J, Ji XL, et al. Identification of CDK2-related immune forecast model and ceRNA in lung adenocarcinoma, a pan-cancer analysis. Front Cell Dev Biol. (2021) 9:682002. doi: 10.3389/fcell.2021.682002
39. Heagerty PJ, Lumley T, and Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. (2000) 56:337–44. doi: 10.1111/j.0006-341x.2000.00337.x
40. Sui Z, Wu X, Du L, Wang H, Yuan L, Zhang JV, et al. Characterization of the immune cell infiltration landscape in esophageal squamous cell carcinoma. Front Oncol. (2022) 12:879326. doi: 10.3389/fonc.2022.879326
41. Xu J, Yang T, Wu F, Chen T, Wang A, and Hou S. A nomogram for predicting prognosis of patients with cervical cerclage. Heliyon. (2023) 9:e21147. doi: 10.1016/j.heliyon.2023.e21147
42. Hänzelmann S, Castelo R, and Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7
43. Liu S, Meng Y, Zhang Y, Qiu L, Wan X, Yang X, et al. Integrative analysis of senescence-related genes identifies robust prognostic clusters with distinct features in hepatocellular carcinoma. J advanced Res. (2025) 69:107–23. doi: 10.1016/j.jare.2024.04.007
44. Robles-Jimenez LE, Aranda-Aguirre E, Castelan-Ortega OA, Shettino-Bermudez BS, Ortiz-Salinas R, Miranda M, et al. Worldwide traceability of antibiotic residues from livestock in wastewater and soil: A systematic review. Animals: an Open Access J MDPI. (2021) 12:60. doi: 10.3390/ani12010060
45. Geeleher P, Cox N, and Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One. (2014) 9:e107468. doi: 10.1371/journal.pone.0107468
46. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. (2021) 184:3573–3587.e3529. doi: 10.1016/j.cell.2021.04.048
47. Yang J, Gong L, Liu Q, Zhao H, Wang Z, Li X, et al. Single-cell RNA-seq reveals developmental deficiencies in both the placentation and the decidualization in women with late-onset preeclampsia. Front Immunol. (2023) 14:1142273. doi: 10.3389/fimmu.2023.1142273
48. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9
49. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402
50. Fagherazzi G, Aguayo GA, Zhang L, Hanaire H, Picard S, Sablone L, et al. Heterogeneity of glycaemic phenotypes in type 1 diabetes. Diabetologia. (2024) 67:1567–81. doi: 10.1007/s00125-024-06179-4
51. Jani Kargar Moghaddam S, Mohammadi Roushandeh A, Hamidi M, Nemati S, Jahanian-Najafabadi A, and Habibi Roudkenar M. Lipocalin-2 upregulation in nasopharyngeal carcinoma: A novel potential diagnostic biomarker. Iranian J Med Sci. (2023) 48:268–76. doi: 10.30476/ijms.2022.93041.2452
52. Horvath LG, Henshall SM, Kench JG, Turner JJ, Golovsky D, Brenner PC, et al. Loss of BMP2, Smad8, and Smad4 expression in prostate cancer progression. Prostate. (2004) 59:234–42. doi: 10.1002/pros.10361
53. Tae BS, Cho S, Kim HC, Kim CH, Kang SH, Lee JG, et al. Decreased expression of bone morphogenetic protein-2 is correlated with biochemical recurrence in prostate cancer: Immunohistochemical analysis. Sci Rep. (2018) 8:10748. doi: 10.1038/s41598-018-28566-9
54. Chianese U, Papulino C, Ali A, Ciardiello F, Cappabianca S, Altucci L, et al. FASN multi-omic characterization reveals metabolic heterogeneity in pancreatic and prostate adenocarcinoma. J Trans Med. (2023) 21:32. doi: 10.1186/s12967-023-03874-5
55. Khan F, Pandey P, Ahmad V, and Upadhyay TK. Moringa oleifera methanolic leaves extract induces apoptosis and G0/G1 cell cycle arrest via downregulation of Hedgehog Signaling Pathway in human prostate PC-3 cancer cells. J Food Biochem. (2020) 44:e13338. doi: 10.1111/jfbc.13338
56. Yoon H, Shaw JL, Haigis MC, and Greka A. Lipid metabolism in sickness and in health: Emerging regulators of lipotoxicity. Mol Cell. (2021) 81:3708–30. doi: 10.1016/j.molcel.2021.08.027
57. Dierge E and Feron O. Dealing with saturated and unsaturated fatty acid metabolism for anticancer therapy. Curr Opin Clin Nutr Metab Care. (2019) 22:427–33. doi: 10.1097/mco.0000000000000601
58. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: Cancer J Clin. (2021) 71:209–49. doi: 10.3322/caac.21660
59. Pishgar F, Ebrahimi H, Saeedi Moghaddam S, Fitzmaurice C, and Amini E. Global, regional and national burden of prostate cancer, 1990 to 2015: results from the global burden of disease study 2015. J Urol. (2018) 199:1224–32. doi: 10.1016/j.juro.2017.10.044
60. Hollenberg SM, Weinberger C, Ong ES, Cerelli G, Oro A, Lebo R, et al. Primary structure and expression of a functional human glucocorticoid receptor cDNA. Nature. (1985) 318:635–41. doi: 10.1038/318635a0
61. Yan M, Wang J, Wang H, Zhou J, Qi H, Naji Y, et al. Knockdown of NR3C1 inhibits the proliferation and migration of clear cell renal cell carcinoma through activating endoplasmic reticulum stress-mitophagy. J Trans Med. (2023) 21:701. doi: 10.1186/s12967-023-04560-2
62. Karvonen H, Arjama M, Kaleva L, Niininen W, Barker H, Koivisto-Korander R, et al. Glucocorticoids induce differentiation and chemoresistance in ovarian cancer by promoting ROR1-mediated stemness. Cell Death Dis. (2020) 11:790. doi: 10.1038/s41419-020-03009-4
63. Obradović MMS, Hamelin B, Manevski N, Couto JP, Sethi A, Coissieux MM, et al. Glucocorticoids promote breast cancer metastasis. Nature. (2019) 567:540–4. doi: 10.1038/s41586-019-1019-4
64. Chien HP, Ueng SH, Chen SC, Chang YS, Lin YC, Lo YF, et al. Expression of ROR1 has prognostic significance in triple negative breast cancer. Virchows Archiv: an Int J Pathol. (2016) 468:589–95. doi: 10.1007/s00428-016-1911-3
65. Zheng Y, Izumi K, Li Y, Ishiguro H, and Miyamoto H. Contrary regulation of bladder cancer cell proliferation and invasion by dexamethasone-mediated glucocorticoid receptor signals. Mol Cancer Ther. (2012) 11:2621–32. doi: 10.1158/1535-7163.Mct-12-0621
66. Prekovic S, Schuurman K, Mayayo-Peralta I, Manjón AG, Buijs M, Yavuz S, et al. Glucocorticoid receptor triggers a reversible drug-tolerant dormancy state with acquired therapeutic vulnerabilities in lung cancer. Nat Commun. (2021) 12:4360. doi: 10.1038/s41467-021-24537-3
67. Pan C, Kang J, Hwang JS, Li J, Boese AC, Wang X, et al. Cisplatin-mediated activation of glucocorticoid receptor induces platinum resistance via MAST1. Nat Commun. (2021) 12:4960. doi: 10.1038/s41467-021-24845-8
68. Arora VK, Schenkein E, Murali R, Subudhi SK, Wongvipat J, Balbas MD, et al. Glucocorticoid receptor confers resistance to antiandrogens by bypassing androgen receptor blockade. Cell. (2013) 155:1309–22. doi: 10.1016/j.cell.2013.11.012
69. Shah N, Wang P, Wongvipat J, Karthaus WR, Abida W, Armenia J, et al. Regulation of the glucocorticoid receptor via a BET-dependent enhancer drives antiandrogen resistance in prostate cancer. eLife. (2017) 6:e27861. doi: 10.7554/eLife.27861
70. Qian C, Yang Q, Rotinen M, Huang R, Kim H, Gallent B, et al. ONECUT2 acts as a lineage plasticity driver in adenocarcinoma as well as neuroendocrine variants of prostate cancer. Nucleic Acids Res. (2024) 52:7740–60. doi: 10.1093/nar/gkae547
71. Jiramongkolchai P, Owens P, and Hong CC. Emerging roles of the bone morphogenetic protein pathway in cancer: potential therapeutic target for kinase inhibition. Biochem Soc Trans. (2016) 44:1117–34. doi: 10.1042/bst20160069
72. Ray A, Singh PN, Sohaskey ML, Harland RM, and Bandyopadhyay A. Precise spatial restriction of BMP signaling is essential for articular cartilage differentiation. Dev (Cambridge England). (2015) 142:1169–79. doi: 10.1242/dev.110940
73. Kobayashi T, Lyons KM, McMahon AP, and Kronenberg HM. BMP signaling stimulates cellular differentiation at multiple steps during cartilage development. Proc Natl Acad Sci United States America. (2005) 102:18023–7. doi: 10.1073/pnas.0503617102
74. Yoon BS and Lyons KM. Multiple functions of BMPs in chondrogenesis. J Cell Biochem. (2004) 93:93–103. doi: 10.1002/jcb.20211
75. Gao H, Chakraborty G, Lee-Lim AP, Mo Q, Decker M, Vonica A, et al. The BMP inhibitor Coco reactivates breast cancer cells at lung metastatic sites. Cell. (2012) 150:764–79. doi: 10.1016/j.cell.2012.06.035
76. Tarragona M, Pavlovic M, Arnal-Estapé A, Urosevic J, Morales M, Guiu M, et al. Identification of NOG as a specific breast cancer bone metastasis-supporting gene. J Biol Chem. (2012) 287:21346–55. doi: 10.1074/jbc.M112.355834
77. Feng H, Deng Z, Peng W, Wei X, Liu J, and Wang T. Circular RNA EPHA3 suppresses progression and metastasis in prostate cancer through the miR-513a-3p/BMP2 axis. J Trans Med. (2023) 21:288. doi: 10.1186/s12967-023-04132-4
78. Li P, Shang Y, Yuan L, Tong J, and Chen Q. Targeting BMP2 for therapeutic strategies against hepatocellular carcinoma. Trans Oncol. (2024) 46:101970. doi: 10.1016/j.tranon.2024.101970
79. Wu G, Huang F, Chen Y, Zhuang Y, Huang Y, and Xie Y. High levels of BMP2 promote liver cancer growth via the activation of myeloid-derived suppressor cells. Front Oncol. (2020) 10:194. doi: 10.3389/fonc.2020.00194
80. Yuan Y, Zhang H, Xia E, Zhao X, Gao Q, Mu H, et al. BMP2 diminishes angiotensin II-induced atrial fibrillation by inhibiting NLRP3 inflammasome signaling in atrial fibroblasts. Biomolecules. (2024) 14:1053. doi: 10.3390/biom14091053
81. Jantsch-Plunger V, Gönczy P, Romano A, Schnabel H, Hamill D, Schnabel R, et al. CYK-4: A Rho family gtpase activating protein (GAP) required for central spindle formation and cytokinesis. J Cell Biol. (2000) 149:1391–404. doi: 10.1083/jcb.149.7.1391
82. Lee JS, Kamijo K, Ohara N, Kitamura T, and Miki T. MgcRacGAP regulates cortical activity through RhoA during cytokinesis. Exp Cell Res. (2004) 293:275–82. doi: 10.1016/j.yexcr.2003.10.015
83. Zhao WM and Fang G. MgcRacGAP controls the assembly of the contractile ring and the initiation of cytokinesis. Proc Natl Acad Sci United States America. (2005) 102:13158–63. doi: 10.1073/pnas.0504145102
84. Chen H, Wu J, Lu L, Hu Z, Li X, Huang L, et al. Identification of hub genes associated with immune infiltration and predict prognosis in hepatocellular carcinoma via bioinformatics approaches. Front Genet. (2020) 11:575762. doi: 10.3389/fgene.2020.575762
85. Eid RA, Soltan MA, Eldeen MA, Shati AA, Dawood SA, Eissa M, et al. Assessment of RACGAP1 as a prognostic and immunological biomarker in multiple human tumors: A multiomics analysis. Int J Mol Sci. (2022) 23:14102. doi: 10.3390/ijms232214102
86. Son JA, Ahn HR, You D, Baek GO, Yoon MG, Yoon JH, et al. Novel gene signatures as prognostic biomarkers for predicting the recurrence of hepatocellular carcinoma. Cancers. (2022) 14:865. doi: 10.3390/cancers14040865
87. Hazar-Rethinam M, de Long LM, Gannon OM, Boros S, Vargas AC, Dzienis M, et al. RacGAP1 is a novel downstream effector of E2F7-dependent resistance to doxorubicin and is prognostic for overall survival in squamous cell carcinoma. Mol Cancer Ther. (2015) 14:1939–50. doi: 10.1158/1535-7163.Mct-15-0076
88. Pan J, Zhang J, Lin J, Cai Y, and Zhao Z. Constructing lactylation-related genes prognostic model to effectively predict the disease-free survival and treatment responsiveness in prostate cancer based on machine learning. Front Genet. (2024) 15:1343140. doi: 10.3389/fgene.2024.1343140
89. Iwasaki A and Medzhitov R. Toll-like receptor control of the adaptive immune responses. Nat Immunol. (2004) 5:987–95. doi: 10.1038/ni1112
90. Chew V and Abastado JP. Immunomodulation of the tumor microenvironment by Toll-like receptor-3 (TLR3) ligands. Oncoimmunology. (2013) 2:e23493. doi: 10.4161/onci.23493
91. Liu R, Liu Y, Huang W, Chen P, and Cheng Y. An anoikis-related signature predicts prognosis and immunotherapy response in gastrointestinal cancers. Front Immunol. (2025) 16:1477913. doi: 10.3389/fimmu.2025.1477913
92. Liu Y, Zhou Q, Ye F, Yang C, and Jiang H. Gut microbiota-derived short-chain fatty acids promote prostate cancer progression via inducing cancer cell autophagy and M2 macrophage polarization. Neoplasia (New York N.Y.). (2023) 43:100928. doi: 10.1016/j.neo.2023.100928
93. Zou X, Guo Y, and Mo Z. TLR3 serves as a novel diagnostic and prognostic biomarker and is closely correlated with immune microenvironment in three types of cancer. Front Genet. (2022) 13:905988. doi: 10.3389/fgene.2022.905988
94. Muresan XM, Slabáková E, Procházková J, Drápela S, Fedr R, Pícková M, et al. Toll-like receptor 3 overexpression induces invasion of prostate cancer cells, whereas its activation triggers apoptosis. Am J Pathol. (2022) 192:1321–35. doi: 10.1016/j.ajpath.2022.05.009
95. Escudero-Lourdes C, Alvarado-Morales I, and Tokar EJ. Stem cells as target for prostate cancer therapy: opportunities and challenges. Stem Cell Rev Rep. (2022) 18:2833–51. doi: 10.1007/s12015-022-10437-6
96. Ilvesaro JM, Merrell MA, Swain TM, Davidson J, Zayzafoon M, Harris KW, et al. Toll like receptor-9 agonists stimulate prostate cancer invasion in vitro. Prostate. (2007) 67:774–81. doi: 10.1002/pros.20562
97. González-Reyes S, Fernández JM, González LO, Aguirre A, Suárez A, González JM, et al. Study of TLR3, TLR4, and TLR9 in prostate carcinomas and their association with biochemical recurrence. Cancer immunol immunother: CII. (2011) 60:217–26. doi: 10.1007/s00262-010-0931-0
98. Smith M, García‑Martínez E, Pitter MR, Fucikova J, Spisek R, Zitvogel L, et al. Trial Watch: Toll-like receptor agonists in cancer immunotherapy. Oncoimmunology. (2018) 7:e1526250. doi: 10.1080/2162402x.2018.1526250
99. Le Naour J, Galluzzi L, Zitvogel L, Kroemer G, and Vacchelli E. Trial watch: TLR3 agonists in cancer therapy. Oncoimmunology. (2020) 9:1771143. doi: 10.1080/2162402x.2020.1771143
100. Seya T, Takeda Y, and Matsumoto M. A Toll-like receptor 3 (TLR3) agonist ARNAX for therapeutic immunotherapy. Advanced Drug delivery Rev. (2019) 147:37–43. doi: 10.1016/j.addr.2019.07.008
101. Stanifer ML, Mukenhirn M, Muenchau S, Pervolaraki K, Kanaya T, Albrecht D, et al. Asymmetric distribution of TLR3 leads to a polarized immune response in human intestinal epithelial cells. Nat Microbiol. (2020) 5:181–91. doi: 10.1038/s41564-019-0594-3
102. Mielcarska MB, Bossowska-Nowicka M, and Toka FN. Cell surface expression of endosomal toll-like receptors-A necessity or a superfluous duplication? Front Immunol. (2020) 11:620972. doi: 10.3389/fimmu.2020.620972
103. Fhu CW and Ali A. Fatty acid synthase: an emerging target in cancer. Molecules (Basel Switzerland). (2020) 25:3935. doi: 10.3390/molecules25173935
104. Gonzalez-Bohorquez D, Gallego López IM, Jaeger BN, Pfammatter S, Bowers M, Semenkovich CF, et al. FASN-dependent de novo lipogenesis is required for brain development. Proc Natl Acad Sci United States America. (2022) 119:e2112040119. doi: 10.1073/pnas.2112040119
105. Lim SA, Wei J, Nguyen TM, Shi H, Su W, Palacios G, et al. Lipid signalling enforces functional specialization of T(reg) cells in tumours. Nature. (2021) 591:306–11. doi: 10.1038/s41586-021-03235-6
106. Ferraro GB, Ali A, Luengo A, Kodack DP, Deik A, Abbott KL, et al. FATTY ACID SYNTHESIS IS REQUIRED FOR BREAST CANCER BRAIN METASTASIS. Nat Cancer. (2021) 2:414–28. doi: 10.1038/s43018-021-00183-y
107. Luo Y, He X, Du Q, Xu L, Xu J, Wang J, et al. Metal-based smart nanosystems in cancer immunotherapy. Explor (Beijing China). (2024) 4:20230134. doi: 10.1002/exp.20230134
108. Huang J, Tsang WY, Fang XN, Zhang Y, Luo J, Gong LQ, et al. FASN inhibition decreases MHC-I degradation and synergizes with PD-L1 checkpoint blockade in hepatocellular carcinoma. Cancer Res. (2024) 84:855–71. doi: 10.1158/0008-5472.Can-23-0966
109. Hinshaw DC and Shevde LA. The tumor microenvironment innately modulates cancer progression. Cancer Res. (2019) 79:4557–66. doi: 10.1158/0008-5472.Can-18-3962
110. Zhang M, Yu L, Sun Y, Hao L, Bai J, Yuan X, et al. Comprehensive analysis of FASN in tumor immune infiltration and prognostic value for immunotherapy and promoter DNA methylation. Int J Mol Sci. (2022) 23:15603. doi: 10.3390/ijms232415603
111. Zorko NA, Makovec A, Elliott A, Kellen S, Lozada JR, Arafa AT, et al. Natural killer cell infiltration in prostate cancers predict improved patient outcomes. Prostate Cancer prostatic Dis. (2025) 28:129–37. doi: 10.1038/s41391-024-00797-0
112. Gannon PO, Poisson AO, Delvoye N, Lapointe R, Mes-Masson AM, and Saad F. Characterization of the intra-prostatic immune cell infiltration in androgen-deprived prostate cancer patients. J Immunol Methods. (2009) 348:9–17. doi: 10.1016/j.jim.2009.06.004
113. Pasero C, Gravis G, Guerin M, Granjeaud S, Thomassin-Piana J, Rocchi P, et al. Inherent and tumor-driven immune tolerance in the prostate microenvironment impairs natural killer cell antitumor activity. Cancer Res. (2016) 76:2153–65. doi: 10.1158/0008-5472.Can-15-1965
114. Lundholm M, Schröder M, Nagaeva O, Baranov V, Widmark A, Mincheva-Nilsson L, et al. Prostate tumor-derived exosomes down-regulate NKG2D expression on natural killer cells and CD8+ T cells: mechanism of immune evasion. PloS One. (2014) 9:e108925. doi: 10.1371/journal.pone.0108925
115. Terrén I, Sandá V, Amarilla-Irusta A, Lopez-Pardo A, Sevilla A, Astarloa-Pando G, et al. IL-12/15/18-induced cell death and mitochondrial dynamics of human NK cells. Front Immunol. (2023) 14:1211839. doi: 10.3389/fimmu.2023.1211839
116. Qian K, Gao S, Jiang Z, Ding Q, and Cheng Z. Recent advances in mitochondria-targeting theranostic agents. Explor (Beijing China). (2024) 4:20230063. doi: 10.1002/exp.20230063
117. Jia S, Liu Z, Zhang S, Liu P, Zhang L, Lee SH, et al. Essential roles of PI(3)K-p110beta in cell growth, metabolism and tumorigenesis. Nature. (2008) 454:776–9. doi: 10.1038/nature07091
118. Liu P, Cheng H, Roberts TM, and Zhao JJ. Targeting the phosphoinositide 3-kinase pathway in cancer. Nat Rev Drug Discov. (2009) 8:627–44. doi: 10.1038/nrd2926
119. Vanhaesebroeck B, Stephens L, and Hawkins P. PI3K signalling: the path to discovery and understanding. Nat Rev Mol Cell Biol. (2012) 13:195–203. doi: 10.1038/nrm3290
120. Choudhury AD, Higano CS, de Bono JS, Cook N, Rathkopf DE, Wisinski KB, et al. A phase I study investigating AZD8186, a potent and selective inhibitor of PI3Kβ/δ, in patients with advanced solid tumors. Clin Cancer Res. (2022) 28:2257–69. doi: 10.1158/1078-0432.Ccr-21-3087
121. Ruiz de Porras V, Bernat-Peguera A, Alcon C, Laguia F, Fernández-Saorin M, Jiménez N, et al. Dual inhibition of MEK and PI3Kβ/δ-a potential therapeutic strategy in PTEN-wild-type docetaxel-resistant metastatic prostate cancer. Front Pharmacol. (2024) 15:1331648. doi: 10.3389/fphar.2024.1331648
122. Palcau AC, Canu V, Donzelli S, Strano S, Pulito C, and Blandino G. CircPVT1: a pivotal circular node intersecting Long Non-Coding-PVT1 and c-MYC oncogenic signals. Mol Cancer. (2022) 21:33. doi: 10.1186/s12943-022-01514-y
123. Cao LL and Kagan JC. Targeting innate immune pathways for cancer immunotherapy. Immunity. (2023) 56:2206–17. doi: 10.1016/j.immuni.2023.07.018
124. Cheng B, Li L, Luo T, Wang Q, Luo Y, Bai S, et al. Single-cell deconvolution algorithms analysis unveils autocrine IL11-mediated resistance to docetaxel in prostate cancer via activation of the JAK1/STAT4 pathway. J Exp Clin Cancer research: CR. (2024) 43:67. doi: 10.1186/s13046-024-02962-8
125. Li M, Chi X, Wang Y, Setrerrahmane S, Xie W, Xu H, et al. Trends in insulin resistance: insights into mechanisms and therapeutic strategy. Signal transduction targeted Ther. (2022) 7:216. doi: 10.1038/s41392-022-01073-0
126. Fan J, To KKW, Chen ZS, and Fu L. ABC transporters affects tumor immune microenvironment to regulate cancer immunotherapy and multidrug resistance. Drug resistance updates: Rev commentaries antimicrobial Anticancer chemother. (2023) 66:100905. doi: 10.1016/j.drup.2022.100905
127. Jiang B, Zhao X, Chen W, Diao W, Ding M, Qin H, et al. Lysosomal protein transmembrane 5 promotes lung-specific metastasis by regulating BMPR1A lysosomal degradation. Nat Commun. (2022) 13:4141. doi: 10.1038/s41467-022-31783-6
128. Yang Y, Ding M, Yin H, Chen W, Shen H, Diao W, et al. GALNT12 suppresses the bone-specific prostate cancer metastasis by activating BMP pathway via the O-glycosylation of BMPR1A. Int J Biol Sci. (2024) 20:1297–313. doi: 10.7150/ijbs.91925
129. Gordon S, Akopyan G, Garban H, and Bonavida B. Transcription factor YY1: structure, function, and therapeutic implications in cancer biology. Oncogene. (2006) 25:1125–42. doi: 10.1038/sj.onc.1209080
130. Wu S, Wang H, Li Y, Xie Y, Huang C, Zhao H, et al. Transcription factor YY1 promotes cell proliferation by directly activating the pentose phosphate pathway. Cancer Res. (2018) 78:4549–62. doi: 10.1158/0008-5472.Can-17-4047
131. Michl P, Ramjaun AR, Pardo OE, Warne PH, Wagner M, Poulsom R, et al. CUTL1 is a target of TGF(beta) signaling that enhances cancer cell motility and invasiveness. Cancer Cell. (2005) 7:521–32. doi: 10.1016/j.ccr.2005.05.018
132. Sharma NV, Pellegrini KL, Ouellet V, Giuste FO, Ramalingam S, Watanabe K, et al. Identification of the transcription factor relationships associated with androgen deprivation therapy response and metastatic progression in prostate cancer. Cancers. (2018) 10:379. doi: 10.3390/cancers10100379
133. Dorris ER, O'Neill A, Treacy A, Klocker H, Teltsh O, Kay E, et al. The transcription factor CUX1 negatively regulates invasion in castrate resistant prostate cancer. Oncotarget. (2020) 11:846–57. doi: 10.18632/oncotarget.27494
134. Whitsett JA and Alenghat T. Respiratory epithelial cells orchestrate pulmonary innate immunity. Nat Immunol. (2015) 16:27–35. doi: 10.1038/ni.3045
135. Hammad H and Lambrecht BN. Barrier epithelial cells and the control of type 2 immunity. Immunity. (2015) 43:29–40. doi: 10.1016/j.immuni.2015.07.007
136. Bansal D, Reimers MA, Knoche EM, and Pachynski RK. Immunotherapy and immunotherapy combinations in metastatic castration-resistant prostate cancer. Cancers. (2021) 13:334. doi: 10.3390/cancers13020334
137. Runcie KD and Dallos MC. Prostate cancer immunotherapy-finally in from the cold? Curr Oncol Rep. (2021) 23:88. doi: 10.1007/s11912-021-01084-0
138. Mitsogiannis I, Tzelves L, Dellis A, Issa H, Papatsoris A, and Moussa M. Prostate cancer immunotherapy. Expert Opin Biol Ther. (2022) 22:577–90. doi: 10.1080/14712598.2022.2027904
139. Zhu W, Zeng H, Huang J, Wu J, Wang Y, Wang Z, et al. Integrated machine learning identifies epithelial cell marker genes for improving outcomes and immunotherapy in prostate cancer. J Trans Med. (2023) 21:782. doi: 10.1186/s12967-023-04633-2
Keywords: prostate cancer, myeloid cell differentiation, prognostic genes, Mendelian randomization, single-cell sequencing analysis, experimental verification
Citation: Chen J, Qiu J, Zhang W, Nie Z, Gao X, Xu G, Kang L and Zhang Z (2025) Mendelian randomization combined with single-cell sequencing analysis revealed prognostic genes related to myeloid cell differentiation in prostate cancer and experimental verification. Front. Immunol. 16:1619194. doi: 10.3389/fimmu.2025.1619194
Received: 27 April 2025; Accepted: 21 August 2025;
Published: 23 September 2025.
Edited by:
Anand Rotte, Arcellx Inc, United StatesReviewed by:
Fuyi Xu, University of Tennessee Health Science Center (UTHSC), United StatesBo-yu Yang, Capital Medical University, China
Copyright © 2025 Chen, Qiu, Zhang, Nie, Gao, Xu, Kang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhiming Zhang, d293Y2h1bmd6em1AMTYzLmNvbQ==