
==== Front
Ann Med
Ann Med
Annals of Medicine
0785-3890
1365-2060
Taylor & Francis

39221762
10.1080/07853890.2024.2398195
2398195
Version of Record
Research Article
Urology
Comprehensive scRNA-seq analysis to identify new markers of M2 macrophages for predicting the prognosis of prostate cancer
Y. Ou et al.
Ou Yitian
Xia Chengxing
Ye Chunwei
Liu Mingming
Jiang Haiyang
Zhu Yong
https://orcid.org/0000-0003-3426-0305
Yang Delin
Urology Department, Kunming Medical University Second Affiliated Hospital, Kunming, Yunnan, China
Supplemental data for this article is available online at https://doi.org/10.1080/07853890.2024.2398195

CONTACT Delin Yang ydelin@163.com Urology Department, Kunming Medical University Second Affiliated Hospital, Dianmian Avenue No. 374, Kunming, Yunnan, 650101, China
2 9 2024
2024
2 9 2024
56 1 239819525 11 2023
6 5 2024
20 5 2024
KnowledgeWorks Global Ltd.30 8 2024
published online in a building issue30 8 2024
© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group
2024
The Author(s)
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

Abstract

Background

Prostate cancer (PCa) has become the highest incidence of malignant tumor among men in the world. Tumor microenvironment (TME) is necessary for tumor growth. M2 macrophages play an important role in many solid tumors. This research aimed at the role of M2 macrophages’ prognosis value in PCa.

Methods

Single-cell RNA-seq (scRNA-seq) data and mRNA expression data were obtained from the Gene Expression Omnibus database (GEO) and The Cancer Genome Atlas (TCGA). Quality control, normalization, reduction, clustering, and cell annotation of scRNA-seq data were preformed using the Seruat package. The sub-populations of the tumor-associated macrophages (TAMs) were analysis and the marker genes of M2 macrophage were selected. Differentially expressed genes (DEGs) in PCa were identified using limma and the immune infiltration was detected using CIBERSORTx. Then, a weighted correlation network analysis (WGCNA) was constructed to identify the M2 macrophage-related modules and genes. Integration of the marker genes of M2 macrophage from scRNA-seq data analysis and hub genes from WGCNA to select the prognostic gene signature based on Univariate and LASSO regression analysis. The risk score was calculated, and the DEGs, biological function, immune characteristics related to risk score were explored. And a predictive nomogram was constructed. CCK8, Transwell, and wound healing were used to verify cell phenotype changes after co-cultured.

Results

A total of 2431 marker genes of M2 macrophage and 650 hub M2 macrophage-related genes were selected based on scRNA-seq data and WGCNA. Then, 113 M2 macrophage-related genes were obtained by overlapping the scRNA-seq data and WGCNA results. Nine M2 macrophage-related genes (SMOC2, PLPP1, HES1, STMN1, GPR160, ABCG1, MAZ, MYC, and EPCAM) were screened as prognostic gene signatures. M2 risk score was calculated, the DEGs, Immune score, stromal score, ESTIMATE score, tumor purity, and immune cell infiltration, immune checkpoint expression, and responses of immunotherapy and chemotherapy were identified. And a predictive nomogram was constructed. CCK8, Transwell invasion, and wound healing further verified that M2 macrophages promoted the proliferation, invasion, and migration of PCa (p < 0.05).

Conclusions

We uncovered that M2 macrophages and relevant genes played key roles in promoting the occurrence, development, and metastases of PCa and played as convincing predictors in PCa.

Keywords

Prostate cancer
tumor microenvironment
tumor-associated macrophages
M2 macrophages
bioinformatics
nomogram
immune checkpoints
co-cultured
Natural Science Foundation of China 81860453 Yunnan Provincial Major Science and Technology Projects 2018ZF009 Top Physician Project of Yunnan 2019 Yunnan Provincial Department of Education 10.13039/501100007846 2022Y195 Yunnan Provincial Science and Technology 202301AT070321 202201AY070001-119 The Natural Science Foundation of China (No. 81860453), the Yunnan Provincial Major Science and Technology Projects (No. 2018ZF009), and the Top Physician Project of Yunnan 2019 has funded this research and the funding receiver is Delin Yang, who conceived this study. The Scientific Research Foundation Project of Yunnan Provincial Department of Education (2022Y195) funded this research and the funding receiver is Yitian Ou, who conducted and wrote this manuscript. The Project of Yunnan Provincial Science and Technology (202301AT070321, 202201AY070001-119) has funded this research and the funding receiver is Chengxing Xia, who contributed to the drafting and data analysis.
==== Body
pmcIntroduction

Prostate cancer is a malignant tumor that occurs on the epithelium of the prostate of which the most common pathological type is prostate adenocarcinoma (PRAD). According to the GLOBOCAN statistics of the World Health Organization in 2021, the incidence of PCa ranks first among male malignant tumors in the world, and its fatality rate is second [1]. The clinical course of PCa patients has great heterogeneity. Localized PCa patients can survive a long time without tumor progression, while invasive and metastatic PCa progress rapidly and lack effective treatment [2,3]. Androgen deprivation therapy (ADT) is an important treatment of PCa. However, 10–20% of the patients will progress to castration-resistant prostate cancer (CRPC) in 5 years and 15–33% of the patients will have tumor metastasis [4]. The above conditions do great harm to the prognosis of patients, so more effective treatment and monitoring methods are urgently needed.

Including tumor and non-tumor components, TME is indispensable for tumor growth. Non-tumor components include infected immune cells, blood vessels, fibrous cells, and other components. The infected immune cells play an important role in the occurrence and progression of tumors. Many scholars have shown that cytotoxic T lymphocytes, mast cells, and natural killer cells (NK cells) greatly affect the progression, metastasis, and drug resistance of PCa [5–8].

