Systematic characterization of neurotransmitter receptor dysregulation identifies a neural-related prognostic signature associated with biochemical recurrence in prostate cancer

Introduction

Prostate cancer (PCa) is the most commonly diagnosed malignancy and the second leading cause of cancer-related mortality among men worldwide (1). In 2026, an estimated 333,830 new PCa cases and 36,320 PCa-related deaths are expected in the United States (1). Currently, radical prostatectomy, radiation therapy, active surveillance, and androgen deprivation therapy represent the first line treatments for localized PCa (2). For patients with advanced disease, therapeutic strategies also include second-generation anti-androgens, radiation therapy, and chemotherapy (2). Despite these advances, many patients eventually develop treatment resistance and experience biochemical recurrence (BCR). BCR occurs in approximately 20–40% of treated patients, and significantly affects both patient survival and quality of life (3). Therefore, early identification of the patients at a high risk of BCR serves an essential role in improving patient prognosis.

The nervous system has recently been recognized as playing an important role in tumor initiation and progression (4), including PCa. Increasing evidence indicates that nerve cells and fibers can infiltrate the tumor microenvironment, and influence diverse biological processes in cancer (5). Central to the intricate relationship between the nervous system and tumors are the interactions between neurotransmitters secreted by neurons and their receptors expressed on cancer cells and other cell types via synapses, which subsequently activate multiple intracellular signaling pathways (6). Neurotransmitter receptors (NTRs), also known as neuroreceptors, are membrane receptor proteins activated by neurotransmitters. Previous studies have shown that gamma-aminobutyric acid (GABA), the major inhibitory neurotransmitter, and its receptor are upregulated in PCa (7). Furthermore, the GABA agonist baclofen has been reported to enhance the invasive potential of PCa cells by promoting EGFR transactivation (8). In addition, Palamiuc and Emerling demonstrated that the glutamate-glutamate receptor axis can activate phosphoinositide 3-kinase through phosphorylation of p110β, thereby facilitating PCa progression (9). Collectively, these findings suggest that NTR genes may provide a promising strategy for identifying novel biomarkers for PCa prognosis and treatment.

Currently, risk stratification for BCR in PCa is primarily based on clinicopathological parameters (10), including prostate-specific antigen (PSA), Gleason score (GS), and pathological T stage. However, classification of BCR risk groups using GS and PSA doubling time, as recommended by the European Association of Urology guidelines, only partially accounts for the variability in clinical outcomes and lacks sufficient accuracy in predicting high-risk BCR (11). Therefore, strategies relying solely on clinicopathological parameters require further refinement.

In this study, we systematically investigated NTR dysregulation in PCa and identified NTR genes associated with BCR. We identified 13 consistently dysregulated NTR (cdNTR) genes, and developed a PCa signature (PCaSig) that was significantly associated with unfavorable BCR outcomes, aggressive clinical and molecular features, and elevated tumor mutation burden (TMB). Single-cell transcriptomic analyses of 167,527 cells derived from 51 prostate specimens obtained from 37 patients further confirmed robust dysregulation of these cdNTR genes, particularly within tumor epithelial compartments. The prognostic performance of PCaSig was subsequently validated in two independent cohorts. By combining PCaSig stratification with TMB, we identified optimal BCR outcomes in patients with both low PCaSig risk and low TMB. Notably, PCaSig was able to stratify PCa patients with distinct BCR outcomes regardless of TMB status and was associated with immune-related characteristics of the tumor microenvironment. Furthermore, we constructed a well-calibrated nomogram integrating PCaSig with GS and pathological T stage, providing a more individualized and accurate approach for BCR risk prediction in PCa patients. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0461/rc).

Methods

