Abstract

Background

   Prostate adenocarcinoma (PRAD) is a leading cause of male mortality,
   with NK/T cell communication being key areas of the research.

Methods

   The Seurat package was utilized to normalize and reduce the
   dimensionality of the single-cell data, and CellMarker 2.0 was employed
   for cell annotation. CellChat was utilized to construct the
   ligand-receptor interaction network of cell subsets. Differentially
   expressed genes (DEGs) were filtered by the limma package. Univariate
   Cox and the LASSO regression in the glmnet package were used to obtain
   biomarkers and develop a risk model. The survminer package was used to
   calculate the optimal threshold for dividing patients into high-risk
   and low-risk groups, and then Kaplan-Meier (KM) survival analysis was
   performed. Single-sample GSEA (ssGSEA), TIMER, and ESTIMATE packages
   were employed for immune infiltration analysis. Pathway analysis was
   conducted for the low- and high-risk groups using GSEA. Immunotherapy
   responses were evaluated by adopting TIDE method. Additional cellular
   validation (quantitative real-time PCR, CCK-8, Transwell, and scratch
   assay) was implemented to confirm the effects of feature genes on PRAD.

Results

   Compared with the benign group, NK/T cells were the cell type with the
   greatest changes in the tumor group, and their communication intensity
   was relatively high among all cell types. A RiskScore model was
   constructed as follows:
   [MATH:
   <mrow><mn>0.579</mn><mo>*</mo><mi>F</mi><mi>O</mi><mi>X</mi><mi>S</mi><
   mi
   mathvariant="italic">1</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.345</mn>
   <mo>*</mo><mi>G</mi><mi>P</mi><mi>C</mi><mi
   mathvariant="italic">6</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.385</mn>
   <mo>*</mo><mi>I</mi><mi>S</mi><mi>Y</mi><mi>N</mi><mi>A</mi><mi
   mathvariant="italic">1</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.418</mn>
   <mo>*</mo><mi>I</mi><mi>T</mi><mi>G</mi><mi>A</mi><mi>X</mi><mtext> </m
   text><mo>+</mo><mtext> </mtext><mn>0.792</mn><mo>*</mo><mi>M</mi><mi>G<
   /mi><mi>A</mi><mi>T</mi><mi
   mathvariant="italic">4</mi><mi>B</mi><mtext> </mtext><mo>+</mo><mtext>
   </mtext><mn>0.368</mn><mo>*</mo><mi>P</mi><mi>R</mi><mi>R</mi><mi
   mathvariant="italic">7</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.458</mn>
   <mo>*</mo><mi>R</mi><mi>E</mi><mi>X</mi><mi>O</mi><mi
   mathvariant="italic">2</mi></mrow> :MATH]
   . Analysis of the differences between the two risk groups showed that
   the level of immune infiltration was higher in the high-risk group, and
   it was significantly enriched in immune-correlated pathways, while the
   low-risk group was mainly enriched in metabolism-related pathways. TIDE
   analysis indicated that the high-risk group had higher immune escape
   potential. The cellular validation assays have revealed the higher
   expression of seven biomarkers in PRAD groups. Further, ISYNA1
   knockdown inhibited the proliferation, migration, and invasion ability
   of PRAD cells.

Conclusion

   The current research reveals key communication genes in PRAD, offering
   new possibilities for the exploration of new therapeutic targets.

   Keywords: NK/T cell communication, tumor microenvironment, prostate
   cancer, immune escape, machine learning

1. Introduction

   Globally, prostate cancer (PCa) is a common cause of death among men.
   Statistics show that in 2022, there were 1,466,680 new PCa cases
   globally, accounting for 7.3% of new cancers and ranking fourth
   ([29]1). More than 95% of PCa cases are adenocarcinomas, with the
   majority originating from acini and a minority originating from ducts
   ([30]2). Its incidence is closely related to location and age ([31]3,
   [32]4). In the early stage, prostate adenocarcinoma (PRAD) has
   relatively mild symptoms, with the most common being dysuria and
   increased frequency of urination ([33]5). In the advanced stage, PRAD
   may present with urinary retention and back pain ([34]5). Currently,
   the main treatment method for PCa is androgen deprivation therapy, but
   it is difficult to completely cure metastatic castration-resistant
   prostate cancer (mCRPC) caused by cancer metastasis ([35]6). The
   research on immunotherapy for PRAD mainly focuses on Sipuleucel-T and
   monoclonal antibodies ([36]7), but experiments have found that their
   efficacy in PRAD patients is not yet significant ([37]8, [38]9). Given
   the long natural course of PRAD and the relatively high mortality rate
   of moderate or high-risk local, locally advanced or metastatic cancers
   ([39]10), it is necessary to establish a prognostic model related to
   PRAD, which helps in the stratified treatment of patients and the
   selection of medical regimens.

   The single-cell RNA sequencing (scRNA-seq) technique has offered
   unprecedented opportunities for studying intercellular interactions
   within the tumor microenvironment (TME) in recent years ([40]11),
   particularly in the research of immune cells such as T cells and
   natural killer (NK) cells. NK/T cells play pivotal roles in tumor
   immune evasion and immune surveillance ([41]12). Studies have found
   that combinations of NK-cell-based immunotherapy and T-cell-based
   immunotherapy may be beneficial for tumor immunotherapy ([42]13).
   Nguyen et al. suggest that further detection of tumour-infiltrating
   lymphocytes in prostate cancer may help to gain insight into the
   regulatory mechanisms of T cells in the TME ([43]14). In addition, NK/T
   cells interact with tumor cells and other immune cells through cellular
   communication ([44]15). Currently, the systematic inference of
   ligand-receptor-mediated intercellular communication has become a
   research hotspot, which is closely linked to tumorigenesis ([45]16,
   [46]17). Understanding these intercellular communication mechanisms can
   uncover new immune escape pathways and provide novel strategies for
   immunotherapy in PRAD.

   In this study, by analyzing single-cell data from PRAD, we identified
   the cell communication receptor-ligand genes between NK/T cells and
   other cell types, and constructed their interaction networks. Through
   bioinformatics methods, we screened for differential genes related to
   the communication receptor-ligands between NK/T cells and other cells.
   This study developed a risk model with these key genes and conducted
   prognostic analysis. Then, we performed in vitro experimental
   verification on the obtained biomarkers and evaluated the immune
   infiltration levels and immune therapy responses of different risk
   groups. These findings not only provide potential biomarkers for early
   detection and prognostic assessment of PRAD but also open up new
   possibilities for future stratified and personalized treatment
   developments in PRAD.