Macrophages are important myeloid cells involved in tumor immunity. Macrophages with different characteristics play different or even opposite roles [9]. The two common types of macrophages are M1 and M2. M1 is known as anti-angiogenesis, promotes chronic inflammation, and inhibits tumor growth. M2 promotes blood vessel growth, organ fibrosis, tumor growth, etc. In TME, M1 and M2 macrophages play opposite roles. The M1/M2 ratio is an important index to evaluate the prognosis of tumors. The macrophages recruited in the TME are called tumor-associated macrophages (TAMs), which mostly show the characteristics of M2 macrophages and make it easier for tumor cells to escape and metastasize.

Single-cell sequencing is a new technique for sequencing transcriptome at the single-cell level. The analysis based on this technique can study the gene expression in a single cell and solve the problems that are difficult to achieve due to the limitations of the current technology, which makes it possible to analyze different components of TME and their interactions. In recent years, the value of single cell sequencing in the study of tumor and non-tumor diseases has been continuously explored [10,11]. The purpose of this study was to use a single cell sequencing technique to explore the prognostic value of M2 macrophages in PCa.

Immune checkpoint refers to the internal regulation mechanism of the immune system, which can maintain its own tolerance and help to avoid collateral damage during the physiological immune response. Current immune checkpoint inhibitors (ICIs) have limited anticancer activity in treating PCa [12]. New checkpoints are constantly being proposed to better monitor and put forward new treatment therapies. One of our objectives was to screen M2-related checkpoints from common immune checkpoints and we got 34 of them. They can be used as follow-up immune targets and ICIs studies.

In this study, we described the immune landscape of PCa from different dimensions. By screening the M2-related genes and scoring each sample, we found that M2 macrophages promote the occurrence and progress of PCa. The M2-related genes with the highest degree of survival were screened, and they were further analyzed by the least absolute shrinkage and selection operator (LASSO) model. We verified that the M2 risk score was an independent predictor for PCa. Besides providing checkpoint targets for the follow-up study, our study provides a nomogram for predicting the PCa prognosis. Then, we analyzed the sensitivity of immunotherapy and chemotherapy between M2 high/low risk groups. To further verify the tumor-promoting effect of M2 TAM, we used PCa cell lines co-cultured with M2 macrophages. In CCK8, Transwell, and wound healing, co-cultured PCa cells showed a higher level of proliferation, invasion, and migration.

Materials and methods

Data acquisition and processing

The data analysis process of this study is shown in Supplementary Figure 1.

The bulk RNA-seq data (TCGA-PARD cohort) and corresponding clinical profiles of PRAD patients were obtained from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/) in Tumor Cancer Genomic Atlas (TCGA), which including 496 tumors and 523 normal samples. Besides, gene profiles and clinical data from the GSE54691 (104 samples) and GSE116918 (248 samples) datasets were collected from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/gds) and used for external validation. Moreover, a single cell transcriptome of 6 PRAD samples (P1-P6) was selected from the GSE137829 dataset in GEO.

scRNA-seq data processing and analysis

In the present study, the Seurat R package (version 4.1.0) was used for scRNA data pre-processing. For the GSE137829 dataset, the data of P2 could not be clustered with others and excluded, the P1, P3, P4, P5, and P6 were combined for subsequent analyses. Cluster analysis was carried out by Harmony and CCA. After integrating, the data were processed according to two allegation parameters: (1) To ensure the gene expression abundance and richness of each cell, the lowest threshold was 100 and the maximum threshold was unlimited. (2) The proportion of mitochondrial genes detected by each cell. The over-expression of mitochondrial genes indicates that the cells are in a ‘dying’ state, and will lead to biased results. The cells with the threshold value of nFeature <7500, the number of genes <5000, and the proportion of mitochondrial genes <15% were enrolled in further analysis. Then, using the ‘NormalizeData’ function, single-cell gene expression data were normalized. The ‘FindVariableFeatures’ function was used for calculating highly variable genes (HVGs). Afterward, Principal component analysis (PCA) was performed, and the top 50 PCs were used for cluster classification. The ‘FindNeighbors’ function was used for clustering based on the initial 50 PCs, further visualized through the Uniform Manifold Approximation and Projection (UMAP).

Cell type identification

The ‘FindAllmarkers’ function was used to identify the DEGs in each cell type. Cell types were recognized according to the marker genes from previous study [13]. Moreover, the subpopulations of monocytic cells also were identified based on this previous study [13]. The cellmarker2.0 was used to label the cell subsets of TAM, which are the subpopulations of monocytic cells. The markers of M1 and M2 were extracted for follow-up analysis. The sample score of marker genes of M1 and M2 was calculated by the PercentageFeatureSet function of the R studio Seurat package. Then the involved pathway analysis was done by GO and Reactome analysis.

Screening of DEGs

The DEGs between cancer and benign tissues were screened by limma package (V3.54.1) with the threshold of |log2 (fold change, FC)| > 1 and adj. p < 0.05.

Immune infiltration analysis

