INTRODUCTION
Metastatic hormone-sensitive prostate cancer (mHSPC) is primarily responsive to androgen deprivation therapy, which has been the mainstay of treatment over past decades. Many treatment options have recently become available for patients with mHSPC, such as early docetaxel (DCT) and androgen receptor (AR)-targeting therapy, including abiraterone acetate, enzalutamide, and apalutamide; however, optimal biomarkers to guide the selection of the best treatment for such patients are lacking [1,2].
Advances in high-throughput sequencing technologies have facilitated the study of multilayered genomic, epigenomic, and transcriptomic cancer aberrations [3-6]. Furthermore, several genome-wide studies aimed at unveiling the underlying molecular mechanisms of prostate cancer (PCa) have been conducted [3-7]. However, although many studies have been performed in patients with PCa or metastatic castration-resistant PCa, the genomic and transcriptomic profiles of mHSPC and its correlation with clinical outcomes remains largely unknown [8-10].
In this study, we conducted a detailed characterization and transcriptomic profiling of mHSPC at the single-cell level. Our analysis identified luminal and mononuclear phagocyte (MNP) subpopulations, which were associated with an aggressive phenotype. Additionally, 2 distinct subtypes of mHSPC were identified using whole transcriptome sequencing (WTS) data, including a set of 87 network genes derived from the analysis of single-cell RNA-sequencing (scRNA-seq) data. We finally developed a 14-gene signature with prognostic potential for predicting patient outcomes and therapeutic responses.
MATERIALS AND METHODS
1. Tissue Collection
For WTS, we retrospectively reviewed the medical records of patients diagnosed with de novo mHSPC between July 2018 and April 2021. A total of 88 patients who received either early DCT or abiraterone acetate and prednisone (AAP) were enrolled. Patients were excluded if no tissue was available (n=26) and if RNA integrity failed (n=2) for WTS (Supplementary Fig. 1). The patients were followed for a total duration of 64.5 months, with a median follow-up period of 37.4 months for AAP and 49.3 months for DCT. scRNA-seq data were prospectively collected from 7 patients at our institution.
2. scRNA-Sequencing and Data Analysis
Twelve primary tissues from 7 patients with mHSPC (n=2), localized PCa (n=7), and matched normal tissues (n=3) were used for 10X Chromium Chip Single-cell library generation and preparation of a single-cell suspension using a 10X Chromium 3’ solution (V2 kit) according to the manufacturer’s instructions, with the aim of capturing 5,000-10,000 cells/channel. Sequencing was performed using an Illumina HiSeq 4000 platform at the Samsung Genomic Institute. Raw sequencing data were processed using Cell Ranger software (10X Genomics, ver. 3.0.2) and aligned to the human reference genome (GRCh38). Doublets and low-quality cells were excluded from subsequent analyses (<500 or >8,000 genes and >15% expression of mitochondrial genes). For the remaining 19,771 cells, gene expression data were normalized and scaled using the NormalizeData and ScaleData functions of the Seurat package, respectively. The batch effects were removed using the RunHarmony function of the Seurat package. To identify cluster marker genes and differentially expressed genes (DEGs, avg_log2FC>0.5 and p_val_adj<0.05) for scRNA-seq, the FindAllMarkers and FindMarkers functions of the Seurat package (ver. 4.1.1) were used with default parameters.
3. RNA-Sequencing Preprocess
RNA sequencing was performed using 52 metastatic mHSPCs and 8 adjacent nontumor tissues. Specifically, total RNA was extracted using the RNeasy Mini Kit (QIAGEN), and its integrity was verified using a 2100 Bioanalyzer (Agilent). Sequencing libraries were constructed using the QuantSeq 3’ Library Prep Kit (Lexogen Inc.) according to the manufacturer’s instructions, and sequenced on an Illumina HiSeq 2000 system (Illumina). The reads were mapped to the GRCh38 human reference genome with STAR (ver. 2.5.2b) using the default parameter values. The number of reads mapped to each gene was calculated using RNAseq with RSEM (ver. 1.3.0). Data processing and analysis were performed using the R/Bioconductor libraries. The transcriptome data were preprocessed by filtering out genes with zero values across samples, and data were normalized by subtracting the average expression values of adjacent nontumor tissues and centering the expression values of each sample and gene. Data processing and analysis were performed using R/Bioconductor libraries.
4. Histological Examination and Immunohistochemistry
Board-certified pathologists at our hospital verified all histological diagnoses and sample adequacy. For immunohistochemical analyses, 4-μm-thick sections were obtained from formalin-fixed paraffin-embedded tissue blocks. Immunohistochemical analysis was performed using primary antibodies against NKX3.1 (Clone24B72D11.1, Cat# 556604, BD Biosciences, 1:25), ARs (Clone: EPR20766, Cat# ab213725, Abcam, 1:250), prostate-specific antigens (PSAs) (clone polyclonal, Cat# NBP1-83966, Novus Biologicals, 1:200), and CDK1 (Clone: EPR165, Cat# ab133327, Abcam, 1:250). Positive and negative internal controls were used for all markers, with the aid of a pathologist (SH) with expertise in genitourinary tumor pathology.
5. Pseudotime Trajectory Analysis
To perform single-cell trajectory analysis for samples from patients with mHSPC (mHSPC_02 and mHSPC_03), the function learn_graph from the Monocle3 package was used [11] (https://cole-trapnell-lab.github.io/monocle3). Luminal cell subclusters with fewer than 50 cells were removed and the starting position was assigned based on the normal luminal subcluster.
6. Cell-Cell Interaction and Copy Number Variation Analysis
Cell communication analysis was performed using the CellChat package with default parameters (https://github.com/sqjin/CellChat) [12]. The copy number variations (CNVs) for each mHSPC samples were estimated using the inferCNV package (ver. 1.15.0) (https://github.com/broadinstitute/inferCNV) and the cells from the 3 normal samples were used as reference. The inferCNV analysis was performed using the parameters “denoise” and default hidden Markov model with a cutoff value of 0.1.
7. Network Analysis
The gene network was constructed using physical interactions obtained using the GeneMANIA plugin of Cytoscape (http://www.cytoscape.org/) [13]
8. Gene Set Enrichment Analysis and Single-Sample GSEA
Gene set enrichment analysis (GSEA) was performed using hallmark gene sets from the Molecular Signatures Database (MSigDB ver. 7.0) [14]. The gene sets used in this study included oncogenes, tumor suppressor genes (TSGs) [15,16], chromosomal instability genes (CIN25, CIN70) [17], human embryonic stem (hES) cells [18], CRPC51 [19] and macrophage migration inhibitory factor (MIF) pathway members including MIF, CD74, CXCR4, and CD44 [20-22]. Single-sample GSEA (ssGSEA) for WTS was computed using the “GSVA” package [23]. ssGSEA for scRNA-seq was computed using the “escape” package (ver. 1.8.0).
9. Statistical Analysis
DEGs were determined using a permutation t-test for the WTS. Gene Ontology (GO) analysis of the DEGs was performed using the DAVID web application (https://david.ncifcrf.gov). The molecular subtyping function from “genefu” package was used to identify PAM50 subtypes, including basal, luminal A, luminal B, HER2, and normal [24]. The “moonBook” package was used to compare the clinicopathological characteristics of the 2 groups. The Kaplan-Meier method was used to estimate the progression-free survival (PFS), disease-free survival (DFS), and overall survival (OS). One-way analysis of variance was performed for all 3 groups. Student t-test was performed for both groups. All statistical analyses were performed using R software.
RESULTS
1. Single-Cell Transcriptomic Profiles of Localized PCa, Normal, and mHSPC
We collected 12 fresh tissue samples to profile mHSPC at the single-cell level (Supplementary Table 1). ScRNA-seq was performed on these tissues to generate data from 19,771 cells after standard processing and quality control (Fig. 1A). Using unsupervised clustering, 20 distinct cell clusters containing both cancerous and normal cells were identified (Fig. 1B). The clusters included various cell types (basal, luminal, fibroblast, and endothelial cells), as well as multiple immune cell clusters (MNP cells, mast cells, natural killer cells [NK], T cells, and B cells) (Fig. 1C and D). The annotations for these clusters were determined by considering canonical marker gene expression (Supplementary Fig. 2A), and examining the marker genes in each cluster. Substantial variations in cell type composition among individual tissues were observed (Supplementary Fig. 2B), and the luminal cluster was found to be the dominant cell type in most tissues (Fig 1E, Supplementary 2B). It must be noted that the proportions of endothelial and T-cell clusters decreased in mHSPC tissues, whereas those of luminal and MNP clusters increased (Fig. 1E), suggesting that luminal and MNP cells may have certain roles in the mHSPC tumor microenvironment.
2. Characteristics of 10 Distinct Luminal Cell Type Subclusters
We performed reclustering to better characterize the luminal cluster, the dominant cell type, which revealed 10 distinct subclusters (Fig. 2A, left); additionally, 10 patient subclusters were identified (Fig. 2A, right). Based on the cell proportion distribution, we classified cells in subclusters L0, L3, and L5, which were predominantly derived from normal tissues, as normal luminal subclusters. The remaining 7 subclusters were characterized by a significant deviation from tumor tissues and were considered tumor luminal subclusters (Fig. 2B). Subcluster L8 exhibited exclusive enrichment in mHSPC tissues and was therefore categorized as the mHSPC luminal subcluster (Fig. 2B). Furthermore, marker genes that were specifically expressed in each luminal subcluster were identified (Fig. 2C). Interestingly, subcluster L5 expressed ACPP and NEFH, a TSG (Fig. 2C). Consequently, we examined the expression levels of oncogenes and TSG sets within each subcluster. Subclusters L3 and L5 exhibited lower oncogene expression than the other subclusters and higher TSG expression, with the exception of subclusters L8 and L9 (Fig. 2D). Furthermore, the biological functions of the DEGs within subclusters L0, L3, and L5 were investigated. Analysis revealed that subclusters L3 and L5, but not L0, were associated with functions related to homeostatic processes such as metal ion homeostasis and response to zinc and copper ions (Fig. 2E).
Hallmark gene sets were used to score the subclusters based on their expression. We found that metabolism-related gene sets—including those involved in xenobiotic and fatty acid metabolism and peroxisomes—were activated in normal and localized tissues, but not in mHSPC tissues (Supplementary Fig. 3A). Furthermore, subclusters L7, L8, and L9 showed lower expression of metabolism-related gene sets than the other subclusters, which correlated with the expression of mitochondrial gene sets (Fig 2C, Supplementary Fig. 3B). We also found that the androgen response was activated in normal and localized tissues but not in mHSPC tissues, and that this gene set was highly expressed in subclusters L3 and L5 (Fig. 2F, top). Specifically, the highest and lowest expression of the E2F target gene set was observed in the mHSPC luminal subcluster (L8) and normal luminal subclusters (L3 and L5), respectively (Fig. 2F, bottom). MYC target gene sets were found to be enriched in localized tumors compared to normal and mHSPC samples. Within the luminal subclusters (L2, L4, and L6), we observed higher expression of MYC targets compared to other subclusters (Supplementary Fig. 3C and D).
3. Presence of MNP Cycling Clusters in Patients With mHSPC Is Associated With Drug Response and Clinical Outcomes
The MNP cluster, which was enriched in mHSPC tissues, was further analyzed. Cluster 12 specifically expressed several proliferation marker genes, including UBE2C, TOP2A, CENPF, MKI67, CDC20, and CDK1. Consequently, we designated this cluster as the MNP cycling cluster based on the marker genes and GO (Fig. 3A, Supplementary Table 2). Moreover, although MNP cycling cells were identified in both patients with mHSPC (Fig. 3B), interpatient heterogeneity was observed as mHSPC_02 exhibited a higher cell proportion the MNP cycling cluster than mHSPC_03 (Fig. 3C). Expectedly, the biological functions associated with the MNP clusters using hallmark gene sets were subsequently assessed, which revealed that pathways such as the G2M checkpoint, E2F targets, and previously defined aggressive-related gene signatures (e.g., CRPC51, CIN25, CIN70, hES1, and hES2) were significantly activated in mHSPC_02 compared to mHSPC_03 (Fig. 3D, Supplementary Fig. 4A and B). However, pathways such as the inflammatory response, interferon gamma/alpha, IL2-STAT5 signaling, and IL6-JAKSTAT3 signaling were significantly activated in mHSPC_03 but not in mHSPC_02 (Fig. 3D, Supplementary Fig. 4A).
Supplementary Table 3 shows the clinical characteristics of the 2 patients with mHSPC. Histopathological examination of mHSPC_02 cells revealed a high mitotic rate, individual tumor cell necrosis, and increased expression of CDK1, an MNP cycling cluster marker gene (Fig. 3E). Immunohistochemistry analysis demonstrated the loss of NKX3.1, PSA, and ARs in mHSPC_02 cells (Supplementary Fig. 5A). The burden of lung metastasis observed through imaging was found to be significantly increased in mHSPC_02 after 2 months of treatment with AAP compared to baseline (Fig. 3F). In contrast, despite the high histological grade of mHSPC_03, there was no evidence of tumor cell necrosis or a high mitotic rate. Furthermore, tumor cells displayed low CDK1 expression, indicating reduced proliferative activity (Fig. 3G), and maintained NKX3.1 expression while exhibiting high levels of PSA and AR expression (Supplementary Fig. 5B). Notably, a partial response to AAP treatment was observed in mHSPC_03 compared to the baseline, as evidenced by the decrease in pelvic lymph nodes and multiple bone metastasis 3 months later (Fig. 3H).
The differential response to drug treatment was further investigated by conducting CellChat analysis to evaluate the cell-to-cell communication network between mHSPC_02 and mHSPC_03 (Fig. 4A). While mHSPC_03 exhibited a higher number of cell-cell interactions than mHSPC_02, the cells from the MNP cycling cluster in mHSPC_02 engaged in cell-to-cell communication with cells from the fibroblast, endo, basal, and T-cell clusters, whereas no such communication networks were observed in the MNP cycling cluster of mHSPC_03. Furthermore, specific reciprocal cell-cell interactions in 2 patients with mHSPC were identified, with the macrophage MIF signaling pathway emerging as the most significant pathway in mHSPC_02 (Fig. 4B). Gene expression level evaluation revealed that mHSPC_02 exhibited higher MIF gene expression than mHSPC_03 (Fig. 4C). Notably, cells in the luminal and MNP cycling clusters displayed higher MIF expression than other cell types (Supplementary Fig 6A and B), while the ligand-receptor pairs MIF−(CD74+CXCR4) and MIF−(CD74+CD44) were identified as the major signaling pathways (Fig. 4D). In contrast, CellChat analysis in mHSPC_03 identified the MHC-1 signaling pathway as the most significant (Supplementary Fig. 6B). Next, we evaluated the expression levels of MHC class 1 genes, including HLA-A, HLA-B, and HLA-C, and found a higher expression of these genes in mHSPC_03 than mHSPC_02 (Supplementary Fig. 6C). NK and T cells also showed a higher expression of MHC class 1 genes than other cell types (Supplementary Fig. 6D). Additionally, the ligand-receptor pairs HLA-A−CD8A, HLA-B−CD8A, and HLA-C−CD8A were identified as the major signaling pathways (Supplementary Fig. 6E). Furthermore, mHSPC_02 had a higher number of copy number alterations across the entire chromosome than mHSPC_03 (Fig. 4E). In particular, the majority of cells from mHSPC_02 exhibited a gain of chromosome 8q.
Next, we aimed to elucidate the potential correlation between intrinsically aggressive cells and disease progression. Notably, pseudotime trajectory analysis revealed that L0 cells, normal luminal subcluster from mHSPC_02 appeared at the end of the trajectory and eventually progressed to the L4 luminal subcluster (Fig. 4F). Therefore, we identified DEGs in the L4 and observed their association with processes such as response to unfolded protein, regulation cell death, and MAPK (mitogenactivated protein kinase) pathway including FOS, JUN, MYC genes (Supplementary Table 4). Conversely, the L0 cells from mHSPC_03 appeared at the end of the trajectory and eventually moved towards the both L3 and L5 normal luminal subclusters (Fig. 4G).
4. Identification of 2 mHSPC Molecular Subtypes Based on the 87-Network Gene Set
To systematically define the gene set associated with mHSPC, we constructed a network of DEGs in the luminal 8 subcluster (n=94), considered mHSPC luminal, and the MNP cycling cluster (n=118). Network analysis using the GeneMANIA plugin in Cytoscape revealed 87 genes associated with physical interactions, including MIF and CD74 genes, which are potentially indicative of an aggressive phenotype (Fig. 5A). The functional and clinical relevance of these 87 genes was further evaluated by performing ssGSEA of the WTS data from 52 mHSPC samples. The clinicopathological characteristics of patients in our cohort are presented in Supplementary Table 5. Two molecular subtypes identified among the mHSPC samples based on their ssGSEA scores were characterized by whether their scores were above or below the average ssGSEA of the 87 network genes. Subtype 1 was observed in 29 samples (56%), whereas subtype 2 was present in 23 samples (44%) (Fig. 5B). Overall, 15 pathway gene sets, including the G2M checkpoint, MYC target, and PI3K/mTOR signaling pathways, were significantly activated in subtype 2 (Fig. 5C). Aggression-related gene signatures such as ONCO, TSG, CIN25, CIN70, hES1, hES2, and CRPC51, and MIF pathway that was associated with mHSPC_02 sample, were investigated to further expand our understanding of the molecular phenotypes within each subtype. We found that subtype 2 exhibited the highest expression of ONCO, CIN25, CIN70, hES1, hES2, CRPC51 and MIF, but not of TSG (Fig. 5D).
Next, we compared radiographic PFS (r-PFS), PSA-PFS, failure-free survival (FFS), and time to castration-resistant PCa (ttCRPC) between the subtypes. Patients with subtype 2 exhibited poorer survival than those with subtype 1, as indicated by the r-PFS (p=0.047), PSA-PFS (p=0.498), FFS (p=0.067), and ttCRPC (p=0.062) values (Fig. 5E). However, no significant differences were observed in clinical characteristics such as PSA, ECOG, Gleason score, and CHAARTED and LATTITUDE risk criteria between subtypes 1 and 2 (Supplementary Table 6), except for the proportion of the PAM50 classification. Although we assumed that the PAM50 subtypes and AR activity in our cohort would validate survival outcomes, these factors were not associated with clinical outcomes in our dataset (Supplementary Figs. 7, 8). The prognostic significance of the 87 gene sets was assessed using ssGSEA scores to categorize the gene sets into high and low groups. High scores were correlated with DFS in both TCGA-PRAD (p=0.0003) and GSE21032 (p=0.0102) datasets (Fig. 5F), but not with OS in the TCGA-PRAD dataset (p=0.096, Supplementary Fig. 9A). Univariate and multivariate analyses of FFS confirmed the prognostic significance of molecular subtype (hazard ratio [HR], 1.97; 95% confidence interval [CI], 0.9356-4.147; p=0.0743; HR, 3.089; 95% CI, 2.9958-5.959; p=0.0092; Supplementary Table 7).
5. Generation and Validation of a 14-Gene Prognostic Signature
We subsequently attempted to narrow down a gene expression signature that can be applied in real-world clinical practice. A total of 14 commonly upregulated genes were identified by comparing the 87 network genes identified from scRNA-seq (n=87), CRPC-like cell marker genes (n=51), and DEGs from subtype 2 (n=421) (Supplementary Fig. 10A). This signature was validated by applying it to over 1,200 PCa samples from 8 independent public datasets using ssGSEA. Distinct expression patterns were observed in each tissue type (Supplementary Fig. 10B); these genes exhibited tumor-specific expression, which gradually increased as the disease progressed, compared to nontumor tissues. These findings confirmed the prediction robustness of the 14-gene signature, which was further validated using 2 independent cohorts (Supplementary Fig. 10C). The signature exhibited excellent predictive performance for survival outcomes, consistently demonstrating that patients with high expression of the 14-gene prognostic signature had significantly shorter survival than those with low expression (Supplementary Fig. 10C), but not OS in TCGA-PRAD (p=0.089) (Supplementary Fig. 9B). Moreover, patients with high signature expression exhibited poorer survival than those in the low-expression group, as evidenced by the r-PFS (p=0.018), PSA-PFS (p=0.399), FFS (p=0.034), and ttCRPC (p=0.028) values (Supplementary Fig. 10D).
6. The 14-Gene Signature Is Associated With Drug Treatment
We observed different patterns in the survival outcomes between the 2 subtypes. Patients in the low-expression group of the signature, characterized by a less aggressive phenotype, showed more favorable survival outcomes following AAP treatment, whereas they had poorer survival outcomes after early DCT treatment (Supplementary Fig. 11). In contrast, survival outcomes depended on the treatment modality in patients with a high signature expression, characterized by an aggressive phenotype. Specifically, patients had more favorable survival outcomes after early DCT treatment and poorer survival outcomes following AAP treatment (Supplementary Fig. 12).
DISCUSSION
To the best of our knowledge, this is the first comprehensive transcriptomic analysis of patients with mHSPC who received early DCT and AAP therapies. The transcriptomic features associated with drug responses and survival outcomes were investigated using scRNA-seq and WTS data. A major highlight of this study was the characterization of luminal subclusters, including normal and tumor luminal cells. In particular, subclusters L3 and L5 exhibited lower enrichment scores for oncogenes and higher enrichment scores for TSGs than the other luminal subclusters. Moreover, biological processes in L3 and L5 cells related to homeostasis and response to ions—specifically metal ion homeostasis, response to zinc ions, and response to copper ions—were activated. Previous studies have reported abnormal zinc levels in various cancers [11,12], while disruption of intracellular copper homeostasis has been shown to inhibit PCa cell growth [13]. We observed an activation of homeostatic processes in normal luminal subclusters, suggesting that cells in L3 and L5 may directly affect prostate Zn and Cu levels through their homeostatic functions. These findings highlight the importance of understanding the molecular characteristics and interactions within luminal subclusters in the context of mHSPC. The elucidation of the role of homeostatic processes and ion responses in these subclusters will provide insights into the potential mechanisms underlying PCa development and progression.
Single-cell transcriptomic profiles and existence of the MNP cycling cluster in 2 patients exhibited significant differences in cell-cell interactions, CNV, and inferred pseudotime trajectory analysis, indicating a distinct molecular signaling pathway that can be utilized in tailored treatment strategies based on these molecular subtypes. Furthermore, gene set analysis of both the WTS and scRNA-seq data identified commonly activated genes associated with cell cycle regulation [14], MYC targets [15], E2F targets [16] and hES cell-related genes [17]. The activation of these genes disrupts tightly coordinated network systems, contributes to cancer development, and indicates aggressive phenotypes and unfavorable clinical outcomes.
Comparison of 3 gene sets revealed a 14-gene signature that is crucial for mHSPC and PCa development. The 14 specific genes were involved in different aspects of tumor progression and appeared to be related to the biological features of aggressive tumors. Being essential genes associated with cell proliferation, MKI67, TOP2A, CDC20, and CDK1 were expressed in nearly all phases of the cell cycle [18-21]. Moreover, TOP2A expression was associated with an aggressive PCa subgroup [22], and NUSAP1 detected in recurrent PCa has been implicated in cancer cell proliferation and invasion [23]. BIRC5 is associated with proliferative activity and is a potential cancer therapeutic target [24]. CENPF interacts with FOXM1 and synergistically promotes tumor growth [25], while high expression of ASPM is associated with the promotion of stemness by augmenting catenin signaling [26]. Gavish et al. [27] recently provided a detailed map of transcriptional intratumor heterogeneity (ITH) and identified the cell cycle as one of the most commonly observed meta-programs across different cancer types. They suggested that cycling cells might suppress other molecular programs and allocate resources towards proliferation. We also observed the involvement of the cell cycle program, which is characteristic of ITH, and found that its activation in patients with mHSPC plays a significant role in promoting the aggressive phenotype.
Our study has several limitations. First, the number of patients with mHSPC in the scRNA-seq data was relatively small. Second, without exome sequencing data, transcriptome analysis alone cannot identify the distinct features of mHSPC. Third, the short-term follow-up period may have underestimated the potential benefits of the treatment. Despite its limitations, this is the first study to provide an indepth evaluation of the biological nature of mHSPC based on single-cell resolution.
CONCLUSIONS
In sum, this study provides valuable insights into the identification of high-risk patients, novel biomarkers, and potential therapeutic targets for individuals with mHSPC. Furthermore, the results in this study can serve as a basis for future investigations aimed at refining prognostic strategies and developing targeted therapies for patients with mHSPC.
PDF Links
PubReader
ePub Link
XML Download
Full text via DOI
Download Citation
Supplement1
Print