Collection of NTR genes

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. NTR genes were obtained from the HGNC gene family database (https://www.genenames.org/), yielding a total of 132 genes. After excluding two pseudogenes, 130 genes were retained for subsequent analyses (Table S1).

Data preparation and processing

The HTSeq-Counts data from The Cancer Genome Atlas Prostate Adenocarcinoma (TCGA-PRAD) project were obtained from the Genomic Data Commons and processed using the functions provided in the R package GDCRNATools (12). The count data of 547 samples were normalized with the Trimmed Mean of M-values (TMM) method implemented in the R package edgeR (13), and lowly-expressed genes with logCPM <0 in more than 50% of samples were filtered out prior to downstream analysis. For the DKFZ-PRAD cohort, the normalized gene expression data were directly downloaded from cBioPortal (https://www.cbioportal.org/). The microarray dataset GSE21034 and GSE54460 were retrieved by the R package GEOquery, and preprocessed with the Robust Multichip Average (RMA) method in the R package oligo (14), including background correction, quantile normalization, and log₂ transformation. When multiple probes corresponded to the same gene ID, only the probe with the highest interquartile range (IQR) was retained to represent the gene expression value.

Acquisition of scRNAseq data for prostate specimen

Single-cell transcriptomic data from the studies of Henry et al. (15), Tuong et al. (16), Song et al. (17), and Chen et al. (18), were downloaded from the GEO database under accession numbers GSE117403, GSE176031, GSE141445, as well as from the Prostate Cell Atlas (https://www.prostatecellatlas.org/). After excluding 1,616 cells annotated as sperm cells due to seminal fluid contamination in the original study, a total of 195,822 single cell transcriptomes were retained, comprising 34 tumor samples and 17 normal samples from 37 patients. Detailed sample information is provided in Table S2.

scRNAseq data preprocessing and quality control

Quality control and downstream analyses were performed using the Seurat package in R. Genes detected in fewer than ten cells were excluded, and cells with low-complexity libraries [defined as those with transcripts detected in fewer than 300 genes or with fewer than 500 unique molecular identifiers (UMIs)] were removed from subsequent analyses. In addition, cells with mitochondrial RNA content exceeding 20% were filtered out. Following initial clustering, putative cell doublets were identified and removed from all clusters based on multiple criteria: (I) library complexity, defined as cells exhibiting outlier characteristics, specifically those with more than 6,000 expressed genes or over 50,000 UMIs; (II) cluster distribution, whereby doublets or multiplets formed distinct clusters with hybrid transcriptional features and aberrantly high gene counts; and (III) cluster marker gene expression, in which cells within a cluster simultaneously expressed canonical markers from distinct lineages (e.g., cells in a T-cell cluster expressing epithelial markers). Canonical marker gene expression patterns were carefully examined using uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding (t-SNE) embeddings, and the above filtering procedures were iteratively repeated to ensure robust exclusion of barcodes associated with cell doublets. The induction of cellular stress during the dissociation of solid tissues into single-cell suspensions is an unavoidable artifact of scRNA-seq experiments and can be mitigated either by removing affected cells prior to sequencing or by applying in silico filtering strategies (19). To identify and remove stressed cells in silico, single-sample gene set enrichment analysis (ssGSEA) was performed using stress-related signature gene sets, and highly stressed cells defined as the top 5% of cells with the highest stress scores, were excluded from downstream analyses (15). After quality control, a total of 167,527 high-quality cells were retained for downstream analyses.

Dimensionality reduction, clustering, and cell-type annotation

Filtered single cells were analyzed using the standard Seurat workflow. Briefly, the filtered gene-cell expression matrix was normalized for sequencing depth by dividing each cell’s total number of UMIs and subsequently log-transformed using the NormalizeData and ScaleData functions. Highly variable genes were identified using the FindVariableGenes function with default parameters. Dimensionality reduction was performed by principal component analysis (PCA) using the RunPCA function. To correct for batch effects and align shared cell types and states across datasets, data integration across experimental batches was carried out using the IntegrateLayers function with the parameter method set to HarmonyIntegration. Cells were clustered using FindNeighbors with the first 20 principal components (PCs) and FindClusters with a resolution of 0.8. For visualization, dimensionality was further reduced using either t-SNE or UMAP via the RunTSNE and RunUMAP functions, respectively, using the same PCs as those employed for clustering. To define major cell types, differentially expressed genes (DEGs) were identified for each cluster using the FindAllMarkers function (genes detected in at least 25% of cells and a log fold-change threshold of 0.25), applying the Wilcoxon rank sum test with Bonferroni correction [false discovery rate (FDR) <0.05]. The top 20 most significant DEGs for each cluster were carefully reviewed. In parallel, feature plots were generated for the top-ranked DEGs and a curated set of canonical cell-type markers, followed by manual inspection. Enrichment of canonical markers (e.g., KRT3/KRT17/KRT19 for epithelial cells, CD3D/CD3E/TRAC for T cells, MS4A1/CD79A/IGHM for B cells, FBLN1/LUM/DCN for fibroblasts, LYZ/CD68/C1QA for myeloid cells, CPA3/KIT/TPSAB1 for mast cells, MYH11/RGS5/TPM2 for smooth muscle cells, and IFI27/ENG/CLDN5 for endothelial cells) within specific clusters was considered strong evidence for corresponding cell-type identifies. These two complementary approaches were integrated to assign major cell types to each cell cluster. Differential gene expression analyses between conditions were calculated using the Wilcoxon rank-sum test with Bonferroni correction, implemented via the FindMarkers function with the following parameters min.pct =0.05 and logfc.threshold =0.25. A minimum of 50 cells per group was required for comparison. Genes with log fold change (FC) >0.25 and FDR <0.05 were considered differentially expressed. All dot plots were generated using the DotPlot function, expression overlays on UMAP or t-SNE embeddings were visualized using the FeaturePlot function, and violin plots were created using the VlnPlot function.

ssGSEA of single cell data

ssGSEA was performed using the escape R package, with the method parameter set of UCell (20,21). Stress-related signature gene sets, and cell type specific signature gene sets, derived from a previous single-cell profiling study of normal prostate tissue, were used to assess cellular stress levels and to validate major cell subtypes (15). Additionally, hallmark gene sets, categorized into multiple biological themes, including epithelial cell identity and differentiation, malignant transformation and invasion, proliferation and cell cycle, metabolic reprogramming, and genomic stability and therapeutic response, were obtained from the Molecular Signatures Database (MSigDB) to characterize the core features of PCa epithelial cells. The associations between the expression of cdNTR genes and the enrichment scores of the selected hallmark gene sets were calculated using Spearman correlation analysis in PCa luminal epithelial cells.

Development and validation of the PCaSig

A flowchart summarizing the analytical workflow of this study is presented in Figure 1A. Differential expression analysis was performed in the TCGA-PRAD discovery cohort, which included 52 tumor and matched normal samples. The limma package was applied (22), and genes with an absolute log₂FC >1.0 and an adjusted P value <0.05 were defined as DEGs. The use of paired-sample comparisons helped mitigate the influence of potential confounding factors. For validation, expression profiles from 29 tumor and matched normal samples in the GSE21034 dataset were analyzed. Univariate Cox proportional hazards regression analysis was conducted to analyze the association between the expression levels of the DEGs and time to BCR, with P<0.05 considered statistically significant. Significant DEGs were further refined using the least absolute shrinkage and selection operator (LASSO) method to construct the PCaSig. The risk score for each patient was calculated as a linear combination of gene expression values, weighted by the model coefficients derived from the training cohort. Patients were then stratified into low- and high-risk groups according to the median risk score. The prognostic value of PCaSig was assessed using Kaplan–Meier method and log-rank test. To further evaluate its predictive performance, PCaSig was applied to two validation cohorts (DKFZ-PRAD and GSE54460), which included 118 and 97 patients, respectively. The predictive performance of the signature was assessed via Harrell’s concordance index (C-index) and time-dependent receiver operating characteristic (ROC) curve analysis.

Figure 1 Strategy for developing PCaSig. (A) Schematic overview of the strategy used to construct and validate the PCaSig signature. (B) Heatmap of the expression profiles of 16 dysregulated NTR genes in 52 paired tumor and normal samples from the TCGA-PRAD cohort. (C) Volcano plot illustrating the 16 dysregulated NTR genes. (D) Venn diagram showing 13 consistently dysregulated NTR genes, with red indicating upregulated genes and blue indicating downregulated genes. AUC, area under the curve; BCR, biochemical recurrence; CI, confidence interval; GS, Gleason score; LASSO, least absolute shrinkage and selection operator; NTR, neurotransmitter receptor; PCaSig, prostate cancer signature; ROC, receiver operating characteristic; T, tumor; TCGA, The Cancer Genome Atlas.

Collection of previously published BCR prognostic signatures

Fifteen previously published multigene signatures for predicting BCR in PCa were collected from the literature. Each signature was associated with a published risk-score formula (Table S3) and was originally developed using gene expression profiles together with regression coefficients derived from univariate Cox regression, multivariate Cox regression, or LASSO-Cox regression analyses. For each signature, the risk score was calculated as a linear combination of gene expression values weighted by their corresponding regression coefficients according to the following formula: , where n represents the total number of genes included in the signature. Detailed information regarding the genes, coefficients, and original publications of all collected signatures is provided in Table S3.

Functional enrichment and pathway activation evaluation

Gene set enrichment analysis (GSEA) was conducted using the clusterProfiler package in R with the gseKEGG function under default parameters, ranking genes by log₂FC. Results were visualized with the gseaplot2 function from the enrichplot package. Pathway activation was assessed by gene set variation analysis (GSVA) using the ssGSEA method implemented in the GSVA package. Differences in pathway activation between the two groups were evaluated using the Wilcoxon rank-sum test, with P values were adjusted for multiple testing with the Benjamini-Hochberg (BH) method. Adjusted P<0.05 was considered statistically significant.

Analysis of the mutational landscape

For the somatic mutation data of the TCGA-PRAD cohort, we utilized the maftools R package for data collation and calculated the TMB for each sample. Hypermutator samples were defined as tumors with mutation counts greater than 1.5 times the IQR above the third quantile (3Q + 1.5 × IQR) and were excluded from further analysis. Additionally, we distinguished high and low TMB subgroups based on median cutoffs and conducted BCR-free survival analyses in conjunction with PCaSig risk stratification.

Immune infiltration analysis

The relative proportions of 22 immune cell types were estimated in the TCGA-PRAD cohort using the CIBERSORT algorithm (23). These immune cells were subsequently aggregated into broader immune categories, including CD8⁺ T cells, CD4⁺ T cells, B cells, NK cells, plasma cells, monocytes, macrophages, dendritic cells, mast cells and neutrophils, following the guidelines proposed by Thorsson et al. (24). Furthermore, Spearman correlation analysis was performed to examine the associations between the PCaSig score and immune checkpoint molecule (ICM) expression, as well as between dysregulated NTR genes and ICMs.

Statistical analysis

Associations between variables were assessed using the Chi-squared test, Fisher’s exact test, Wilcoxon rank-sum test, or Kruskal-Wallis test, as appropriate. Survival analyses were conducted using the Kaplan-Meier method, with P values determined by the log-rank test. Hazard ratios (HRs) were estimated through univariate and multivariate Cox proportional hazards models. The nomogram and corresponding calibration maps were constructed using the R package rms. Calibration curves, generated by bootstrap resampling, were used to assess the concordance between observed and predicted survival. The C-index was calculated with the Hmisc package. Time-dependent ROC analyses and corresponding area under the curve (AUC) values were obtained using the timeROC package. Decision curve analysis was performed using the dcurve R package to determine the clinical utility by quantifying the net benefits at different threshold probabilities. All statistical tests were two-tailed, performed in R, and P<0.05 was considered statistically significant.

Results

Transcriptomic dysregulation of NTR genes in PCa

To comprehensively characterize the transcriptomic dysregulation of NTR genes in PCa, we performed differential expression analysis on 52 paired tumor and adjacent normal samples from the TCGA-PRAD discovery cohort, identifying 16 dysregulated NTR genes based on the criteria of adjusted P<0.05 and |log₂FC| >1.0. Among these genes, P2RX5, ADRB1, GABRB3, CHRNA5, ADRA2A, GABRG3, CHRNA2, ADRB2, GRIN3A, and CHRM3 were upregulated, whereas ADRA1D, ADRA1A, P2RX1, P2RX2, GABRE, and GABRP were downregulated (Figure 1B,1C). To validate these results, we analyzed an independent validation cohort (GSE21034) comprising 29 paired tumor and normal samples. In this cohort, 18 NTR genes were identified as dysregulated, 13 of which overlapped with those detected in the discovery cohort (Figure 1D). Although β-adrenergic receptor signaling has been well established in PCa progression, the roles of other neural signaling pathways remain largely unexplored (25,26). For example, studies in PCa mouse models and preliminary clinical trials have shown that ADRB2, identified here as upregulated, is an integral component of a highly redundant signaling network that promotes PCa progression and therapeutic resistance (25). Thus, these 13 overlapping genes were defined as cdNTR genes and serve as a starting point for future functional investigations of neural regulation in PCa.

Cell type specific expression and single cell validation of cdNTR genes

To explore the expression patterns of cdNTR genes, we analyzed 167,527 single-cell transcriptomic profiles derived from 51 prostate specimens obtained from 37 patients. Graph-based clustering combined with canonical cell marker-based annotation identified eight major cell types (Figure 2A,2B; Figure S1A), including epithelial cells (n=117,247), fibroblasts (n=13,144), smooth muscle cells (n=11,572), T cells (n=9,242), endothelial cells (n=1,332), myeloid cells (n=4,867), mast cells (n=2,035) and B cells (n=1,223). Epithelial cells exhibited pronounced transcriptomic heterogeneity and were therefore subjected to subclustering analysis, which resolved four major epithelial subtypes, including luminal, basal, hillock, and club cells (Figure 2A,2B; Figure S1A,S1B). These subtypes were defined based on established marker genes and ssGSEA using signature gene sets developed from a previous single-cell profiling study of normal prostate tissue (15). The cell identities identified here were highly consistent with those reported in prior studies (16-18). A dot plot summarizing the expression of 13 cdNTR genes revealed that these genes were expressed in distinct cell types, generally at relatively low levels (Figures 2C,2D; Figure S1C). With the exception of ADRB2, which was broadly expressed across multiple cell types, most cdNTR genes exhibited preferential expression in one or two specific cell populations. Notably, luminal epithelial cells expressed 8 of the 13 cdNTR genes, indicating a strong enrichment of neural signaling components within this compartment.

Figure 2 cdNTR genes exhibit cell-type specific expression changes. (A) UMAP visualization of filtered single cells, colored by annotated cell types. (B) UMAP visualization showing the expression of canonical marker genes used to define the identified cell types. (C) Dot plots illustrating the proportion of each cell type expressing cdNTR genes (indicated by dot size) and their scaled expression levels (indicated by dot color). Dots enclosed by black dashed boxes indicate those analyzed in this study. (D) UMAP visualization showing the expression patterns of representative cdNTR genes in luminal epithelial cells. (E) Dot plots depicting the expression of eight cdNTR genes that are differentially expressed between tumor and normal luminal epithelial cells. Dot size represents the percentage of cells expressing the gene within each condition, while dot color indicates the average expression level. See Table S3 for details. (F) Heatmap showing the Spearman correlation coefficients between cdNTR gene expression and the enrichment scores of selected hallmark gene sets in PCa luminal epithelial cells. cdNTR, consistently dysregulated neurotransmitter receptor; FDR, false discovery rate; PCa, prostate cancer; UMAP, uniform manifold approximation and projection.

Given that the 13 cdNTR genes were initially identified at the bulk transcriptomic level based on differential expression between paired tumor and normal tissues, we further investigated their expression dynamics at single-cell resolution. We found that 9 of the 13 cdNTR genes exhibited significant expression differences between tumor and normal tissues within specific cell types, particularly in luminal epithelial cells, largely consistent with the bulk-level findings (Figure 2E; Figure S1D,S1E; Table S4). Notably, all eight cdNTR genes expressed in luminal epithelial cells were significantly upregulated in tumor tissue (Figure 2E). In addition, CHRM3 was upregulated in endothelial cells, whereas P2RX5 showed increased expression in tumor-derived T cells (Figure S1D,S1E). In contrast, the four cdNTR genes identified as downregulated in bulk tissue did not display sufficiently robust cell-type specific downregulation, likely reflecting a combination of biological and technical factors (see Discussion). Moreover, we found that the expression of cdNTR genes was significantly positively associated with the core characteristics of PCa luminal epithelial cells (Figure 2F), including cell identity, malignant transformation and invasion, and proliferation and cell cycle, supporting their functional roles in PCa. For example, the expression of adrenergic receptors ADRB1 and ADRB2 was broadly positively correlated with hallmarks such as androgen response, epithelial-mesenchymal transition, P53 pathway, hypoxia, mTORC1 signaling, and PI3K-AKT-mTOR signaling, consistent with previous findings that βadrenergic receptor signaling activates oncogenic pathways and contributes to PCa progression (7,27). Other associations were also observed, such as the positive correlation between CHRM3 expression and P53 pathway, and that between CHRNA2 expression and fatty acid metabolism. Collectively, these analyses demonstrate that cdNTR genes exhibit pronounced transcriptomic upregulation at single-cell resolution, particularly in luminal epithelial cells, and show strong concordance with bulk transcriptomic data, supporting their potential functional relevance in PCa progression.

Development and validation of the PCaSig using the elastic-net algorithm

Among the 13 cdNTR genes, seven BCR-related genes were identified in the TCGA-PRAD cohort using univariate Cox regression analysis (Figure 3A). Specifically, P2RX5, ADRA2A and GABRE were associated with HRs >1, whereas CHRM3, P2RX1, ADRA1A and CHRNA2 showed HRs <1. Subsequently, the elastic-net algorithm, which has been shown to be more robust than other algorithms was applied to remove redundancy and select the most informative prognostic markers for PCa. Under the minimum criteria (Figure 3B), all seven BCR-related genes were retained, enabling the construction of a seven-gene BCR prognostic signature (PCaSig). The risk score for each patient was calculated as follows: risk score = (0.053 × P2RX5 expression level) + (0.121 × ADRA2A expression level) + (0.185 × GABRE expression level) + (−0.103 × CHRM3 expression level) + (−0.108 × P2RX1 expression level) + (−0.112 × ADRA1A expression level) + (−0.148 × CHRNA2 expression level).

Figure 3 Performance of the PCaSig. (A) Forest plot showing seven genes significantly associated with BCR based on univariate Cox analysis. (B) Diagram of 10-fold cross-validated error under the value of tuning parameter lambda in the LASSO Cox regression model. The red dashed vertical line indicates the optimal lambda value that minimized the mean error. (C) Kaplan-Meier survival curves of the PCaSig in the discovery cohort for predicting BCR-free survival. (D) Time-dependent ROC curves for 1-, 3-, and 5-year BCR prediction based on PCaSig in the discovery cohort. (E) Kaplan-Meier survival curves of PCaSig in the validation cohort DKFZ-PRAD for predicting BCR-free survival. (F) Time-dependent ROC curve for 1-, 3-, and 5-year BCR prediction based on PCaSig in the validation cohort DKFZ-PRAD. (G) Kaplan-Meier survival curves of PCaSig in the second validation cohort GSE54460 for predicting BCR-free survival. (H) Time-dependent ROC curve for 1-, 3-, and 5-year BCR prediction based on PCaSig in the validation cohort GSE54460. AUC, area under the curve; BCR, biochemical recurrence; CI, confidence interval; C-index, concordance index; FPR, false positive rate; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; PCaSig, prostate cancer signature; ROC, receiver operating characteristic; TPR, true positive rate.

Patients were stratified into high and low-risk groups based on the median risk score. This stratification was significantly associated with BCR, with an HR of 3.35 [95% confidence interval (CI): 2.15–5.22; P<0.001; Figure 3C]. The prognostic performance of PCaSig yielded AUC values of 0.67 (95% CI: 0.59–0.75), 0.73 (95% CI: 0.66–0.79), and 0.71 (95% CI: 0.62–0.80) for 1-, 3-, and 5-year BCR-free survival, respectively (Figure 3D), with a C-index of 0.68 (95% CI: 0.63–0.74). In the DKFZ-PRAD validation cohort, PCaSig consistently identified patients at high risk of BCR (Figure 3E), achieving AUC values of 0.83 (95% CI: 0.74–0.93), 0.83 (95% CI: 0.73–0.93), and 0.88 (95% CI: 0.77–0.99) for 1-, 3-, and 5-year BCR-free survival (Figure 3F), and a C-index of 0.76 (95% CI: 0.67–0.84). Similar results were also observed in the GSE54460 validation cohort (Figures 3G,3H). In a meta-analysis integrating all training and validation cohorts, PCaSig showed a strong association with outcome (P=1.18×10⁻⁵⁹; HR =3.90; 95% CI: 3.31–4.60). Importantly, PCaSig remained independently associated with BCR in multivariate Cox proportional hazards models after adjustment for age, GS, PSA level, pathologic T stage, pathologic N stage, and residual tumor (P=0.002; Table S5). Moreover, we compared PCaSig with 15 previously published multigene prognostic signatures associated with diverse biological processes (Table S3), including ferroptosis, metabolism, cellular senescence, epithelial-mesenchymal transition, metastasis, and cancer-associated fibroblasts. The results demonstrated that PCaSig exhibited superior or comparable performance to the other 15 signatures in predicting 1-, 3-, and 5-year BCR-free survival. In particular, PCaSig achieved the best predictive accuracy for 3- and 5-year BCR-free survival (Table S6). Notably, 10 of the 15 comparator signatures were also developed using TCGA-PRAD as the training cohort, further supporting the robustness of the comparison. Collectively, these results demonstrate the robust prognostic performance, reproducibility, and concordance of PCaSig across multiple independent cohorts.

PCaSig correlates with aggressive clinical and molecular features

Exploring the clinicopathological and biological underpinnings of the PCaSig, we found that patients with high PCaSig scores exhibited a higher proportion of residual tumor (P=5.44×10⁻⁵; Figure 4A), more advanced pathologic T stage (P=1.62×10⁻⁷; Figure 4B), positive N stage (P=6.59×10⁻⁴; Figure 4C), elevated PSA levels (P=0.02; Figure 4D), and higher GS (P=1.11×10⁻¹¹; Figure 4E). Consistent results were observed in the DKFZ-PRAD validation cohort, in which high PCaSig scores were also associated with advanced pathologic T stage (P=2.17×10⁻⁸; Figure 4F), elevated PSA levels (P=2.23×10⁻⁶; Figure 4G), and higher GS (P=6.39×10⁻⁶; Figure 4H). Collectively, these results indicate that aggressive clinical characteristics, including advanced tumor stage and elevated PSA levels, are more prevalent in PCa patients with high PCaSig scores.

Figure 4 Clinicopathological and molecular characteristics of PCa stratified by PCaSig risk scores. (A-C) Bar charts showing the distribution of residual tumor status (A), pathological T stage (B), and N stage (C) in the discovery cohort. (D) Bar chart showing the distribution of pathological T stage in the validation cohort DKFZ-PRAD. (E,F) Box plots comparing PSA levels (E) and GSs (F) between high- and low-PCaSig groups in the discovery cohort. (G,H) Box plots comparing PSA levels (G) and GSs (H) between high- and low-PCaSig groups in the validation cohort DKFZ-PRAD. (I) Gene set enrichment analysis illustrating significantly upregulated and downregulated pathways in tumors with high PCaSig scores. (J) Box plots comparing representative pathway activities between high- and low-PCaSig groups. The center line represents the median, and the box bounds indicate the 25th and 75th percentiles. FDR, false discovery rate; GS, Gleason score; N, node; PCa, prostate cancer; PCaSig, prostate cancer signature; PSA, prostate-specific antigen; T, tumor.

GSEA showed that tumors with high PCaSig scores exhibited significant suppression of neurotransmitter-related signaling pathways (Figure 4I; Table S7), including cholinergic synapse [normalized enrichment score (NES) =−1.62, FDR =7.58×10⁻³], dopaminergic synapse (NES =−1.51, FDR =1.40×10⁻²), axon guidance (NES=−1.81, FDR =4.41×10⁻⁵), and neuroactive ligand signaling (NES =−1.45, FDR =1.74×10⁻²). In contrast, these tumors showed activation of cell cycle-related and immune-related pathways (Figure 4I; Table S7), such as the cell cycle (NES =1.58, FDR =1.62×10⁻³), cytokine-cytokine receptor interaction (NES =1.51, FDR =2.15×10⁻³), and antigen processing and presentation (NES =1.74, FDR =4.72×10⁻³). Consistently, GSVA of tumor hallmark pathways further supported these findings (Figure 4J; Table S8), demonstrating significant activation of cell cycle-associated processes, including E2F targets (FDR =1.48×10⁻⁶) and the G2M checkpoint (FDR =5.50×10⁻⁴), whereas androgen response signaling was markedly suppressed (FDR =1.37×10⁻¹²). Together, these results indicate that PCa with high PCaSig scores is characterized by suppressed neurotransmitter-related signaling and enhanced activation of cell cycle-related and immune-related pathways.

Correlation of PCaSig risk scores with TMB and mutation landscapes

Somatic mutations have been shown to play a critical role in cancer evolution. We found that the PCaSig risk score was positively correlated with TMB (R=0.15, P=7.30×10⁻⁴; Figure 5A), and that patients in the high PCaSig-risk group exhibited significantly higher TMB levels (P=2.02×10⁻³; Figure 5B). Kaplan-Meier analysis demonstrated that patients with low TMB had more favorable BCR-free survival (P=0.02; Figure 5C). Furthermore, combined prognostic stratification revealed that patients with both low PCaSig risk and low TMB experienced the most favorable clinical outcomes (P<0.001; Figure 5D). Importantly, PCaSig further stratified PCa patients into distinct BCR risk groups regardless of TMB status, whereas TMB alone failed to discriminate BCR outcomes within either the high- or low-PCaSig risk groups. Given that TMB has been associated with immunotherapy response in PCa (28), these findings imply that incorporating PCaSig may provide complementary information beyond TMB for the assessment of immunotherapy-related characteristics. However, validation in independent immunotherapy-treated cohorts will be required to determine whether PCaSig has predictive value for immunotherapy response in PCa.

Figure 5 Mutant landscape of PCa. (A) Comparison of TMB between high and low PCaSig risk groups. (B) Correlation of PCaSig risk score with TMB. (C) KM curves for TMB stratification. (D) KM curves for TMB-PCaSig risk stratification. (E) Overall mutation profile. BCR, biochemical recurrence; KM, Kaplan-Meier; PCa, prostate cancer; PCaSig, prostate cancer signature; TMB, tumor mutation burden.

Comparison of the somatic mutation landscapes between PCaSig risk groups revealed that the most frequently mutated genes were largely consistent across groups (Figure 5E). Notably, mutations in TP53 and CTNNB1 were significantly enriched in the high PCaSig-risk group (FDR =4.42×10⁻³, FDR =8.53×10⁻⁴, respectively), with all ten patients with CTNNB1 mutations belonging to this high-risk group. The enrichment of TP53 and CTNNB1 mutations in PCa patients with high PCaSig scores, indicative of aggressive disease, aligns with previous studies linking these mutations to PCa progression (29,30). In addition, STAB2, PCDHB7 and LRP1B exhibited increased mutation rates in the high PCaSig-risk group (FDR =3.57×10⁻³, FDR =1.48×10⁻², FDR =2.32×10⁻², respectively; data not shown). In contrast, no mutations were significantly enriched in the low-risk group.

Immune microenvironment underpinnings of PCaSig

To elucidate the immune microenvironmental basis of PCaSig, we first estimated the relative abundance of infiltrating immune cells using the CIBERSORT algorithm (23), and analyzed their correlations with PCaSig risk scores. Significant differences in immune cell infiltration were observed between the high- and low-risk groups, particularly for CD8⁺ T cells, CD4⁺ T cells, B cells and mast cells (Figure 6A). Beyond immune cell infiltration, we further assessed the expression of ICMs. Multiple ICMs, including CD27, CTLA4, PDCD1, and TIGIT, exhibited markedly different expression levels between the two groups (Figure 6B). In addition, correlation analyses between the 13 cdNTR genes and ICMs revealed significant associations (Figure 6C). Specifically, CHRNA2 and ADRB1 were positively correlated with multiple ICMs, whereas P2RX5 and GABRP exhibited negative correlations. Other cdNTRs showed more complex interaction patterns. Collectively, these findings suggest that neural elements may infiltrate the tumor microenvironment and modulate antitumor immune responses through NTRs, consistent with previous observations reported by Winkler et al. (5).

Figure 6 Immune microenvironment characteristics associated with PCaSig. (A) Violin plots showing differences in immune cell infiltration between high- and low-PCaSig risk groups. (B) Violin plots showing differential expression of ICMs between PCaSig risk groups. (C) Heatmap depicting the correlation coefficients between the 13 consistently dysregulated NTR genes and the ICMs shown in (B). Statistical significance: *, P<0.05; **, P<0.01; ***, P<0.001. ICM, immune checkpoint molecule; NTR, neurotransmitter receptor; PCaSig, prostate cancer signature.

Incremental predictive value of PCaSig in a clinical nomogram

A clinical nomogram was initially constructed based on the GS and pathologic T stage (Figure 7A). The AUC values for predicting 1-, 3- and 5-year BCR-free survival were 0.74 (95% CI: 0.67–0.81; Figure 7B), 0.73 (95% CI: 0.67–0.80; Figure 7C) and 0.77 (95% CI: 0.69–0.85; Figure 7D), respectively, yielding a C-index of 0.71 (95% CI: 0.66–0.76). Subsequently, an integrated nomogram incorporating the PCaSig risk score together with GS and pathologic T stage was developed (Figure 7E). This integrated model demonstrated improved predictive performance, with AUC values of 0.75 (95% CI: 0.68–0.81; Figure 7B) for 1-year, 0.77 (95% CI: 0.71–0.83; Figure 7C) for 3-year, and 0.79 (95% CI: 0.71–0.87; Figure 7D) for 5-year BCR-free survival, resulting in a higher C-index of 0.78 (95% CI: 0.69–0.79). Time-dependent ROC analyses further demonstrated that the integrated nomogram achieved superior sensitivity and specificity compared with any single prognostic factor, including PCaSig, GS, and pathologic T stage, for predicting 1-year (Figure 7B, AUC: 0.67, 0.71, and 0.70, respectively), 3-year (Figure 7C, AUC: 0.73, 0.71, and 0.68, respectively) and 5-year (Figure 7D, AUC: 0.71, 0.71, and 0.74, respectively) BCR-free survival. Calibration curves demonstrated good agreement between the predicted and observed BCR-free survival probabilities (Figure 7F). Furthermore, decision curve analysis indicated that the integrated nomogram provided the greatest clinical net benefit across a wide range of threshold probabilities compared with GS, pathologic T stage, and the initial clinical nomogram (Figure 7G). Collectively, these results indicate that the nomogram integrating PCaSig with GS and pathologic T stage offers superior prognostic accuracy compared with either the clinical nomogram alone or any individual clinicopathological factor.

Figure 7 Nomogram construction and performance evaluation. (A) Clinical nomogram constructed based on pathological T stage and GS for predicting 1-, 3-, and 5-year BCR-free survival. (B-D) Time-dependent ROC curves evaluating the predictive performance of pathological T stage, GS, PCaSig, the clinical nomogram, and the integrated nomogram (PCaSig combined with GS and pathological T stage) for 1- (B), 3- (C), and 5-year (D) BCR-free survival. (E) Integrated nomogram incorporating PCaSig, GS and T stage for predicting 1-, 3-, and 5-year BCR-free survival. (F) Calibration curves assessing the accuracy of the integrated nomogram for predicting 1-, 3-, and 5-year BCR-free survival. The dotted line represents the ideal reference, while the violet, green, and pink solid lines represent the observed calibration curves for 1-, 3-, and 5-year BCR-free survival predictions, respectively. (G) Decision curve analysis assesses the clinical utility of the integrated nomogram and clinical factors for 3-year risk assessment. BCR, biochemical recurrence; FPR, false positive rate; GS, Gleason score; PCaSig, prostate cancer signature; ROC, receiver operating characteristic; T, tumor; TPR, true positive rate.

Discussion

Accurate prediction of BCR is essential for optimizing postoperative management and reducing PCa-related mortality among patients who experience recurrence following radical prostatectomy. Although numerous transcriptomic biomarkers associated with PCa progression and prognosis have been reported, none have yet been translated into routine clinical practice, underscoring the need for more robust and biologically informative predictive models.

The nervous system has historically received relatively little attention in cancer research. However, with advancements in tumor biology, it is increasingly recognized as playing important roles in tumor initiation and progression (4). NTRs, key membrane proteins mediating interactions between nerves and cancer cells, are emerging as potential molecular targets for cancer therapy (31). While, the spectrum of NTR dysregulation has been investigated in colorectal cancer (CRC) (32), brain tumor (33), and hepatocellular carcinoma (HCC) (34), it has not been systematically analyzed in PCa. In this study, using paired tumor-normal samples analyses in discovery and validation datasets, we identified 13 cdNTR genes, including four adrenergic receptor genes (ADRB1, ADRB2, ADRA2A, and ADRA1A), three cholinergic receptor genes (CHRNA5, CHRNA2, and CHRM3), two purinergic receptor genes (P2RX1 and P2RX5), and four GABA receptor genes (GABRB3, GABRG3, GABRE, and GABRP). Adrenergic signaling represents one of the best-characterized neural regulatory mechanisms in PCa (35). Sympathetic nerve-derived catecholamines activate adrenergic receptors and promote tumor growth, angiogenesis, invasion, metastasis, and treatment resistance (26). In particular, β-adrenergic receptor signaling has been shown to drive PCa progression through activation of Sonic Hedgehog-Gli1 signaling and other oncogenic pathways (27). Consistently, three adrenergic receptor genes (ADRB1, ADRB2, and ADRA2A) were upregulated in our study, suggesting enhanced sensitivity of tumor cells to sympathetic neural inputs. Cholinergic signaling has also been implicated in PCa progression. Previous studies have demonstrated that CHRM3 activation promotes tumor growth and castration resistance through CaM/CaMKK-Akt and FAK-YAP signaling pathways (36,37). Moreover, autonomic nerve-derived acetylcholine can directly regulate prostate epithelial and tumor cell behavior (38). The upregulation of CHRM3, CHRNA5, and CHRNA2 observed in our study further supports the importance of cholinergic signaling in PCa biology. Purinergic receptors function as sensors of extracellular ATP and mediate communication between stressed cells, immune cells, and tumor cells (39). Growing evidence suggests that purinergic signaling regulates tumor proliferation, inflammation, and antitumor immunity (40). Interestingly, we observed opposite dysregulation patterns of the P2RX family members, with P2RX5 being upregulated and P2RX1 being downregulated, implying receptor-specific functions of purinergic signaling during PCa progression. Finally, accumulating evidence indicates that GABAergic signaling participates in cancer development through the regulation of proliferation, migration, invasion, and metabolic adaptation (41). In our study, GABRB3 and GABRG3 were upregulated, whereas GABRE and GABRP were downregulated, suggesting substantial remodeling of GABA receptor composition in PCa. Such alterations may influence both tumor cell behavior and neuroepithelial communication within the tumor microenvironment. In addition, cdNTR genes showed significant associations with hallmark programs that define the biological characteristics of malignant prostate luminal epithelial cells. In particular, adrenergic receptor genes ADRB1 and ADRB2 were positively correlated with androgen response, epithelial-mesenchymal transition, PI3K-AKT-mTOR signaling, hypoxia, and P53 pathway activity, suggesting potential links between neural signaling and oncogenic processes in PCa. Similar associations were observed for CHRM3 and CHRNA2. These findings indicate that dysregulated NTR genes may not simply represent lineage-specific markers but may actively participate in biological programs associated with tumor progression, providing additional mechanistic support for their inclusion in PCaSig.

Comparing 13 identified cdNTR genes with the 29 NTR-related DEGs reported by Zhang et al. in CRC (10 upregulated and 19 downregulated genes) (32), we found that only CHRNA5 and ADRB1 were shared and consistently dysregulated in the same direction, both showing increased expression in tumor tissues. Four additional cdNTR genes (GABRE, GABRP, ADRA2A, and ADRB2) were also reported by Zhang et al.; however, they exhibited opposite directions of dysregulation. Similarly, Belotti et al. identified 50 differentially expressed NTR genes in brain tumors, including low-grade glioma and glioblastoma multiforme (33). Among these genes, only CHRNA5 overlapped with our findings and showed increased expression in glioblastoma multiforme. Six additional cdNTR genes (ADRA2A, ADRB1, CHRM3, CHRNA2, GABRB3, and GABRG3) were also reported as dysregulated in brain tumors but displayed expression patterns opposite to those observed in our study. In HCC, Wang et al. identified seven NTR genes with increased expression and one NTR gene, GRIN2B, with decreased expression in high-risk HCC tumors (34); however, none of these genes overlapped with our identified cdNTR genes. Notably, we also observed limited overlap in dysregulated NTR genes among the studies by Zhang et al., Belotti et al., and Wang et al. Detailed information on the dysregulated NTR genes reported in these studies is provided in Table S9. Overall, the limited overlap and inconsistent dysregulation patterns of NTR genes across cancer types suggest that NTR-mediated signaling may have diverse and context-dependent functions in tumor development and progression. This observation is consistent with a recent study reporting distinct and even opposite prognostic implications of neuroregulatory subtypes across different cancers (42). These observations further highlight the necessity of conducting cancer-type-specific investigations of NTR dysregulation and support the rationale for a systematic characterization of neuroregulatory alterations in PCa.

Several factors may explain the discordance between bulk and single-cell transcriptomic analyses, particularly for genes identified as downregulated in tumors. Bulk RNA-seq captures averaged expression signals across all cells in a tissue, integrating both transcriptional regulation within individual cell types and changes in cellular composition. In contrast, scRNA-seq enables cell-type resolved analysis but is inherently limited by sparse gene detection and reduced sensitivity, especially for genes expressed in a small fraction of cells (43). In our study, the 13 cdNTR genes exhibited strong cell-type specificity, with most being expressed in only one or two cell types and in a relatively small proportion of cells. The nine genes upregulated in tumors were consistently validated at the single-cell level, suggesting that their increased bulk expression reflects genuine transcriptional activation within specific tumor-associated cell populations. However, the four genes downregulated in bulk data did not display sufficiently robust cell-type specific downregulation. This discrepancy likely reflects a combination of biological and technical factors. Apparent downregulation in bulk tissue may primarily result from a reduced abundance of the cell populations expressing these genes, rather than decreased expression per cell. Additionally, low expression levels and dropout events in scRNA-seq reduce statistical power to detect subtle decreases in gene expression. Finally, increased transcriptional heterogeneity within tumor cells may further obscure consistent downregulation signals at the single-cell level. These observations highlight the importance of integrative analyses combining bulk and single-cell data and underscore the need for cautious interpretation of cell-type specific downregulation inferred solely from bulk transcriptomic analyses.

Currently, prediction of BCR in PCa relies primarily on clinicopathological parameters, including PSA, GS, and pathological T stage. However, these parameters are insufficient for precise risk stratification and guiding treatment decisions. In this study, we developed PCaSig, a transcriptomic signature that was significantly associated with unfavorable BCR outcomes and aggressive clinical and molecular features, in both the discovery and validation cohorts. Importantly, PCaSig signature could be used to stratify PCa patients with distinct BCR outcomes, regardless of whether they fell into the high or low TMB groups, suggesting that PCaSig provides prognostic information beyond TMB alone. Given its significant associations with immune-related characteristics, PCaSig may also have potential relevance for immunotherapy-related patient stratification. However, immunotherapy-treated cohorts and prospective studies will be required to determine whether PCaSig has predictive value for immunotherapy response, either alone or in combination with TMB. We further established a well-calibrated predictive nomogram that integrated PCaSig with GS and pathological T stage. While the incorporation of PCaSig led to only modest improvements in time-dependent AUC values, it substantially improved the overall discriminative ability of the model, as reflected by an increased C-index (from 0.71 to 0.78), and additionally demonstrated greater clinical net benefit in decision curve analysis. These findings suggest that PCaSig provides prognostic information complementary to established clinicopathological factors rather than replacing them. Clinically, PCaSig may be particularly valuable for refining risk stratification among patients with similar GSs and pathological stages but differing risks of BCR. Such molecular stratification could potentially assist in postoperative surveillance planning and patient selection for early intervention strategies. However, prospective validation and clinical utility studies will be required before incorporation of PCaSig into routine clinical practice.

Notably, the signature assessed in this study has certain limitations regarding its clinical application. First, the construction of PCaSig was based on PCa 495 patients from TCGA, and validated in two independent cohorts involving 215 patients. Further evaluation in larger, independent PCa cohorts is warranted to confirm its robustness as a predictor of BCR. Second, the predictive accuracy of this model was assessed using retrospective data; thus, prospective, multicenter clinical studies are required to validate its clinical utility. In addition, the current model is only applicable to patients undergoing radical prostatectomy and may not be directly generalizable to those receiving radical radiotherapy. Beyond these clinical limitations, this study also has biological constraints. Experimental studies are needed to confirm the expression of the dysregulated NTR genes identified here and to validate their protein-level expression in PCa tissues. Moreover, future in vitro and in vivo investigations are necessary to functionally characterize the regulatory roles of genes within PCaSig, particularly in the context of neuron-cancer interactions.

Conclusions

In summary, to the best of our knowledge, this study is the first to systematically characterize NTR transcript dysregulation in PCa at both bulk and single cell resolution, and to investigate its association with BCR, leading to the development of a machine-learning-based BCR risk prediction signature, PCaSig. High PCaSig scores were significantly associated with unfavorable BCR outcomes as well as aggressive clinical and molecular features across multiple cohorts. Our findings further suggest that neural signaling may influence tumor immunity through NTR-mediated mechanisms, highlighting a potential link between NTR-mediated signaling and the tumor immune microenvironment. The integrated nomogram combining PCaSig with GS and pathological T stage provides improved BCR risk prediction and may facilitate more personalized therapeutic decision-making for patients with PCa.

Acknowledgments

We extend our gratitude to all the researchers whose contributions to the public dataset were essential for this study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0461/rc

Peer Review File: Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0461/prf

Funding: This work was supported by the Fujian Province Natural Science Foundation (No. 2025J08060 to S.L.), the Hong Kong Scholars Program (No. XJ2023053 to Y.Z.), the China Postdoctoral Science Foundation (No. 2023T160175 to Y.Z. and No. 2021MD703831 to Y.Z.), the Heilongjiang Postdoctoral Young Talents Program (No. LBH-TZ2114 to Y.Z.), the 2025 Undergraduate Education and Teaching Reform Research Project of Fuzhou University (No. JG36250035 to S.L.), and the Fujian Provincial Joint Fund for Scientific and Technological Innovation Project (No. 2024Y9623 to J.C.).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0461/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Kratzer TB, Wagle NS, et al. Cancer statistics, 2026. CA Cancer J Clin 2026;76:e70043. [Crossref] [PubMed]
2. Wei JT, Barocas D, Carlsson S, et al. Early Detection of Prostate Cancer: AUA/SUO Guideline Part I: Prostate Cancer Screening. J Urol 2023;210:46-53. [Crossref] [PubMed]
3. Van den Broeck T, van den Bergh RCN, Arfi N, et al. Prognostic Value of Biochemical Recurrence Following Treatment with Curative Intent for Prostate Cancer: A Systematic Review. Eur Urol 2019;75:967-87. [Crossref] [PubMed]
4. Monje M, Borniger JC, D'Silva NJ, et al. Roadmap for the Emerging Field of Cancer Neuroscience. Cell 2020;181:219-22. [Crossref] [PubMed]
5. Winkler F, Venkatesh HS, Amit M, et al. Cancer neuroscience: State of the field, emerging directions. Cell 2023;186:1689-707. [Crossref] [PubMed]
6. Hanahan D, Monje M. Cancer hallmarks intersect with neuroscience in the tumor microenvironment. Cancer Cell 2023;41:573-80. [Crossref] [PubMed]
7. Li TJ, Jiang J, Tang YL, et al. Insights into the leveraging of GABAergic signaling in cancer therapy. Cancer Med 2023;12:14498-510. [Crossref] [PubMed]
8. Xia S, He C, Zhu Y, et al. GABA(B)R-Induced EGFR Transactivation Promotes Migration of Human Prostate Cancer Cells. Mol Pharmacol 2017;92:265-77. [Crossref] [PubMed]
9. Palamiuc L, Emerling BM. PSMA brings new flavors to PI3K signaling: A role for glutamate in prostate cancer. J Exp Med 2018;215:17-9. [Crossref] [PubMed]
10. Shore ND, Moul JW, Pienta KJ, et al. Biochemical recurrence in patients with prostate cancer after primary definitive therapy: treatment based on risk stratification. Prostate Cancer Prostatic Dis 2024;27:192-201. [Crossref] [PubMed]
11. Tilki D, Preisser F, Graefen M, et al. External Validation of the European Association of Urology Biochemical Recurrence Risk Groups to Predict Metastasis and Mortality After Radical Prostatectomy in a European Cohort. Eur Urol 2019;75:896-900. [Crossref] [PubMed]
12. Li R, Qu H, Wang S, et al. GDCRNATools: an R/Bioconductor package for integrative analysis of lncRNA, miRNA and mRNA data in GDC. Bioinformatics 2018;34:2515-7. [Crossref] [PubMed]
13. Chen Y, Chen L, Lun ATL, et al. edgeR v4: powerful differential analysis of sequencing data with expanded functionality and improved support for small counts and larger datasets. Nucleic Acids Res 2025;53:gkaf018. [Crossref] [PubMed]
14. Carvalho BS, Irizarry RA. A framework for oligonucleotide microarray preprocessing. Bioinformatics 2010;26:2363-7. [Crossref] [PubMed]
15. Henry GH, Malewska A, Joseph DB, et al. A Cellular Anatomy of the Normal Adult Human Prostate and Prostatic Urethra. Cell Rep 2018;25:3530-3542.e5. [Crossref] [PubMed]
16. Tuong ZK, Loudon KW, Berry B, et al. Resolving the immune landscape of human prostate at a single-cell level in health and cancer. Cell Rep 2021;37:110132. [Crossref] [PubMed]
17. Song H, Weinstein HNW, Allegakoen P, et al. Single-cell analysis of human primary prostate cancer reveals the heterogeneity of tumor-associated epithelial cell states. Nat Commun 2022;13:141. [Crossref] [PubMed]
18. Chen S, Zhu G, Yang Y, et al. Single-cell analysis reveals transcriptomic remodellings in distinct cell types that contribute to human prostate cancer progression. Nat Cell Biol 2021;23:87-98. [Crossref] [PubMed]
19. van den Brink SC, Sage F, Vértesy Á, et al. Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations. Nat Methods 2017;14:935-6. [Crossref] [PubMed]
20. Borcherding N, Vishwakarma A, Voigt AP, et al. Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun Biol 2021;4:122. [Crossref] [PubMed]
21. Andreatta M, Carmona SJ. UCell: Robust and scalable single-cell gene signature scoring. Comput Struct Biotechnol J 2021;19:3796-8. [Crossref] [PubMed]
22. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
23. Chen B, Khodadoust MS, Liu CL, et al. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol 2018;1711:243-59. [Crossref] [PubMed]
24. Thorsson V, Gibbs DL, Brown SD, et al. The Immune Landscape of Cancer. Immunity 2018;48:812-830.e14. [Crossref] [PubMed]
25. Kulik G. ADRB2-Targeting Therapies for Prostate Cancer. Cancers (Basel) 2019;11:358. [Crossref] [PubMed]
26. Gazova S, Klena L, Galvankova K, et al. Role of adrenergic receptors and their blocking in cancer research. Biomed Pharmacother 2025;192:118637. [Crossref] [PubMed]
27. Zhang M, Wang Q, Sun X, et al. β2 -adrenergic receptor signaling drives prostate cancer progression by targeting the Sonic hedgehog-Gli1 signaling activation. Prostate 2020;80:1328-40. [Crossref] [PubMed]
28. Zang PD, Chawla NS, Barragan-Carrillo R, et al. Tumor Mutational Burden in Metastatic Castration-Resistant Prostate Cancer and Response to Checkpoint Inhibition. JAMA Oncol 2024;10:531-2. [Crossref] [PubMed]
29. Nientiedt C, Budczies J, Endris V, et al. Mutations in TP53 or DNA damage repair genes define poor prognostic subgroups in primary prostate cancer. Urol Oncol 2022;40:8.e11-8.
30. Yeh Y, Guo Q, Connelly Z, et al. Wnt/Beta-Catenin Signaling and Prostate Cancer Therapy Resistance. Adv Exp Med Biol 2019;1210:351-78. [Crossref] [PubMed]
31. Jiang SH, Hu LP, Wang X, et al. Neurotransmitters: emerging targets in cancer. Oncogene 2020;39:503-15. [Crossref] [PubMed]
32. Zhang L, Deng Y, Yang J, et al. Neurotransmitter receptor-related gene signature as potential prognostic and therapeutic biomarkers in colorectal cancer. Front Cell Dev Biol 2023;11:1202193. [Crossref] [PubMed]
33. Belotti Y, Tolomeo S, Yu R, et al. Prognostic Neurotransmitter Receptors Genes Are Associated with Immune Response, Inflammation and Cancer Hallmarks in Brain Tumors. Cancers (Basel) 2022;14:2544. [Crossref] [PubMed]
34. Wang X, Li Y, Shi Y, et al. Comprehensive analysis to identify the neurotransmitter receptor-related genes as prognostic and therapeutic biomarkers in hepatocellular carcinoma. Front Cell Dev Biol 2022;10:887076. [Crossref] [PubMed]
35. Wang X, Shi M, Tian J, et al. Cancer and neurotransmitter receptors. Chin Med J (Engl) 2025;138:1540-58. [Crossref] [PubMed]
36. Wang N, Yao M, Xu J, et al. Autocrine Activation of CHRM3 Promotes Prostate Cancer Growth and Castration Resistance via CaM/CaMKK-Mediated Phosphorylation of Akt. Clin Cancer Res 2015;21:4676-85. [Crossref] [PubMed]
37. Goto Y, Ando T, Izumi H, et al. Muscarinic receptors promote castration-resistant growth of prostate cancer through a FAK-YAP signaling axis. Oncogene 2020;39:4014-27. [Crossref] [PubMed]
38. Liu Z, Peng Q, Wang Y, et al. Neuroscience in prostate cancer. Prostate Cancer Prostatic Dis 2026;29:478-87. [Crossref] [PubMed]
39. Wang X, Jiang SH, Ma M, et al. P2 Purinergic Receptors in Tumor Immunity. Cancer Res 2025;85:3826-41. [Crossref] [PubMed]
40. Wang Z, Zhu S, Tan S, et al. The P2 purinoceptors in prostate cancer. Purinergic Signal 2023;19:255-63. [Crossref] [PubMed]
41. Yang Y, Ren L, Li W, et al. GABAergic signaling as a potential therapeutic target in cancers. Biomed Pharmacother 2023;161:114410. [Crossref] [PubMed]
42. Luo S, Qiao S, Liu L, et al. Pan-cancer neurotransmitter receptor alterations define neuroregulatory subtypes with prognostic significance. Cell Rep 2026;45:117340. [Crossref] [PubMed]
43. Fan J, Slowikowski K, Zhang F. Single-cell transcriptomics in cancer: computational challenges and opportunities. Exp Mol Med 2020;52:1452-65. [Crossref] [PubMed]