CIBERSORTx (https://cibersortx.stanford.edu/) was used to calculate the immune infiltration. The differences in the fractions of immune cells were detected using the Kruskal-Wallis test.

WGCNA analysis and screening of M2 key genes

Weighted gene co-expression network analysis (WGCNA) was used to screen the gene modules and key genes related to M2 TAM. The expression matrix of cancer tissue samples in TCGA-PRAD was extracted. According to the median score of M2 TAM immune infiltration, the cancer samples were divided into high and low-score groups. The limma package in R Studio was used to differentiate genes between the high/low score group of M2 TAM immune infiltration. The screening conditions were p < 0.05 and LogFC > 0. This method was used to find the positive correlation module which was significantly related to M2 TAM. In WGCNA, the first 25% genes were retained, and the soft threshold was screened by power threshold, according to the standard of hybrid dynamic tree cutting. Each gene module contains at least 100 genes. We drew the heat map of the correlation between the module and the shape and selected the most significant module through the correlation coefficient and significant p-value of the module. Module membership (MM) and Gene significance (GS) were screened. The genes in these modules were selected for follow-up analysis. The genes selected by single cell M2 macrophage subset markers and WGCNA were intersected by a Venn map and used as M2 macrophage-related genes for follow-up research.

Calculation of the M2 risk score and the effectiveness analysis

Based on the screened M2 macrophage-related hub genes, the markers obtained by VENN map and scRNA-seq analysis were intersected. The univariate Cox regression model was analyzed by using the coxph function of R studio. p < 0.1 was used to obtain the prognostic chemokine-related genes. Then, LASSO regression analysis was performed to reduce overfitting and to construct a prognostic gene signature using the glmnet package in R. Afterward, the M2-risk score was calculated according to the gene expression level and its regression coefficient.

The formula is as follows: Score=∑i=0nβi×xi

βi: weight coefficient of each gene and xi: expression of each gene.

Based on the signature score, cox regression analysis was carried out, and the over-fitted redundant factors were removed by LASSO regression analysis. Genes related to prognosis were obtained.

Then, all patients, in the training set TCGA-PRAD and external validation set GSE54691, and GSE116918, were distributed into M2-high/low-risk group. The Kaplan–Meier (KM) curve with the log-rank test was conducted using the survminer R package. The time-dependent receiver operating characteristic (ROC) curves were drawn using the timeROC package in R, and the area under the curve (AUC) value was used to estimate the predictive ability of the risk model.

Construction of M2 prognosis-related genes nomogram

Clinical information, such as Age, pathological T stage, and other information as well as M2 risk score were incorporated into univariate and multivariate Cox models to identify the independent predictors. Then, a nomogram was constructed based on the independent predictors and took the overall survival (OS) as the observation index. According to the contribution of each influencing factor to the survival risk in the model, the accuracy was tested by a calibration curve. We used the decision curve analysis (DCA) to explore the maximum net income threshold of different variables to the survival outcome, to test the validity and practicability of the line chart model, and then draw the ROC curves of 2, 3, and 4 years to verify its reliability.

HALLMARK pathway scoring and enrichment of M2 risk groups

The DEGs between M2 high and low-risk score groups were identified using the limma R package with the cutoff values of adj. p < 0.05. Then, the GSVA package was used to identify the significantly enriched pathways between M2 high and low-risk score groups based on the HALLMARK pathways.

Immune characteristics in M2 risk groups

ESTIMATE R package was used to score immune infiltration. ESTIMATE was based on the ssGSEA algorithm to determine the Immune score, ESTIMATE score, Strome score, and tumor purity of each ­sample. The significant differences between M2 high and low-risk score groups were detected using the Kruskal-Wallis test. Moreover, COBERSORT methods in the IOBR package in R were also used to analyze the immune cell infiltrating between M2 high and low-risk score groups, and the different fractions of immune cells between M2 high and low-risk score groups were detected using the Kruskal–Wallis test.

Immune checkpoints, immunotherapeutic, and chemotherapeutic responses in M2 risk groups

Tumor Immune Dysfunction and Exclusion (TIDE) score was calculated for each sample in high and low-risk groups on the TIDE website (http://tide.dfci.harvard.edu/). The immunophenoscore (IPS) also was obtained without bias by detecting the expression of four categories of immunogenicity-determining genes, including effector cells, immunosuppressive cells, MHC molecules, and immunomodulators. TIDE was negatively associated with the immunotherapeutic response, but IPS was positively associated with the immunotherapeutic response. Moreover, the expression of 34 immune checkpoints between M2 high and low-risk groups was detected. The responses of patients to immunotherapy between high and low-risk groups were analyzed using the Submap algorithm. oncoPredict R package was employed to screen the antitumor drugs whose sensitivity was associated with prognostic genes. The half-maximal inhibitory concentration (IC50) of different drugs was calculated based on the Genomics of Drug Sensitivity in Cancer database (GDSC2, https://www.cancerrxgene.org/), and the differences between high and low-risk groups were detected using the Kruskal–Wallis test.

Cell line co-culture

Prostate Cancer cell lines LnCAP (ZQ0039), Du145 (ZQ0037), and PC3 (ZQ0040) were provided by Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd. THP-1 (CL-0233) cell line was provided by Wuhan Pricella Biotechnology Co., Ltd. We used THP-1 cells to induce M0 macrophages by adding PMA (100 ng/ml, 48 h) and continued to add IL-4 (20 ng/ml, 48 h) and IL-13 (20 ng/ml, 48 h) to induce M2 macrophages. Cells were washed with PBS and collected with a scraper. M2 macrophages were cultured in upper Transwell chambers (0.4 μm). Then PCa cell lines were put it into lower chambers, respectively. Cells were used in the subsequent experiments after incubation for 48 h.

Cell phenotypic experiment

CCK-8 assay

Cell supernatant was removed from the culture plate, and 100ul basic medium was added to each well. After setting up the blank control, we added 100ul basic medium and 10ul cck-8 solution to each hole. After 2 h of incubation in the dark, the OD value of each hole was detected at 450 nm by an enzyme labeling instrument. Cell viability% = (experimental OD-blank hole OD)/(control hole OD-blank hole OD)*100%.

Wound healing assay

The cells were digested and collected by trypsin, and the cell concentration was adjusted to 2.5 × 104/ml. The cells were inoculated into the Ibidi scratch plug-in. The left and right holes of the plug-in were inoculated with 100ul cell suspension. The plug-in were removed when the overnight culture reached the fusion degree of 95%. After washing, the photos were taken at 0, 12, and 24 h.

Transwell invasion

Matrigel matrix glue was melted at 4 °C. The gun heads, Transwell chambers, and 24-hole plates should be precooled at 4 °C. The whole process of glue laying should be operated on ice. Then 50ul Matrigel was added to the Transwell film (8 μm). The Matrigel should be evenly spread on the Transwell film and dried for 30 min at 37 °C. Cell concentration was adjusted to 1 × 105/ml with the culture supernatant, and the 100ul cell suspension was inoculated into the Transwell-Matrigel chambers. The three holes in each group were repeated and incubated for 24 h. Transwell chambers were put into 1 ml 4% polymethyl fixative. After removing the upper chamber liquid, we dripped 200ul 4% polymethyl fixative to the upper chamber for 30 min. The fixing solution was removed and the chambers were put into 500ul crystal violet dye solution for 15 min. After 3-times washing by PBS, the results were photoed under the microscope.

Results

scRNA-seq reading, quality control, and hypervariable gene recognition

After preprocessing the scRNA-seq of GSE137829, we found that the data of P2 and other samples could not be clustered, so it was excluded. After combining P1, P3, P4, P5, and P6, a total of 25,313 cells and 30,074 genes (Supplementary Figures 2A–C) were obtained. After quality control, we gained 23,866 cells and 30,074 genes (Supplementary Figures 2D–F). After that, the relationship between gene proportion and sequence number of single-cell mitochondria, and the relationship between sequence number and gene number were drawn (Supplementary Figures 2G,H).

Single-cell principal component analysis, UMAP dimension, and cell population annotation

After normalization, principal component analysis was carried out. Samples and genes (Supplementary Figures 3A,B) representing the characteristics of the data set were selected from the transcriptome high-throughput sequencing by dimension, and the 1st to 50th principal components were selected for follow-up analysis.

Based on the first 50 main components, the adjacent cells were determined by Find Neighbors function. The cells were grouped by FindClusters, and 34 cell groups were obtained. Then cluster analysis (Supplementary Figures 3C,D) was carried out by UMAP.

Bubble diagram shows the expression of cell marker genes in different cell groups. Using the above methods to annotate the cell population. The markers used are basically the same as the original data set (Supplementary Figure 4A). We got a result of seven cell types (Supplementary Figure 4B): Basal/intermediate (KRT19, KRT18, and KRT8), Endothelial (ENG, VWF, and PECAM1), Fibroblast (ACTA2), Luminal (KRT18, KRT8, and AR), Mast Cell (TPSB2, TPSAB1, and MS4A2), Monolytic (FCGR3A, CSF1R, CD14, LYZ, CD163, and CD68), T Cell (CD7, CD3G, CD3D, CD3E, and CD2). A small cell group was classified as unknown because the expression characteristics were not obvious. We selected the most significantly expressed markers, and drew their expression to determine the different cell groups (Supplementary Figures 4C,D).

Monolytic subgroup analysis

After Extracting the monolytic cells, we selected the three main myeloid immune cells: monocytes (FCN1, RETN, and S100A8), dendritic cells (HLA-DOB, CLNK, and IDO1), and TAM cells (C1Q1A, C1Q1C, and C1Q1B) for annotation (Supplementary Figures 5A,B). Then we drew the marker expression highlight map and violin map (Supplementary Figures 5C,D).

TAM subgroup analysis

We further explored the differences between M1 (S100A8, FAM26F, and FOS) and M2 (TGFB1, MS4A4A, FCGR3A, FOLR2, and CD163) in TAM subsets. Using markers in the cellmarker2.0 database, TAM was divided into M1 and M2 macrophage subsets, and the results were shown in Supplementary Figures 6A–F. Cellular markers were extracted for follow-up analysis (2431 genes).

Go and Reactome pathway enrichment

The marker genes of M1 and M2 macrophages were analyzed by GO and Reactome pathway enrichment. In GO Biological Process (GOBP) analysis, M1 marker genes were enriched to myeloid leukocyte activation, regulation of immune effector process, and other pathways (Supplementary Figure 7A). M2 marker genes were enriched to positive regulation of cytokine production, leukocyte mediated immunity pathway (Supplementary Figure 7E). In GO Cellular Component (GOCC) analysis, M1 and M2 marker genes were enriched to secretory granule membrane, tertiary granule, and other pathways (Supplementary Figures 7B,F). In GO Molecular Function (GOMF) analysis, M1 and M2 marker genes were enriched to immune receptor activity, MHC protein complex binding, and other pathways (Supplementary Figures 7C,G). In the Reactome pathway enrichment analysis, M1 marker genes were enriched to Neutrophil degranulation and interleukin-10 (IL-10) signaling pathway (Supplementary Figure 7D). M2 marker genes were enriched to Neutrophil degranulation, Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell pathway (Supplementary Figure 7H).

TCGA DEGs screening

Four hundred and ninety-six tumor samples and 52 normal samples of TCGA-PRAD were analyzed to screen DEGs. The condition was that |LogFC| > 1 and q < 0.05. 524 DEGS were screened, of which 159 were up-regulated and 365 were down-regulated (Supplementary Figure 8).

CIBERSORTx of immune cells infiltration in PCa TME

The DEGs were used to calculate the immune infiltration (Supplementary Figure 9A) using CIBERSORTx. The difference in immune infiltration between cancer and benign tissues (Supplementary Figure 9B) was shown by a box map. The results showed a significant difference in immune infiltration between cancer and benign tissues among all immune cells except T.cells.gamma.delta.

Screening of M2 related modules and key genes by WGCNA

The DEGs were screened between high and low M2 immune infiltration groups, and a total of 4519 up-regulated genes were obtained in high groups (Supplementary Figure 10A). After WGCNA analysis (Supplementary Figures 10B–D), three modules with the highest correlation with M2 macrophages (Brown, turquoise, green. Supplementary Figures 10E–H) were selected, and 650 hub genes from these three modules were extracted for follow-up research.

Construction and efficacy evaluation of prognostic model based on M2 marker genes

Two thousand four hundred and thirty-one M2 macrophage markers obtained from scRNA-seq data analyses were intersected with 650 M2 hub genes obtained from bulk RNA-seq data analyses. A total of 113 intersected genes were obtained (Supplementary Figure 11A). A total of 12 M2 related genes were obtained by COX regression analysis. After removing redundant factors by LASSO, nine prognostic related genes were selected: SMOC2, PLPP1, HES1, STMN1, GPR160, ABCG1, MAZ, MYC, and EPCAM (Supplementary Figures 11B–D).

Effectiveness analysis of M2 risk scoring model of PRAD

In the training set TCGA-PRAD, the prognosis of M2 high-risk group was worse than low risk-group (Figure 1A, p = 0.0462). The ROC curves of 2, 3, and 4 years were drawn, and the corresponding AUC values were 0.648, 0.806, and 0.860, respectively (Figure 1B). The risk score (Figure 1C), survival time distribution (Figure 1D), and the expression of prognostic factors (Figure 1E) in the high/low score group were shown by Figure 1, indicating that the model has a good prognostic value for the training set.

Figure 1. M2 risk scoring model’s prognostic value of training set TCGA-PRAD. (A) KM curve of TCGA training set. The survival prognosis of high-score group was significantly lower than that of low-score group (p = 0.0462). (B) ROC curves of 2, 3, and 4 years. (C) Risk score results. (D) Survival time distribution. (E) Expression of prognosis-related genes.

In the verification set GSE54691, the prognosis of M2 high-risk group was worse than that of low-risk group (Figure 2A, p = 0.0063). The ROC curves 2, 3, and 4 years were drawn, and the corresponding AUC values were 0.653, 0.680, and 0.684, respectively (Figure 2B). The risk score (Figure 2C), survival time distribution (Figure 2D), and the expression of prognostic factors (Figure 2E) in the high/low score group are shown in Figure 2, indicating that the model has a certain prognostic value for this verification set.

Figure 2. M2 risk scoring model’s prognostic value of verification set GSE54691. (A) KM curve of GSE54691 verification set. The survival prognosis of high-score group was significantly lower than that of low-score group (p = 0.0063). (B) ROC curves of 2, 3, and 4 years. (C) Risk score results. (D) Survival time distribution. (E) Expression of prognosis-related genes.

In the verification set GSE116918, the prognosis of M2 high-risk group was worse than that of low-risk group (Figure 3A, p = 0.02). The ROC curves of 2, 3, and 4 years were drawn, and the corresponding AUC values were 0.942, 0.809, and 0.658, respectively (Figure 3B). The risk score (Figure 3C), survival time distribution (Figure 3D), and the expression of prognostic factors (Figure 3E) in the high/low score group are shown in Figure 3, indicating that the model has a good prognostic value for this verification set.

Figure 3. M2 risk scoring model’s prognostic value of verification set GSE116918. (A) KM curve of GSE116918 verification set. The survival prognosis of high-score group was significantly lower than that of low-score group (p = 0.02). (B) ROC curves of 2, 3, and 4 years. (C) Risk score results. (D) Survival time distribution. (E) Expression of prognosis-related genes.

The univariate and multivariate analysis of age, pathological stage, and M2 risk score, the results were as shown in Table 1. Age is not an independent predictor. T staging only showed a significant difference in univariate analysis (p = 0.0096). M2 risk score was an independent predictor of PCa prognosis in both univariate (p = 9.10E-08) and multivariate analysis (p = 0.0098).

Table 1. Univariate and multivariate analysis of age, T stage, M2 risk score.

 	Univariate cox	Multivariate cox	
p-Value	HR	p-Value	HR	
Age ≥ 60	0.78	1.2 (0.34–4.3)	0.93	0.92 (0.15–5.7)	
T stage	0.0096	6.6 (1.6–28)	0.24	3 (0.48–19)	
M2 risk score	9.10E-08	4.8 (2.7–8.5)	0.0098	4.1 (1.4–12)	

Age, T stage, and M2 risk score were used to construct a nomogram (Figure 4A) for predicting OS PRAD patients. The calibration curve showed that the predicted OS was in good consistency with the actual survival time, indicating that the nomogram had a high accuracy (Figure 4B). The reliability was verified by ROC, and the AUC values of 2, 3, and 4 years were 0.703, 0.810, and 0.867, respectively (Figure 4C). DCA results showed that among all the risk factors, the model had the highest prognostic net income value, indicating that the model had good validity and practicability in predicting OS (Figure 4D).

Figure 4. PRAD OS predicting nomogram of M2 risk score. (A) Predictive OS nomogram of PRAD patients. (B) The consistency curve of actual OS and predicted OS in 2, 3, and 4 years. (C) ROC curves of 2, 3, and 4 years and corresponding AUC. (D) The DCA results showed the net income values of different prognostic factors, and complex represents the complex of the factors in the model.

DEGs screening and pathway analysis of M2 high and low-risk groups

After screening the DEGs between M2 high/low risk group, 8629 DEGs were obtained, including 4826 up-regulated genes and 3803 down-regulated genes (Supplementary Figure 12A). The HALLMARK pathway scores of 8629 DEGs in each sample were calculated, and the pathway heat map (Supplementary Figure 12B) was drawn. The results showed that there was a significant statistical difference between M2 high/low-risk group in 41 pathways.

Immune infiltration analysis between M2 high/low-risk groups

The results showed significant differences in ESTIMATE Score, Immune Score, Stromal Score, Tumor Purity between M2 high/low-risk groups (p < 0.01, Figures 5A,B). CIBERSORT was used to calculate the infiltration degree between M2 high/low-risk groups. The results showed that there were significant differences in naive B cells, M2 macrophages, and CD4+ memory T cells (p < 0.05. Figures 5C,D).

Figure 5. Immune cells infiltration analysis. (A) Heat map (pheatmap, ver. 1.0.12) and (B) Box diagram of ESTIMATE score, immune score, stromal score, and tumor purity between M2 high/low-risk groups. (C) Heat map (pheatmap, ver. 1.0.12) and (D) Box diagram of immune infiltration between M2 high/low-risk groups.

Immune checkpoints, immunotherapeutic and chemotherapeutic responses in M2 high/low-risk groups

We further screened 73 immune checkpoints and found significant differences of 34 immune checkpoints between high/low M2 groups (Figures 6A,B). The immunotherapeutic responses results (Figure 6C) showed that the expression patterns of patients with high M2 score were more similar to those with CTLA4-nOR inhibitor response (p = 0.07). The expression pattern of M2 low score group was more similar to that of PD1-R inhibitor response group (p = 0.15). Figure 6D shows the differential chemotherapeutic drugs between M2 high/low group: Vinblastine, Cytarabine, Gefitinib, Vorinostat, Doramapimod, Wee1.Inhibitor, Camptothecin, Cisplatin, Docetaxel, Navitoclax, Nilotinib, Olaparib.

Figure 6. Immune checkpoints, immunotherapeutic and chemotherapeutic responses analysis. (A) Box diagram and (B) Heat map (pheatmap, ver. 1.0.12) of differentially expressed immune check points between M2 high/low risk groups. (C) Submap analysis showed immunotherapy sensitivity in high/low M2 groups. (D) Differential chemotherapeutic drugs between M2 high/low group.

Promoting effect of M2 TAM on prostate cancer

In order to further verify the promoting effect of M2 TAM on PCa. We verified the proliferation, migration, and invasion ability of PCa cells before and after co-culture by CCK8, wound healing, and Transwell, respectively. The results showed that the ability of proliferation and invasion in all M2 co-cultured PCa cell lines was higher than that in normally cultured cell lines (p < 0.01, Figures 7A–C,H–K). The migration ability of DU145 and PC3 cell lines was enhanced after co-culture (p < 0.05, Figures 7D–G), but there was no significant change in LnCAP cell line.

Figure 7. M2 TAMs Promote PCa cell proliferation, migration, and invasion in vitro. PCa cell lines (LnCAP, DU145, PC3) were non-contact co-cultured with M2 macrophages. (A–C) Cell proliferation was quantified by CCK-8 assay. (D–G) Cell migration was detected by wound healing assay. (H–K) Cell invasion was measured by Transwell assay. *p < 0.05, **p < 0.01, ***p < 0.001. Each experiment was repeated three times.

Discussion

PCa has become a major threat to the health of male. The progression and prognosis of PCa vary greatly due to the spatial and clonal heterogeneity [14]. The development of single-cell technology provides the possibility for more accurate study of PCa transcriptome information. We depicted the immune landscape of PCa. Focusing on the analysis of the function of M2 macrophages in PCa, we screened out M2-related DEGs and the mainly involved pathways, such as Wnt pathway.

TME contains a variety of tumor and non-tumor components. These components interact with each other and form a dynamic system [15]. The interaction of various components in PCa TME has also become a hotspot in urology in recent years. Kwon et al. [16] summarized the interaction of cytotoxic T cells, B cells, M2 macrophages, and other immune cells through different cytokines in PCa. Bahmad et al. [17] concluded that mutual interference between epithelial cells and cellular stroma plays an indispensable role in progression and metastasis in PCa. Using single cell sequencing and analysis, Chen et al. [18] found that monocytes, dendritic cells affect tumor progression, and TME promoted the metastasis and spread of PCa when it was formed early. In this study, we found that M2 macrophage infiltration was related to the prognosis of PCa, and higher M2 score corresponds to worse prognosis.

When TAMs are gathered to TME, they release a variety of cytokines and promote tumor formation. Studies [19] have shown that after tumor formation, TAM can promote tumor growth, progression, and metastasis by ① Secreting vascular endothelial growth factor (VEGF) and epidermal growth factor (EGF) to promote tumor microvascular growth. ② Secreting interleukin-1 (IL-1), colony-stimulating factor-1 (CSF-1), matrix metalloproteinases (MMPs) to promote tumor invasion and metastasis. ③ Secreting interleukin-10 (IL-10), prostaglandin E2 (PGE2), transforming growth factor-β (TGF-β) to promote tumor immune escape. In our study, we found that M2 marker genes were related to positive regulation of cytokine production, leukocyte mediated immunity pathway, Neutrophil degranulation, Immunoregulatory interactions. Our follow-up study will focus on these pathways to explore the interaction between M2 macrophages and tumors.

Meanwhile, the study of regulating the polarization of TAM is also the key of tumor immunity. Li et al. [20] have found that MS4A4A promoted M2 polarization of TAM by activating PI3K/AKT pathway and JAK/STAT6 pathway. Cheng et al. [21] have found that tumor necrosis factor α-induced protein 8-like 1 (TIPE1) promotes the activation of PI3K/Akt pathway in TAM by directly binding and regulating the metabolism of phosphatidylinositol 4-diphosphate (PIP2) and phosphatidylinositol 3-pyrrosine 5-trisphosphate (PIP3). In this study, we found that MYC and MAZ genes are M2 key genes related to prognosis. In our follow-up study, we screened the prognosis-related gene MS4A6A in PCa TME and found that MS4A6A was associated with MYC and MAZ gene promoter through double luciferase report experiment. We will continue to study the mechanism of MS4A6A-MYC/MAZ in TAM polarization.

In clinic, the prognosis of patients with PCa is affected by great individual differences. Some patients can survive with tumors for a long time without obvious progress, while others will enter the stage of CRPC even after proper treatment. The existing clinical prognosis detection methods still have many deficiency. With the development of genomics and transcriptome research, more subtypes of PCa have been revealed, which provides the possibility for new immune checkpoints affecting the progression and metastasis of PCa [22,23]. He et al. [12] concluded that MMR gene deficiency can enhance PRAD anti-tumor immune response, showing more tumor infiltrating lymphocytes. Some mCRPC patients will have CDK12 aberration, and the inactivation of CDK12 will lead to focal tandem repetition, which will enhance tumor immune response. Abida et al. [24] treated 11 mCRPC patients with metastatic high micro-satellite instability or mismatch repair defect with immune checkpoint inhibitors, of which 6 (54.5%) had PSA decreased by more than 50%. In our study, nine prognostic related genes were screened, and a nomogram prognostic prediction model was constructed, which showed a good ability to predict prognosis in both training and verification sets. The predicting results were consistent with clinical prognosis. We further confirmed that M2 score was an independent predictor of PCa prognosis.

Different from traditional therapies, tumor immunotherapy takes effect by activating immunity. Nowadays, there are common immune checkpoint inhibitor targets, such as PD-1/PD-L1, CTLA-4. Studies have shown that CTLA-4 inhibitor Epimazumab has a certain effect on early PCa patients, but no improvement in OS in CRPC patients [25]. Yang et al. Found that PD-1/PD-L1 inhibitors are poor-effected when used alone in lung cancer [26]. In our study, we found that the expression pattern of M2 high score group was more similar to CTLA4 inhibitor response group while the expression pattern of M2 low score group was more similar to PD1-R inhibitor response group. This may indicate that different immunotherapy regimens will be more effective for different M2 risk groups. Moreover, we found 34 immune checkpoints between high/low M2 groups, which may provide targets for follow-up research.

Among the differential chemotherapeutic drugs between high/low M2 groups, some have commonly used in PCa chemotherapy, such as Gefitinib, Vorinostat, Vinblastine, Cisplatin, Docetaxel, etc. [12]. In our follow-up study, we will further clarify whether better efficacy can be achieved by using different chemotherapeutic drugs for different M2 risk groups.

This study also has some limitations. Due to the high cost of single cell sequencing technology, it has not been widely popularized, but with the continuous improvement of single cell sequencing, the technical limitations will continue to be overcome. This study screened the M2 key gene in prostate cancer TME and showed a good prognostic value through the nomogram, and further verified the promoting effect of M2 macrophages on prostate cancer in vitro.

Conclusions

In prostate cancer, it is further verified that M2 macrophages promote the occurrence, development, and metastases of PCa. A Higher M2 score indicates worse prognosis. M2 macrophage-related differential genes are mainly involved in tumor-related pathways, such as Wnt pathway. Risk score model and nomogram have good ability to predict the survival and prognosis of PCa. M2 risk score was an independent predictor of PCa prognosis.

Ethical approval

The data used in this study are sourced from publicly available databases and do not involve the direct collection of data from human participants. According to the guidelines of the local ethics review board, research using publicly available data is exempt from requiring ethical approval. Therefore, no additional ethical approval is needed for this study.

Supplementary Material

Clean version：Supplementary Figures.docx

Authors contributions

The study was conceived by Delin Yang. Yitian Ou conducted the research and wrote the manuscript. The literature search and data collection were done by Yitian Ou, Chengxing Xia, Chunwei Ye, and Haiyang Jiang. Yitian Ou, Mingming Liu, and Yong Zhu all contributed to the manuscript’s drafting and data interpretation. The essay was co-authored by all writers, who all gave their approval to the final edition.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All datasets are open for public use. The TCGA-PRAD dataset can be downloaded at https://portal.gdc.cancer.gov/projects/TCGA-PRAD. GSE54691 dataset can be downloaded at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54691. GSE116918 dataset can be downloaded at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116918. GSE137829 dataset can be downloaded at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE137829.
==== Refs
References

1 Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71 (3 ):209–249. doi: 10.3322/caac.21660.33538338
2 Sartor O. Localized prostate cancer – then and now. N Engl J Med. 2023;388 (17 ):1617–1618. doi: 10.1056/NEJMe2300807.36912567
3 Achard V, Putora PM, Omlin A, et al. Metastatic prostate cancer: treatment options. Oncology. 2022;100 (1 ):48–59. doi: 10.1159/000519861.34781285
4 Mansinho A, Macedo D, Fernandes I, et al. Castration-resistant prostate cancer: mechanisms, targets and treatment. Adv Exp Med Biol. 2018;1096 :117–133. doi: 10.1007/978-3-319-99286-0_7.30324351
5 Elsässer-Beile U, Przytulski B, Gierschner D, et al. Comparison of the activation status of tumor infiltrating and peripheral lymphocytes of patients with adenocarcinomas and benign hyperplasia of the prostate. Prostate. 2000;45 (1 ):1–7. doi: 10.1002/1097-0045(20000915)45:1<1::AID-PROS1>3.0.CO;2-V.10960837
6 Hu S, Li L, Yeh S, et al. Infiltrating T cells promote prostate cancer metastasis via modulation of FGF11→miRNA-541→androgen receptor (AR)→MMP9 signaling. Mol Oncol. 2015;9 (1 ):44–57. doi: 10.1016/j.molonc.2014.07.013.25135278
7 Cortesi F, Delfanti G, Grilli A, et al. Bimodal CD40/Fas-dependent crosstalk between iNKT cells and tumor-associated macrophages impairs prostate cancer progression. Cell Rep. 2018;22 (11 ):3006–3020. doi: 10.1016/j.celrep.2018.02.058.29539427
8 Xie H, Li C, Dang Q, et al. Infiltrating mast cells increase prostate cancer chemotherapy and radiotherapy resistances via modulation of p38/p53/p21 and ATM signals. Oncotarget. 2016;7 (2 ):1341–1353. doi: 10.18632/oncotarget.6372.26625310
9 Chávez-Galán L, Olleros ML, Vesin D, et al. Much more than M1 and M2 macrophages, there are also CD169(+) and TCR(+) macrophages. Front Immunol. 2015;6 :263.26074923
10 Yang Q, Zhang H, Wei T, et al. Single-cell RNA sequencing reveals the heterogeneity of tumor-associated macrophage in non-small cell lung cancer and differences between sexes. Front Immunol. 2021;12 :756722. doi: 10.3389/fimmu.2021.756722.34804043
11 Krasniewski LK, Chakraborty P, Cui CY, et al. Single-cell analysis of skeletal muscle macrophages reveals age-associated functional subpopulations. Elife. 2022;11 :e77974. doi: 10.7554/eLife.77974.36259488
12 He Y, Xu W, Xiao YT, et al. Targeting signaling pathways in prostate cancer: mechanisms and clinical trials. Signal Transduct Target Ther. 2022;7 (1 ):198. doi: 10.1038/s41392-022-01042-7.35750683
13 Dong B, Miao J, Wang Y, et al. Single-cell analysis supports a luminal-neuroendocrine transdifferentiation in human prostate cancer. Commun Biol. 2020;3 (1 ):778. doi: 10.1038/s42003-020-01476-1.33328604
14 Espiritu SMG, Liu LY, Rubanova Y, et al. The evolutionary landscape of localized prostate cancers drives clinical aggression. Cell. 2018;173 (4 ):1003–1013.e15. doi: 10.1016/j.cell.2018.03.029.29681457
15 Spano D, Zollo M. Tumor microenvironment: a main actor in the metastasis process. Clin Exp Metastasis. 2012;29 (4 ):381–395. doi: 10.1007/s10585-012-9457-5.22322279
16 Kwon JTW, Bryant RJ, Parkes EE. The tumor microenvironment and immune responses in prostate cancer patients. Endocr Relat Cancer. 2021;28 (8 ):T95–T107. doi: 10.1530/ERC-21-0149.34128831
17 Bahmad HF, Jalloul M, Azar J, et al. Tumor microenvironment in prostate cancer: toward identification of novel molecular biomarkers for diagnosis, prognosis, and therapy development. Front Genet. 2021;12 :652747. doi: 10.3389/fgene.2021.652747.33841508
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 (1 ):87–98. doi: 10.1038/s41556-020-00613-6.33420488
19 Mantovani A, Allavena P, Marchesi F, et al. Macrophages as tools and targets in cancer therapy. Nat Rev Drug Discov. 2022;21 (11 ):799–820. doi: 10.1038/s41573-022-00520-5.35974096
20 Li Y, Shen Z, Chai Z, et al. Targeting MS4A4A on tumour-associated macrophages restores CD8+ T-cell-mediated antitumour immunity. Gut. 2023;72 (12 ):2307–2320. doi: 10.1136/gutjnl-2022-329147.37507218
21 Cheng Y, Bai F, Ren X, et al. Phosphoinositide-binding protein TIPE1 promotes alternative activation of macrophages and tumor progression via PIP3/Akt/TGFβ axis. Cancer Res. 2022;82 (8 ):1603–1616. doi: 10.1158/0008-5472.CAN-21-0003.35135809
22 Chen S, Huang V, Xu X, et al. Widespread and functional RNA circularization in localized prostate cancer. Cell. 2019;176 (4 ):831–843.e22. doi: 10.1016/j.cell.2019.01.025.30735634
23 You S, Knudsen BS, Erho N, et al. Integrated classification of prostate cancer reveals a novel luminal subtype with poor outcome. Cancer Res. 2016;76 (17 ):4948–4958. doi: 10.1158/0008-5472.CAN-16-0902.27302169
24 Abida W, Cheng ML, Armenia J, et al. Analysis of the prevalence of microsatellite instability in prostate cancer and response to immune checkpoint blockade. JAMA Oncol. 2019;5 (4 ):471–478. doi: 10.1001/jamaoncol.2018.5801.30589920
25 Beer TM, Kwon ED, Drake CG, et al. Randomized, double-blind, phase III trial of ipilimumab versus placebo in asymptomatic or minimally symptomatic patients with metastatic chemotherapy-naive castration-resistant prostate cancer. J Clin Oncol. 2017;35 (1 ):40–47. doi: 10.1200/JCO.2016.69.1584.28034081
26 Yang JC, Shepherd FA, Kim DW, et al. Osimertinib plus durvalumab versus osimertinib monotherapy in EGFR T790M-positive NSCLC following previous EGFR TKI therapy: CAURAL brief report. J Thorac Oncol. 2019;14 (5 ):933–939. doi: 10.1016/j.jtho.2019.02.001.30763730