2. Material and methods

2.1. Collection and processing of data

   The gene expression data and related clinical data of PRAD patients
   were obtained from the The Cancer Genome Atlas (TCGA,
   [47]https://portal.gdc.cancer.gov) database ([48]18). The FPKM values
   of TCGA-PRAD RNA sequencing (RNA-Seq) data was transformed into TPM and
   then log2-converted. After retaining the samples with complete
   progression-free interval (PFI), 495 PRAD samples and 52 adjacent
   normal control samples were recruited. The expression data and clinical
   data (MSKCC, Cancer Cell 2010) of PRAD were extracted from the
   cBioPortal website ([49]https://www.cbioportal.org/). After screening,
   a sum of 131 tumor samples were included. From Gene Expression Omnibus
   (GEO, [50]https://www.ncbi.nlm.nih.gov/geo/), we sourced four benign
   and four PRAD single-cell samples in [51]GSE193337. The library
   construction platform was 10x Genomics, and the sequencing platform was
   Illumina HiSeq 4000.

2.2. ScRNA-seq analysis

   The Seurat object was created by using the CreateSeuratObject function
   in the Seurat package ([52]19). Cells with the number of retained genes
   ranging from 200 to 8000 and the quantity of mitochondrial genes less
   than 20% were reserved. Next, the NormalizeData function ([53]19) was
   used for normalization. After the principal component analysis (PCA)
   dimensionality reduction, batch effect was eliminated by the harmony
   package ([54]20). Then, the RunUMAP function ([55]19) was utilized for
   dimensionality reduction by uniform manifold approximation and
   projection (UMAP). Finally, cells were clustered by the FindNeighbors
   and FindClusters functions ([56]19) with the parameters of dims = 1:30
   and resolution = 0.1. According to the marker genes provided by the
   CellMarker2.0 database, the cell types were annotated ([57]21).

2.3. Cell communication

   CellChat ([58]22) was used to develop the ligand-receptor interaction
   network of cell subpopulations. The netVisual_bubble and
   netVisual_circle packages ([59]22) were utilized to display the bubble
   plots of receptors and ligands between NK/T cells and other cells and
   the number of communications between them.

2.4. Screening of differently expressed genes differentially expressed genes
(DEGs)

   Differential analysis between the tumor and control sample groups in
   the TCGA cohort was performed in the limma package ([60]23). The
   screening criteria were that |log 2fold change (FC)| > log2(1.5) and p
   < 0.05, thus obtaining the DEGs. Then, the intersection was taken
   between the mRNAs related to the receptor and ligand genes of the
   communication between NK/T cells and other cells (cor > 0.5, p < 0.01)
   and the TCGA DEGs to obtain the mRNAs related to the communication of
   NK/T cells.

2.5. Development of the risk model and validation

   The TCGA-PRAD samples were assigned randomly at the ratio of 5:5 into
   training set and test set. The survival package ([61]24) was utilized
   to conduct univariate Cox regression analysis to determine the
   prognostic relevance of the intersection genes and p < 0.05 was
   selected for filtering. Subsequently, the glmnet package ([62]25) was
   employed to carry out LASSO COX regression analysis and stepwise
   regression to further narrow down the gene scope. Through multivariate
   analysis, the key genes and their corresponding coefficients were
   calculated, and the RiskScore for patients was calculated based on the
   formula:
   [MATH:
   <mrow><mtext>RiskScore </mtext><mo>=</mo><mtext> Σβi </mtext><mo>×</mo>
   <mtext> Expi</mtext></mrow> :MATH]

   where βi refers to the regression coefficient of each key gene, and
   Expi represents the expression of the corresponding gene. Then, the
   survminer package ([63]26) was employed to calculate the optimal
   threshold to classify the patients into low- and high-risk groups. KM
   survival analysis was performed, and receiver operating characteristic
   (ROC) curve model was constructed to predict the prognostic performance
   of the TCGA training set and test set. The area under the curve (AUC)
   value is a robust metric for assessing the classification accuracy of
   the RiskScore model. By comparing the AUC values to the assessment
   criteria, we can determine the effectiveness of the model in
   distinguishing between high- and low-risk patients. To better validate
   whether the model was robust, the same method was employed to the MSKCC
   dataset for verification.

2.6. Analysis of the tumor immune microenvironment (TIME)

   To analyze the association between RiskScore and the TIME, the
   infiltration status of immune cells was calculated using different
   methods. The ssGSEA function of GSVA ([64]27, [65]28) was utilized to
   compute the scores of 28 types of tumor infiltrating immune cells
   (TIICs) (1). The TIMER online tool ([66]http://cistrome.org/TIMER) was
   utilized to calculate six immune scores, and the ESTIMATE package
   ([67]29) was used to score the immune cells.

2.7. Pathway enrichment analysis

   The gene set enrichment analysis (GSEA) ([68]30) was employed for
   pathway analysis to explore the pathways of diverse biological
   processes in the two risk groups. The candidate gene sets from the
   Kyoto Encyclopedia of Genes and Genomes (KEGG) database were used to
   conduct GSEA ([69]31).

2.8. Evaluation of immunotherapy

   The TIDE method was adopted to evaluate immunotherapy response. The
   standardized transcriptome data were uploaded into the TIDE website
   ([70]http://tide.dfci.harvard.edu/) to calculate the TIDE score, with a
   higher TIDE prediction score indicating greater immune escape
   possibility and lower immunotherapy efficacy.

2.9. Cell culture

   In this study, the cell lines, namely human normal prostate epithelial
   cells PNT1A (cat.no RE59812) and human PCa cells DU145 (cat.no YS101C),
   were sourced from Shanghai Yaji Biotechnology Co., Ltd (Shanghai,
   China) ([71]http://www.yajimall.com/). All cells were cultured in
   Roswell Park Memorial Institute (RPMI) 1640 (11875093, Gibco, Grand
   Island, NY, USA) with the supplementation of 10% fetal bovine serum
   (FBS) (S9020, Solarbio Lifesciences, Beijing, China). All cells were
   incubated in the incubator at 37°C with 5% CO[2]. The cells underwent
   short tandem repeat (STR) identification, and the results of mycoplasma
   detection for these cells were found to be negative.

2.10. Quantitative real-time PCR

   Following the instructions, total RNA was isolated from PNT1A and DU145
   cells utilizing the TriZol total RNA extraction kit (15596026CN,
   Invitrogen, Carlsbad, CA, USA). Subsequently, the concentration of the
   isolated RNA was measured. Then, complementary DNA was synthesized by
   reverse transcription with a relevant assay kit (D7178S, Beyotime,
   Shanghai, China). After that, SYBR Green qPCR Mix (D7260, Beyotime,
   Shanghai, China) was used for the PCR assay according to the protocols.
   The conditions for qPCR included an initial step at 94°C for 30
   seconds, succeeded by 40 cycles consisting of 5 seconds at 94°C and 30
   seconds at 60°C. Finally, the relative level was calculated by the
   2^-ΔΔCT method ([72]32), with GAPDH as a reference gene. Based on
   National Center for Biotechnology Information (NCBI) sequences, the
   qRT-PCR primers used in this study were designed using Primer Premier 6
   software. The primer sequences were presented in [73]Supplementary
   Table 1 .

2.11. Cell transfection

   For the liposome transfection, the small interfering RNA against ISYNA1
   (si-ISYNA1) and the negative control small interfering RNA (si-NC) were
   all purchased from Merck KGaA (Shanghai, China) and transfected into
   DU145 cells utilizing lipofectamine 2000 transfection reagent
   (11668027, Invitrogen, Carlsbad, CA, USA) as per the manuals. The
   sequences applied for the transfection were
   5’-GCACCCATCATGCTGGACCTA-3’(si-ISYNA1#1) and
   5’-CAGAAGAATGGTACAAATCCAAG-3’(si-ISYNA1#2).

2.12. Cell viability

   DU145 cells in the logarithmic growth phase were plated into a 96-well
   dish at a density of 1 × 10^4 cells per well and incubated at 37°C with
   5% CO[2] for durations of 0, 24, 48, or 72 hours. Following this, 10 μL
   of CCK-8 reagent was added to each well, and the samples were incubated
   at 37°C for 2 hours. For the construction of the CCK-8 curve,
   absorbance measurements were taken at 450 nm, which served as the
   y-axis, while time was plotted on the x-axis.

2.13. Cell invasion assay

   For the cell invasion assay, 1 × 10^5 DU145 cells were suspended in 200
   μL of serum-free medium and cultured in the upper chamber of the
   Transwell (3422, Corning, Inc., Corning, NY, USA) coated with Matrigel
   (C0372, Beyotime, China), while 700 μL of medium containing 10% bovine
   serum were supplemented to the lower Transwell chamber. After
   incubation, the chamber was taken out and gently rinsed with phosphate
   buffered saline (PBS) buffer (P1010, Solarbio Lifesciences, Beijing,
   China) to remove the non-invasive cells. Subsequently, 4%
   paraformaldehyde (P1110, Solarbio Lifesciences, Beijing, China) was
   employed for fixing the cells, which were dyed by 0.1% crystal violet
   (G1063, Solarbio Lifesciences, Beijing, China) at room temperature for
   20 minutes. The invasive cells were observed with an optical microscope
   (DP27, Olympus, Tokyo, Japan) in three randomly selected fields of view
   ([74]33). Finally, to ensure the accuracy and reproducibility of cell
   counting, we used ImageJ software for automated cell counting.

2.14. Cell migration assay

   The DU145 cells (5 × 10^5 cells/well), after being transfected, were
   grown in a 6-well plate containing media devoid of serum. Upon reaching
   full confluence, an artificial wound was introduced into the monolayer
   using a 200-μL sterile pipette tip. Following a 48-hour incubation
   period, the cells were photographed under an inverted optical
   microscope (DP27, Olympus, Japan). Subsequently, the percentage of
   wound closure (%) was calculated to reflect the migration ability of
   the PRAD cells, ensuring uniqueness in reporting ([75]34).

2.15. Statistical analysis

   All the statistical data were analyzed in R language (version 3.6.0).
   Experimental data were analyzed using GraphPad Prism 8.0 software. The
   Wilcoxon rank-sum test was used to compare the differences between
   two-group continuous variables, the Spearman method was employed to
   calculate the correlations, the log-rank test was utilized to compare
   the survival differences among patients in each grouping. Unpaired
   t-test and one-way analysis of variance were applied for the comparison
   on the experimental data, and p < 0.05 signified a statistical
   significance.

3. Results

3.1. Single-cell atlas of PRAD

   Single-cell clustering and annotation analysis were conducted on
   samples from benign and tumor tissues of PRAD, identifying a total of
   seven cell subpopulations: NK/T cells, Macrophage cells, Epithelial
   cells, Fibroblast cells, Mast cells, Endothelial cells, and B cells
   ([76] Figure 1A ). [77]Figures 1B, C depict the marker genes for each
   of these cell subpopulations. Among them, high expression of CD3D and
   NKG7 was observed in NK/T cells; EPCAM in epithelial cells; VWF in
   endothelial cells; CD163 in macrophages; COL1A1 and ACTA2 in
   fibroblasts; CD79A in B cells; and TPSAB1 and CPA3 genes in mast cells.
   Analysis of the distribution of various cell types between the benign
   and tumor groups revealed that, compared to the benign group, NK/T
   cells exhibited the highest proportional increase in the tumor group
   ([78] Figure 1D ).

Figure 1.

   [79]Figure 1
   [80]Open in a new tab

   Single-cell atlas of PRAD. (A) UMAP dimensionality reduction plot for
   single-cell clustering and annotation of benign and tumor samples from
   PRAD. (B) Bubble plot for annotating the expression of marker genes in
   subpopulations. (C) Violin plot illustrating the expression of marker
   genes. (D) Proportion of each cell subpopulation within each sample
   between the two groups.

3.2. Signal communication among various cell types

   Analysis of the number of communications among various cell types
   revealed that B cells had the least number of communications among all
   cell types ([81] Figure 2A ). Examination of the intensity of
   communications among cells showed that NK/T cells exhibited high
   communication intensity, indicating their active role in the TME ([82]
   Figure 2B ). To visually present this point more clearly, the
   communication intensity of each cell type was further detailed, and the
   results demonstrated that the communication intensity between NK/T
   cells and Macrophage cells, Epithelial cells, as well as Endothelial
   cells was particularly prominent ([83] Figure 2C ). These findings
   showed that NK/T cells played a pivotal role in immune regulation
   within the TME and have close interactions with other key cell types.
   [84]Figure 2D further illustrates the receptors and ligands involved in
   cellular communications.

Figure 2.

   [85]Figure 2
   [86]Open in a new tab

   Analysis of signaling communication between NK/T Cells and other cell
   types. (A) Number of communications between cells. (B, C) Intensity of
   communications between cells. (D) Ligand-receptor mediated cellular
   communications between NK/T cells and other cell types.

3.3. Identification of DEGs related to communication receptors and ligands

   Differential analysis was conducted between tumor and control sample
   groups in the TCGA cohort, resulting in 1951 downregulated genes and
   971 upregulated genes ([87] Figures 3A, B ). Next, these DEGs were
   intersected with mRNAs linked to communication receptors and ligands
   between NK/T cells and other cells (correlation coefficient > 0.5, p <
   0.01). This intersection identified a total of 2026 overlapping genes
   that are associated with communication between NK/T cells ([88]
   Figure 3C ).

Figure 3.

   [89]Figure 3
   [90]Open in a new tab

   Identification of DEGs. (A) Volcano plot on the distribution of DEGs
   between PRAD and control samples in TCGA. (B) Heatmap illustrating the
   expression of the top 50 DEGs in TCGA. (C) Venn diagram showing the
   intersection between mRNAs related to communication receptors and
   ligands between NK/T cells and other cells, as well as the DEGs in
   TCGA.

3.4. Model construction and validation in PRAD

   The TCGA-PRAD samples were assigned at the ratio of 5:5 into a training
   set and a testing set to construct the model. Using the data from the
   training set, the intersected genes were subjected to univariate Cox
   proportional hazards regression. Subsequently, LASSO COX and stepwise
   regression were applied to further shrink the range of genes, and seven
   genes independently associated with prognosis were screened out for
   constructing the RiskScore model ([91] Figures 4A, B ). The formula for
   this RiskScore model is ([92] Figure 4C ):

Figure 4.

   [93]Figure 4
   [94]Open in a new tab

   Construction and validation of the risk model. (A, B) Changes in the
   regression coefficients of gene features in the LASSO regression model
   and the optimal penalty parameter (λ) determined by cross-validation.
   (C) Genes and their coefficients in the risk model. (D-G) 1 - 5-year
   ROC curves of the risk model and differences in PFI between the
   high-risk and low-risk groups in the training set (D), test set (E),
   TCGA cohort (F), and MSKCC cohort (G). (H) Expression of seven
   prognostic genes in the TCGA cohort. (I) Expression of seven prognostic
   genes in the MSKCC cohort. ns denotes p > 0.05, **p < 0.01, ****p <
   0.0001. * means vs. High-risk.
   [MATH:
   <mrow><mtext>RiskScore </mtext><mo>=</mo><mtext> </mtext><mn>0.579</mn>
   <mtext> </mtext><mo>*</mo><mtext> </mtext><mi>F</mi><mi>O</mi><mi>X</mi
   ><mi>S</mi><mi
   mathvariant="italic">1</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.345</mn>
   <mtext> </mtext><mo>*</mo><mtext> </mtext><mi>G</mi><mi>P</mi><mi>C</mi
   ><mi
   mathvariant="italic">6</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.385</mn>
   <mtext> </mtext><mo>*</mo><mtext> </mtext><mi>I</mi><mi>S</mi><mi>Y</mi
   ><mi>N</mi><mi>A</mi><mi
   mathvariant="italic">1</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.418</mn>
   <mtext> </mtext><mo>*</mo><mtext> </mtext><mi>I</mi><mi>T</mi><mi>G</mi
   ><mi>A</mi><mi>X</mi><mtext> </mtext><mo>+</mo><mtext> </mtext><mn>0.79
   2</mn><mtext> </mtext><mo>*</mo><mtext> </mtext><mi>M</mi><mi>G</mi><mi
   >A</mi><mi>T</mi><mi
   mathvariant="italic">4</mi><mi>B</mi><mtext> </mtext><mo>+</mo><mtext>
   </mtext><mn>0.368</mn><mtext> </mtext><mo>*</mo><mtext> </mtext><mi>P</
   mi><mi>R</mi><mi>R</mi><mi
   mathvariant="italic">7</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.458</mn>
   <mtext> </mtext><mo>*</mo><mtext> </mtext><mi>R</mi><mi>E</mi><mi>X</mi
   ><mi>O</mi><mi mathvariant="italic">2</mi></mrow> :MATH]

   Patients were classified into high- and low-risk group by the optimal
   cut-off value of the RiskScore. The ROC curve demonstrated that the
   TCGA cohort, training set, and testing set all exhibited relatively
   high AUC values at various time points, which verified its good
   classification accuracy. Moreover, compared with those with a low-risk
   score, the PFI of patients with a high-risk score was significantly
   lower ([95] Figures 4D-F , p < 0.001). The same method was used to
   validate the MSKCC dataset, and the model also had a relatively high
   AUC value. There was a significant difference in prognosis between two
   risk groups ([96] Figure 4G ). By analyzing the expressions of
   prognostic genes in the TCGA cohort and the MSKCC cohort, it was found
   that, except for the REXO2 gene which showed no significant difference
   between the two risk groups in the MSKCC dataset, the expressions of
   the other genes in the low-risk group were all remarkably lower ([97]
   Figures 4H-I , p < 0.05).

3.5. Correlation between the RiskScore and TIME

   To analyze the relationship between RiskScore and TIME, different
   methods were used to calculate the infiltration of immune cells.
   Analyzing the results of immune cell infiltration assessment by
   ESTIMATE, it was found that the low-risk group had lower StromalScore,
   ImmuneScore and ESTIMATEScore than the high-risk group, indicating that
   the immune infiltration level was lower in the low-risk group ([98]
   Figure 5A , p < 0.05). Based on TIMER, the immune cell scores of the
   TCGA dataset were calculated, and the results showed that the scores of
   Neutrophil, B_cell, Dendritic, Macrophage, CD4_Tcell were all
   remarkably lower in the low-risk group ([99] Figure 5B , p < 0.01).
   Using ssGSEA functional analysis to score 28 types of immune cells
   ([100]35), it could be seen that the infiltration scores of multiple
   immune cells such as myeloid derived suppressor cells (MDSC), activated
   CD8 T cell, NKT cell, regulatory T cell (Treg), NK cell in the low-risk
   group were all significantly lower ([101] Figure 5C , p < 0.05).

Figure 5.

   [102]Figure 5
   [103]Open in a new tab

   Differences in immune infiltration between the high-risk and low-risk
   groups. Differences in ESTIMATE (A), Timer (B), and ssGSEA (C) scores
   between the high-risk and low-risk groups in the TCGA cohort. *p <
   0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. ns (abbreviation of
   "non-significant") means p > 0.05.

3.6. Differences in enriched pathways between high-risk and low-risk groups

   The results of KEGG pathway enrichment analysis demonstrated that the
   high-risk group was significantly enriched in primary immunodeficiency,
   malaria, complement and coagulation cascades, ECM−receptor interaction,
   staphylococcus aureus infection pathways ([104] Figure 6A ). The
   low-risk group was significantly enriched in arginine and proline
   metabolism, beta−Alanine metabolism, Butanoate metabolism, mineral
   absorption and valine, leucine and isoleucine degradation pathways
   ([105] Figure 6B ).

Figure 6.

   [106]Figure 6
   [107]Open in a new tab

   KEGG enrichment pathway analysis. (A, B) KEGG enrichment pathway
   analysis of high-risk and low-risk groups.

3.7. Characterization of immune escape potential and immunosuppression in
PRAD patients from different risk groups

   TIDE was employed to estimate the potential clinical immunotherapy
   benefit to the two risk groups. The results demonstrated that in the
   TCGA cohort, the high-risk group had the highest TIDE score, showing a
   higher likelihood of immune escape in the high-risk group and less
   immunotherapy benefit ([108] Figure 7 , p < 0.05). Moreover, in the
   high-risk group, the scores of dysfunction, exclusion, MDSC, and cancer
   associated fibroblast (CAF) were significantly higher than those in the
   low-risk group, which revealed the significant immunosuppressive
   features and complex mechanisms in the TME ([109] Figure 7 , p < 0.05).

Figure 7.

   [110]Figure 7
   [111]Open in a new tab

   Differences in TIDE scores between high-risk and low-risk groups in
   TCGA.

3.8. In vitro experiments to validate the function of key genes in PRAD cells
and the potential role of ISYNA1

   Subsequently, the roles of seven characteristic genes in PRAD were
   validated by in vitro experiments. The results showed that all six
   genes (GPC6, ISYNA1, ITGAX, MGAT4B, PRR7, and REXO2), except FOXS1,
   were significantly highly expressed in PRAD cells ([112] Figure 8A , p
   < 0.01). Since ISYNA1 is highly expressed in DU145 cells and to date,
   there are few studies on ISYNA1 in tumor research, especially PRAD. For
   this reason, we chose to silence ISYNA1 in order to validate its
   potential effect on PRAD cell development ([113] Figure 8B , p < 0.01).
   Subsequently, we observed a significant decrease in the proliferation,
   migration and invasion ability of DU145 cells after silencing ISYNA1
   relative to controls ([114] Figures 8C-E , p < 0.001). These results
   suggest the potential impact of our NK/T cell-related key gene-based
   screen on PRAD cell genesis and development.

Figure 8.

   [115]Figure 8
   [116]Open in a new tab

   Cell validation results. (A) Quantified mRNA level of seven feature
   genes in PRAD cell line DU145 and the control cell line PNT1A. (B) The
   qRT-PCR to verify the transfection efficiency of ISYNA1 between the
   si-negative control (si-NC) and si-ISYNA1 groups. (C) Assessment of the
   proliferation level of DU145 cells after ISYNA1 silencing based on the
   CCK-8 assay. (D) The results of Transwell assay on the invasion of
   DU145 cells after the silencing of ISYNA1. (E) The results of scratch
   assay on the migration of DU145 cells after the silencing of ISYNA1.
   **p < 0.01, ***p < 0.001, ****p < 0.0001, and ns stands for no
   significant difference.

4. Discussion

   According to global data in 2022, although the mortality rate of PRAD
   ranks only eighth, its incidence rate ranks fourth ([117]36), and it
   has a high recurrence rate, with about one-third of men experiencing
   recurrence after local treatment, ultimately leading to the generation
   of malignant cells ([118]37). NK cells are cytotoxic immune cells
   capable of killing cancer cells in the innate immune system ([119]38).
   Research has found that cellular communication functions importantly in
   the TME ([120]39), particularly involving NK/T cells, which are
   associated with the immune evasion mechanisms of cancers ([121]40). In
   this study, by analyzing PRAD-related single-cell data, the interaction
   between NK/T cells and other cells was identified, and a risk model
   related to communication receptors and ligands was established. These
   results provided new insights into the immunotherapy of PRAD and
   opening up new potential pathways for personalized cancer treatment.

   The risk model in this study includes a total of seven related
   molecular markers, namely FOXS1, GPC6, ISYNA1, ITGAX, MGAT4B, PRR7, and
   REXO2. As a type of DNA-binding proteins, Forkhead-box (FOX) family
   proteins participate in cell growth, embryogenesis, differentiation,
   and lifespan ([122]41). FOXS1 has a high expression in pan-cancers, and
   the gene is linked to a worse surviva. Its overexpression is related to
   high infiltration levels of tumor-associated macrophages (TAM) and
   Tregs ([123]42), which is consistent with the result of a higher
   infiltration level of Tregs in the high-risk group in this study. GPC6
   belongs to the glypican gene family, and its function depends on
   interactions with cytoplasmic and/or extracellular ligands ([124]43).
   Chen et al.’s research has confirmed that GPC2 in the GPC family can
   promote the progression of PRAD ([125]44). ISYNA1 encodes an enzyme
   essential for inositol biosynthesis. Silencing ISYNA1 can inhibit the
   growth of tumor cells ([126]45). Jia et al.’s research has demonstrated
   that it can serve as a biomarker related to pan-cancer immunomodulation
   prognosis ([127]46). The in vitro experiments in this study also
   confirmed that the knockout of ISYNA1 would influence both the
   migration and invasion abilities of PRAD cells. ITGAX usually functions
   as a receptor for the extracellular matrix and is associated with tumor
   angiogenesis ([128]47, [129]48). Williams et al.’s research has
   confirmed that it can act as a susceptibility gene for aggressive PRAD
   ([130]49). The MGAT4B protein belongs to the glycosyltransferase family
   and can recognize and bind to both donor and receptor substrates
   ([131]50). In liver cancer cells, studies have found that the
   expression of MGAT4B is upregulated and it is related to immune evasion
   ([132]51). PRR7 is a transmembrane adaptor protein (TRAP) and an
   important organizer and regulator of immune receptor-mediated signal
   transduction ([133]52). It mediates T cell receptor signaling by
   assembling the membrane-proximal signalosome ([134]53). REXO2 is a
   homologue of the prokaryotic oligoribonuclease existing in human
   mitochondria and cytoplasm ([135]54). Previous experiments have
   confirmed that REXO2 can serve as a biomarker for PRAD ([136]55). The
   above evidence demonstrates the possible roles that these cell
   communication ligand-related genes may play in cancer and also implies
   their potential as biomarkers for PRAD.

   Analysis on the correlation between the RiskScore and the immune
   microenvironment revealed that most of the immune scores in the
   high-risk group were higher. However, despite the higher scores of NK
   cells, the scores of Tregs in the high-risk group were also higher.
   Tregs are a type of T lymphocyte with immunosuppressive functions and
   functions crucially in the process of the immune system eliminating
   tumor cells as well as in immune self-tolerance ([137]56). Tregs play
   an immunosuppressive role in the TIME. Studies have confirmed that
   Tregs can promote tumor immune evasion through the activation of
   transforming growth factor-β (TGF-β) mediated by integrin αvβ8
   ([138]57, [139]58). Therefore, although there is a strong immune
   response in the high-risk group, the high scores of Tregs indicate that
   tumors may “escape” immune attacks by inducing an immunosuppressive
   environment. In addition, the biomarkers FOXS1 and PRR7 obtained in
   this study are both closely related to Treg signal transduction or
   infiltration ([140]42, [141]53), which further implies that high-risk
   patients affect the TME through the immune evasion mechanism.

   The results of pathway enrichment showed that the pathways in the
   low-risk group were mainly enriched in metabolic pathways, whereas the
   high-risk group were mainly enriched in immune-related pathways. This
   reflects that the high-risk group is linked to immune dysfunction,
   while the low-risk group is in a relatively healthy state. The TIDE
   analysis revealed that the exclusion score of the high-risk group was
   high, suggesting that effector immune cells have difficulty effectively
   infiltrating into the interior of the tumor, which may be related to
   physical barriers or chemical repellent factors in the TME. Meanwhile,
   the score of MDSC in the high-risk group was also relatively high. MDSC
   is a heterogeneous group of immature bone marrow cells with the
   potential to effectively target T cells and suppress autoimmune
   responses ([142]59, [143]60). Studies have already demonstrated that
   MDSC can influence the suppression of specific T cell responses in
   myeloma by inducing the development of Tregs in the TME ([144]61).
   Combined with the results of immune infiltration, the high score of
   MDSC indicates that these cells may suppress effector immune responses
   by influencing Tregs, thereby helping the tumor achieve immune evasion.
   The increase in the score of CAF further reflects the activity of CAFs.
   CAF activated by TGF-β is one of the most abundant stromal cell types
   in almost all solid tumors, and its special role is widely recognized
   ([145]62). Moreover, TGF-β is closely related to Tregs, NK cells, and
   immune evasion ([146]63, [147]64). In conclusion, these characteristics
   of the high-risk group jointly reveal the tumor’s strategy of enhancing
   its survival ability through multi-level immune suppression mechanisms,
   providing important implications for in-depth research on tumor
   immunotherapy.

   The limitations of this study require further consideration and
   refinement. Firstly, the sample size and representativeness of this
   study may not fully encapsulate the heterogeneity and diversity of PRAD
   within the general population, potentially limiting the
   generalizability of the study results. To address this issue, future
   research should include larger and more diversified cohorts of PRAD
   samples, encompassing different stages and subtypes of the disease as
   well as various ethnic groups. Additionally, the data verification in
   this study mainly relies on bioinformatics techniques and in vitro
   experiments, lacking the support of in vivo experiments. To verify our
   research results more comprehensively, the construction of relevant
   animal models is an essential step. Furthermore, the clinical
   application of the identified biomarkers and risk models needs to be
   further explored and validated in clinical trials.

5. Conclusion

   The present study successfully developed a risk model closely related
   to communication receptors and ligands, marking significant progress in
   exploring the immune mechanisms of PRAD. The establishment of this
   model not only deepened our understanding of the complex immune
   microenvironment of PRAD, but also provided strong theoretical support
   and practical tools for the formulation of stratified and personalized
   strategies for PRAD. Furthermore, we confirmed a crucial role of NK/T
   cells in the development of PRAD, offering new targets and ideas for
   the progression of novel immunotherapies. In the future, we anticipate
   enhancing the antitumor activity of NK/T cells by modulating the
   expression or function of communication receptors and ligands, or
   designing more precise targeted drugs to inhibit the immune escape
   mechanisms of tumor cells, thereby further improving the treatment
   outcomes of PRAD and bringing more effective therapy and better quality
   of life to patients.

5.1. Scope statement

   Compared with the benign group, NK/T cells were the cell type with the
   greatest changes in the tumor group, and their communication intensity
   was relatively high among all cell types. A RiskScore model was
   constructed as follows:
   [MATH:
   <mrow><mn>0.579</mn><mo>*</mo><mi>F</mi><mi>O</mi><mi>X</mi><mi>S</mi><
   mi
   mathvariant="italic">1</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.345</mn>
   <mo>*</mo><mi>G</mi><mi>P</mi><mi>C</mi><mi
   mathvariant="italic">6</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.385</mn>
   <mo>*</mo><mi>I</mi><mi>S</mi><mi>Y</mi><mi>N</mi><mi>A</mi><mi
   mathvariant="italic">1</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.418</mn>
   <mo>*</mo><mi>I</mi><mi>T</mi><mi>G</mi><mi>A</mi><mi>X</mi><mtext> </m
   text><mo>+</mo><mtext> </mtext><mn>0.792</mn><mo>*</mo><mi>M</mi><mi>G<
   /mi><mi>A</mi><mi>T</mi><mi
   mathvariant="italic">4</mi><mi>B</mi><mtext> </mtext><mo>+</mo><mtext>
   </mtext><mn>0.368</mn><mo>*</mo><mi>P</mi><mi>R</mi><mi>R</mi><mi
   mathvariant="italic">7</mi><mo> </mo><mo>+</mo><mo> </mo><mn>0.458</mn>
   <mo>*</mo><mi>R</mi><mi>E</mi><mi>X</mi><mi>O</mi><mi
   mathvariant="italic">2</mi></mrow> :MATH]
   . Analysis of the differences between the two risk groups showed that
   the level of immune infiltration was higher in the high-risk group, and
   it was significantly enriched in immune-correlated pathways, while the
   low-risk group was mainly enriched in metabolism-related pathways. TIDE
   analysis indicated that the high-risk group had higher immune escape
   potential. The cellular validation assays have revealed the higher
   expression of seven biomarkers in PRAD groups. Further, ISYNA1
   knockdown inhibited the migratory and invasive capabilities of PRAD
   cells. The current research reveals key communication genes in PRAD,
   offering new possibilities for the exploration of new therapeutic
   targets.

Glossary

   PCa
          prostate cancer

   PRAD
          prostate adenocarcinoma

   mCRPC
          metastatic castration-resistant prostate cancer

   scRNA-seq
          single-cell RNA sequencing

   TME
          tumor microenvironment

   NK
          natural killer

   TCGA
          The Cancer Genome Atlas

   FPKM
          ragments per kilobase of exon model per million mapped fragments

   RNA-Seq
          RNA sequencing

   TPM
          transcripts per million

   PFI
          progression-free interval

   GEO
          Gene Expression Omnibus

   PCA
          principal component analysis

   UMAP
          uniform manifold approximation and projection

   DEGs
          differentially expressed gene

   LASSO
          least absolute shrinkage and selection operator

   KM
          Kaplan-Meier

   ROC
          receiver operating characteristic

   TIME
          tumor immune microenvironment

   RPMI
          Roswell Park Memorial Institute

   FBS
          fetal bovine serum

   qRT-PCR
          quantitative real-time PCR

   NCBI
          National Center for Biotechnology Information

   PBS
          phosphate buffered saline

   ssGSEA
          single sample gene set enrichment analysis

   TIIC
          tumor infiltrating immune cell

   KEGG
          kyoto encyclopedia of genes and genomes

   TIDE
          tumor immune dysfunction and exclusion

   AUC
          Area Under the Curve

   CAF
          cancer associated fibroblast

   FOX
          forkhead-box

   TAM
          tumor-associated macrophages

   Treg
          regulatory T cell

   TRAP
          transmembrane adaptor protein

   TGF-β
          transforming growth factor-β.

Funding Statement

   The author(s) declare that financial support was received for the
   research and/or publication of this article. This study was supported
   by National Natural Science Foundation of China (81802580).

Data availability statement

   The datasets presented in this study can be found in online
   repositories. The names of the repository/repositories and accession
   number(s) can be found in the article/[148] Supplementary Material .

Ethics statement

   Ethical approval was not required for the studies on humans in
   accordance with the local legislation and institutional requirements
   because only commercially available established cell lines were used.

Author contributions

   KZ: Data curation, Formal Analysis, Methodology, Project
   administration, Validation, Visualization, Writing – original draft,
   Writing – review & editing. HX: Conceptualization, Data curation,
   Formal Analysis, Methodology, Project administration, Resources,
   Software, Visualization, Writing – original draft, Writing – review &
   editing. FZ: Formal Analysis, Investigation, Methodology, Resources,
   Supervision, Visualization, Writing – review & editing. YH:
   Conceptualization, Data curation, Formal Analysis, Funding acquisition,
   Methodology, Project administration, Supervision, Visualization,
   Writing – original draft, Writing – review & editing.

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.

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:
   [149]https://www.frontiersin.org/articles/10.3389/fimmu.2025.1564784/fu
   ll#supplementary-material
   [150]Table1.docx^ (17.8KB, docx)

References