Integrated single-cell transcriptomic analyses identify a novel lineage plasticity-related cancer cell type involved in prostate cancer progression

Abstract

Background: Cancer cell plasticity is the ability of neoplastic cells to alter their identity and acquire new biological properties under microenvironmental pressures. In prostate cancer (PCa), lineage plasticity often results in therapy resistance and trans-differentiation to neuroendocrine (NE) lineage. However, identifying the cancer cells harboring lineage plasticity-related status remains challenging.

Methods: Based on 13 multi-center human PCa bulk transcriptomic cohorts (samples = 3314) and 9 bulk transcriptomic datasets derived from PCa experimental models, we established an integrated lineage plasticity-related gene signature, termed LPSig. Leveraging this gene signature, AUCell enrichment analysis was applied to identify the cell population with high lineage plasticity from a comprehensive single-cell RNA-sequencing (scRNA-seq) meta-atlas assembled by us, which consisted of 10 public human PCa scRNA-seq datasets (samples = 93, cells = 222,529). Moreover, additional scRNA-seq dataset of human PCa, multiplex immunohistochemistry staining for human PCa tissues, in vitro and in vivo functional experiments, as well as qPCR and Western blot analyses were employed to validate our findings.

Findings: We found that LPSig could finely capture the dynamics of tumor lineage plasticity throughout the progression of PCa, accurately estimating the status of lineage plasticity. Based on LPSig, we identified a previously undefined minority population of lineage plasticity-related PCa cells (LPCs) from the human PCa scRNA-seq meta-atlas assembled by this study. Furthermore, in-depth dissection revealed pivotal roles of LPCs in trans-differentiation, tumor recurrence, and poor patient survival during PCa progression. Furthermore, we identified HMMR as a representative cell surface marker for LPCs, which was validated using additional scRNA-seq datasets and multiplexed immunohistochemistry. Moreover, HMMR was transcriptionally inhibited by androgen receptor (AR), and was required for the aggressive adenocarcinoma features and NE phenotype.

Interpretation: Our study uncovers a novel population of lineage plasticity-related cells with low AR activity, stemness-like traits, and elevated HMMR expression, that may facilitate poor prognosis in PCa.

Funding: This work was supported by National Key R&D Program of China (2022YFA0807000), National Natural Science Foundation of China (82160584), Advanced Prostate Cancer Diagnosis and Treatment Technology Innovation Team of Kunming Medical University (CXTD202216), and Reserve Talents of Young and Middle-aged Academic Leaders in Yunnan Province (202105AC160013).

Keywords: HMMR; LPCs; Lineage plasticity; Neuroendocrine prostate cancer; Single-cell.

Fig. 1

Overall methodology of this study and meta-analysis across multiple published scRNA-seq datasets of PCa. A and B. The experimental scheme for discovering (A) and validating (B) lineage plasticity-related gene signature (LPSig) and lineage plasticity-related epithelial cells (LPCs) in PCa.

Fig. 2

Identification of an integrated lineage plasticity-related gene signature (LPSig) in PCaProfiler bulk RNA-seq dataset. A. Heatmap showing the Spearman’s correlation between module eigengenes and lineage plasticity related traits. SCC, Spearman’s correlation coefficient. p-value was adjusted by false discovery rate (FDR). B. Univariate Cox regression analysis for progression-free survival (PFS) in TCGA PRAD (tumor = 499) cohort or overall survival (OS) in SU2C CRPC/Met (tumor = 106) cohort. Data are presented as log10 hazard ratio (HR) ± 95% confidence interval (CI). p-value was adjusted by FDR. The proportional hazard assumption was validated for all variables using Schoenfeld residuals test. C. The high Spearman’s correlation between gene significance and module membership (MM) in the pink module (left panel) or tan module (right panel). Dots within the red rectangle were defined as lineage plasticity-related genes, with both high gene significant and MM. D. Pearson’s correlation analysis between LPSig scores and the activities of multiple lineage plasticity-related features in 12 bulk transcriptomic cohorts of PCa. Positively correlated features are depicted in yellow while negatively correlated features are depicted in blue. The size of each circle is proportional to the magnitude of the absolute value of the correlation coefficient, with the total area of all circles normalized to 1. E. Violin plot showing the LPSig scores between mutation (Mut) and wildtype (WT) alleles in RB1 (Mut, n = 13; WT, n = 315; left panel) or TP53 (Mut, n = 119; WT, n = 209; right panel) in SU2C CRPC/Met cohort. F. Comparison of LPSig scores between wild-type LNCaP cells (LNCaP^WT, n = 2) and their counterparts with a dual knockout of TP53 and RB1 (LNCaP^DKO, n = 5), as well as between LNCaP WT (n = 3) and the enzalutamide-resistant LNCaP cells (LNCaP-Resist, n = 3). G. The distribution of LPSig scores across various PDX tumor groups with differing androgen receptor (AR) and neuroendocrine (NE) statuses from GSE199596 dataset (n = 112). ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. Student’s t test (E, F). One-way ANOVA followed by Bonferroni post hoc test (G).

Fig. 3

Identification and characterization of lineage plasticity-related epithelial cells (LPCs) in a large scRNA-seq meta-atlas of human PCa. A. UMAP showing the distribution of seven major cell types (left panel) and LPSig scores (right panel) in scRNA-seq meta-atlas of PCa (total cells = 222,529). B. Re-cluster analysis of 89,055 epithelial cells. C. Violin plots showing the expression levels of representative gene markers associated with luminal, proliferation, stemness, and lineage plasticity as well as NE features across different epithelial clusters. D. GSEA analysis revealing activated pathways in different epithelial cell types. NES, normalized enrichment score. FDR, false discovery rate adjusted p-value. E. Violin plots showing the enrichment scores for AR signaling, stemness, and NEPAL signature in various epithelial clusters. Error bars indicate the SD. ∗∗∗∗p < 0.0001. Kruskal–Wallis tests followed by Dunn’s post hoc tests with Bonferroni correction. F. UMAP plot showing four clusters of LPCs and NE cells. Cells are colored by cell subtypes (top left), dedifferentiation scores inferred by CytoTRACE (top right), and pseudotime inferred by monocle3 (bottom panel). G. Spearman’s correlation analysis between pseudotime scores and luminal feature, NEPAL signature, AR activity, and CytoTRACE dedifferentiation scores. The smooth curves in each plot, generated using locally weighted smoothing (LOESS) method, illustrated the relationship between pseudotime and the respective feature scores. The shaded areas of the curve represented ±95% confidence interval (CI).

Fig. 4

LPCs are correlated with PCa progression and poor disease outcomes. A. Tissue prevalence of each epithelial cluster estimated by Ro/e analysis. N, normal prostate; Adj, adjacent non-tumor tissue; Pri, primary HSPC; mCRPC, metastatic CRPC. +++, Ro/e > 3; ++, 2 < Ro/e ≤ 3; +, 1 < Ro/e ≤ 2; +/−, 0.1 ≤ Ro/e ≤ 1; −, Ro/e < 0.1; in which Ro/e denotes the ratio of observed to expected cell number. One cluster was identified as being enriched in a specific group if Ro/e > 1. B. Bar plot (left panel) and box plot (right panel) showing the distribution of LPC frequency in each sample across different subtypes of human PCa. Kruskal–Wallis tests followed by Dunn’s post hoc tests with Bonferroni correction. C. Correlation analysis between the frequency of LPCs and NE cells. The smooth curve, generated using linear model, and the shaded areas of the curve represented ±95% CI. D. Forest plot showing Univariate Cox regression analyses of LPC gene signature for biochemical recurrence (BCR) in primary PCa cohorts or overall survival (OS) in CRPC/Met cohorts. Data are presented as hazard ratio (HR) ± 95% confidence interval (CI). The p-value is FDR adjusted and adjust p-value < 0.05 was considered statistically significant. The proportional hazard assumption was validated using Schoenfeld residuals test. E. Kaplan–Meier survival analyses based on the estimated LPC frequency using CIBERSORTx algorithm in TCGA PRAD (n = 464) and SU2C CRPC/Met (n = 106) cohorts. F. Scatter diagram showing cell surface markers ordered by their Avg_log2FC in LPCs versus other epithelial cells. G. Violin plot showing the expression level of HMMR across different cell types in the scRNA-seq meta-atlas of human PCa. H. Representative images of mIHC staining for HMMR (green), CHGA (red), AR (pink), and DAPI (blue) in adjacent normal tissue, tumor samples with varying Gleason scores (GS), and CRPC samples. Scale bar, 40 μm. I. Dot plots comparing the AR expression between HMMR-positive and HMMR-negative tumor cells across different PCa subtypes. Each dot represents individual AR-positive cell with HMMR negative (Neg) or positive (Pos) expression. ∗∗∗p < 0.001. Wilcoxon rank-sum test. J. H&E staining and heatmaps of the spatial distribution of AR signaling activity or HMMR expression from the GSE230282 (n = 1).

Fig. 5

HMMR is transcriptionally repressed by AR. A. Dot plots showing the Pearson’s correlation between the expression of HMMR and indicated genes or other pathway activities in multiple cohorts of human PCa. B and C. Bar plots showing the changes in RNA expression of HMMR, KLK2, KLK3, and NKX3-1 in LNCaP cells (n = 3 for each group) treated with R1881 (B) or enzalutamide (C) compared to control groups at indicated time points. Error bars denote SD. ∗p < 0.05. Student’s t-test with Bonferroni correction for multiple comparisons. D. Western plot analyses of PSA (encoded by KLK3) and HMMR in LNCaP cells treated with R1881 (left panel), enzalutamide (middle panel) or DMSO at indicated time points, as well as AR and HMMR in LNCaP, DU145, and 22Rv1 cells (right panel). PCa cells were synchronized in the G1/S-phase by double thymidine (2 mM) block before treatment. E. Scatter diagram showing AR direct target genes ordered by their transcriptomic Log Fold Change (Log FC) in LNCaP cells treated with 5α-dihydrotestosterone (DHT) compared to control groups in GSE7868 dataset. AR direct target genes were identified by AR ChIP-seq analysis (binding peaks number ≥ 1). F. Genomic browser representation showing AR binding in HMMR intron 1, which encompasses an ARE, from LNCaP cells (GSE161167) treated with R1881 compared to control cells (n = 3 for each group). The nucleotides identical to the canonical ARE are underlined. JASPAR relative score represents the similarity level between the indicated sequence and canonical ARE, as predicted using the JASPAR web server. nARE, negative androgen response element. G. ChIP-qPCR of AR occupancy at sequences of the HMMR nARE-centric intron 1 and transcription start site (TSS), along with the KLK3 promoter region as a positive control, in LNCaP cells (n = 3 for each group) treated with R1881 for 48 h. Error bars denote SD. ∗p < 0.05. Student’s t-test with Bonferroni correction.

Fig. 6

Knockdown of HMMR represses the aggressive behavior of advanced PCa cells. A. Volcano plot showing differentially expressed genes (DEGs) between samples with high HMMR expression (HMMR^high, n = 164) and samples with low HMMR expression (HMMR^low, n = 164) in the SU2C CRPC/Met cohort. Up-regulated DEGs (adjusted p-value < 0.01 and LogFC ≥ 0.4) are depicted in red, while down-regulated DEGs (adjusted p-value < 0.01 and LogFC ≤ −0.4) are depicted in blue. Representative genes for indicated signaling or cell types are represented in different colors: AR targets are highlighted in green, NEPC drivers in orange, markers of LPCs in lightcyan, proliferation-related markers in purple. B. Bubble diagram showing pathways of activation or inhibition between HMMR^high and HMMR^low samples in SU2C CRPC/Met cohort through GSEA analysis. C. Bar plots showing the changes in RNA expression of HMMR, KLK2, KLK3, NKX3-1, EZH2, and AURKA in 22Rv1 cells (n = 3 for each group) following transfection with HMMR siRNA (siHMMR) compared to control (siNC). D. Bar plots showing the changes in RNA expression of NEPC differentiation markers CHGA and SYP in 22Rv1 cells (n = 3 for each group) following transfection with siHMMR compared to siNC. E. Western plot analyses of HMMR, AR, PSA (encoded by KLK3), and CHGA in 22Rv1 cells following transfection with siHMMR compared to siNC. F. Representative images (upper panel) and quantification results (lower panel) of mIHC staining for HMMR (green), AR (pink), CHGA (red), and DAPI (blue) in 22Rv1 cells treated with siHMMR or siNC (n = 4 for each group). Scale bar, 20 μm. G. Cell counting assays of DU145 (left panel) and 22Rv1 cells (right panel) treated with siHMMR or siNC (n = 4 for each group) at indicated time points. H. Representative images of colony forming in PCa cells transfected with siNC or siHMMR (n = 3 for each group). Colonies were visualized after 2 weeks. I and J. Transwell assays of migration and invasion in DU145 (I) and 22Rv1 cells (J) transfected with siHMMR or siNC (n = 3 for each group). After 24 h, the number of cells crossing the chambers was assessed. Scale bars, 50 μm. Error bars denote SD. ∗p < 0.05; ∗∗p < 0.01. One-way ANOVA were used to conduct difference comparisons of three groups (C, D, F–J).

Fig. 7

HMMR silencing suppresses NEPC tumor growth in mice. A. Western plot analyses of HMMR, AR, PSA (encoded by KLK3), and CHGA in 22Rv1 cells with long-term HMMR silencing using shRNA lentivirus (shHMMR) at 4 weeks, compared to the control group (shNC). B. Representative images (left panel) and quantification results (right panel) of mIHC staining for HMMR (green), AR (pink), CHGA (red), and DAPI (blue) in 22Rv1 cells infected with shHMMR or shNC lentivirus at 4 weeks (n = 4). Scale bar, 40 μm. C. 22Rv1 tumor growth curves of nude mice inoculated with shRNA-mediated HMMR silencing 22Rv1 cells or control cells (n = 6 for each group). D. Representative image of subcutaneous xenograft assay using shRNA-mediated HMMR silencing 22Rv1 cells and respectively control cells (n = 6 for each group). E and F. Representative images (E) and quantification results (F) of IHC staining for HMMR, Ki67, cleaved Caspase-3 (c-Caspase 3), AR, PSA, and CHGA in 22Rv1 xenografts tumor samples (n = 3 for each group). Scale bar, 20 μm. Error bars denote SD. ∗p < 0.05; ∗∗p < 0.01. Student’s t test (B, C, and F).

Fig. 8

Schematic summary of the findings from the present study.