Human prostate cancer bone metastases have an actionable immunosuppressive microenvironment

Abstract

Bone metastases are devastating complications of cancer. They are particularly common in prostate cancer (PCa), represent incurable disease, and are refractory to immunotherapy. We seek to define distinct features of the bone marrow (BM) microenvironment by analyzing single cells from bone metastatic prostate tumors, involved BM, uninvolved BM, and BM from cancer-free, orthopedic patients, and healthy individuals. Metastatic PCa is associated with multifaceted immune distortion, specifically exhaustion of distinct T cell subsets, appearance of macrophages with states specific to PCa bone metastases. The chemokine CCL20 is notably overexpressed by myeloid cells, as is its cognate CCR6 receptor on T cells. Disruption of the CCL20-CCR6 axis in mice with syngeneic PCa bone metastases restores T cell reactivity and significantly prolongs animal survival. Comparative high-resolution analysis of PCa bone metastases shows a targeted approach for relieving local immunosuppression for therapeutic effect.

Keywords: T cell exhaustion; bone marrow; bone metastasis; human; mouse model; prostate cancer bone metastasis; single-cell RNA sequencing; single-cell landscape.

Figure 1.. Transcriptional impact of PCa metastases on the BM.

(A) Schematic illustration of study design. (B) H&E-stained tissue sections of metastatic tumor from a representative patient with bone metastatic PCa (100x -top, 400x -bottom). (C) Sagittal T1 MRI imaging of the thoracic spine for each of those three cases showing tumor masses with spinal cord compression. (top) The T10 vertebral body (arrow) is structurally compromised by tumor. (bottom) An extradural tumor mass in the dorsal spinal canal at the T5 level (arrow). (D,E) The analysis integrates scRNA-seq from three fractions (Tumor, Involved, Distal) of 9 metastatic PCa patients (rows), and 7 Benign BM controls. The datasets were integrated to establish joint embedding and annotation (E). Projections of individual samples are shown in D. (F) Marker genes for major cell populations. (G) The difference in frequency of major subpopulations between Benign controls and cancer patients (all three fractions). (H) Compositional differences between Tumor and Involved/Distal fractions.

Figure 2.. Common axes of intratumoral heterogeneity are seen within the metastatic tumor populations.

(A) Transcriptional heterogeneity of tumor cells from different samples is visualized on a tSNE embedding. (B) Integration of tumor cells from different patients using Conos (same colors as in A). (C) Expression of high-loading genes for four independent components (ICs) is shown based on the ICA analysis of the integrated tumor cell populations, illustrating re-current aspects of transcriptional heterogeneity between patients (red -high expression; blue -low). (D) Expression of high-loading genes identified by ICA are shown on the combined Conos embedding B.

Figure 3.. Overview of compositional and transcriptional shifts associated with presence of metastases in the BM.

(A) Compositional shifts for the most affected cell types are shown as boxplots comparing proportion of major cell populations between Benign controls and different cancer fractions (see Fig. S3A for all populations). (*p<0.05, ****p<0.0001) for select pairs, based on a Wilcoxon rank sum test. (B) Average combined expression distances between different patients. (C) The boxplots show normalized expression distances (y axis) show the extent of transcriptional variation between Involved and Benign fractions relative to expression variation within the fractions. Significance is shown based on empirical randomization of sample group assignment. Boxplots in A-C include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range. (D,E) Overview of the GO BP category enrichment in the top 300 up-regulated (D) and down-regulated (E) genes in the Involved vs. Benign differential expression analysis performed for each major cell subpopulation. Gene sets showing significant enrichment (FDR <0.05) were clustered into 20 clusters according to the similarity of the participating genes.

Figure 4.. Myeloid populations shift towards inflammatory monocytes and repressive macrophages in the Tumor fraction.

(A) Detailed annotation of the myeloid subpopulations is shown on a myeloid-specific joint embedding (left), together with select gene markers (right). (B) Changes in the composition of the myeloid compartment between sample fractions is visualized as cell density on the joint embedding. (C,D) Average expression of monocyte (C) and macrophage (D) signature genes (see STAR Methods, GSE5099) in different myeloid subpopulations shown as boxplot. Statistical significance assessed by Wilcoxon test, reported at 99% reproducibility power (see STAR Methods). Boxplots include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range. (E) Barplot showing flow cytometry validation of the monocytes/macrophages abundance from three independent patients. (*p<0.05) by Wilcoxon rank sum test. (F) Heatmap shows average expression of select genes from different categories (rows) across different cell populations (top color bar, colors matching panel A) for each patient (columns). (G,H) Average expression of the M2 macrophage signature gene panel (G), and monocyte inflammatory gene panel (H) across different monocyte populations shown as boxplot. Boxplots include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range.

Figure 5.. Comparisons of myeloid states with other cancers.

(A) Joint alignment of myeloid cells from different tumor types including other bone metastases (2 lung cancer and 1 melanoma bone metastases. The same contours outline the major myeloid cell populations. (B,C) Average expression of monocytes inflammatory signatures (B) and M2 macrophage signatures (C) in different cancer types. Boxplots include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range. (D,E) Heatmap shows average expression of select genes from different categories (rows) across different cancer types (top color bar, colors matching panel A) for each patient (columns), respectively for TIM (D) and TAM (E). (F,G) Overall expression similarity of TIM (F) and TAM (G) populations found in different samples is visualized using multidimensional scaling.

Figure 6.. Tumor fraction shows increased abundance of exhausted cytotoxic populations and helper T cells.

(A) Detailed annotation of the T cell subpopulations, shown on a T-cell specific joint embedding. (B) Expression of key genes from different categories (rows) in CTL-1 and CTL-2 populations across different sample types. (C) Shifts in the T cell populations visualized as cell density on the joint embedding. (D) Boxplots of cell frequency changes of different subpopulations between sample fractions. Boxplots include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range. (E) Scatter plot showing correlation of cytotoxicity and exhaustion scores for CTL-2 in Benign and Tumor fraction. Statistical significance was assessed using Pearson correlation test. (F) Barplot showing flow cytometry validation of the NK cell abundance from three independent patients. (*p<0.05) by Wilcoxon rank sum test. (G) The abundance (x axis) of TAM populations is associated with the increased exhaustion of the CTL-2 populations (y axis) (upper). Correlation of TAM signature score and exhaustion signature score is shown as scatter plot for TCGA PRAD bulk RNA-seq data (lower). (H) The number of known cognate receptor-ligand pairs for which the ligand is expressed in the myeloid subpopulation (rows), and the ligand is expressed in the T cell population (columns). (I) Receptor-ligand channels connecting TIM and TAM populations to the T cells, based on additional filtering criteria (see STAR Methods). (J,K) Cluster-average expression of CCL20 (J) and CCR6 (K) is shown for different cell populations. Statistical significance was assessed using Wilcoxon test. Boxplots include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range.

Figure 7.. Disruption of the CCL20-CCR6 signaling axis relieves T cell exhaustion and improves survival.

(A) Representative images of bone metastases localization as shown by bioluminescence optical imaging in the mandible or in long bones one week post injection of RM1-BoM3 cells. (B) CCR6 knockout results in a statistically significant survival benefit in syngeneic mouse model of bone metastases (n=9 WT and n=13 CCR6-KO). Two independent experiments. (C) Block of CCL20 by antibody results in a statistically significant survival benefit in syngeneic mouse model of bone metastases (n=22 IgG and n=20 CCL20Ab). Two independent experiments. (D) CD3+ T cells were isolated from the BM and analyzed by single-cell RNA-seq (n=2 mice per group). Changes in the composition of the T cell compartment is visualized as cell density on the joint embedding. (E) Abundance of Exhausted CTLs populations is shown as a percentage of all CT in mice with BM metastases. (F) Abundance of IFNg+ CTLs is shown as a percentage of all CTL in mice with BM metastases. (G) Dilution of CFSE in sorted CTLs was analyzed by flow cytometry after 96 hours of stimulation and compared to baseline levels as an indication of T cell proliferation. (n= 4 naïve, n= 5 WT, n= 5 CCR6-KO). Boxplots in E-G include centerline(median), 25th to 75th percentiles; and whiskers are highest and lowest values no greater than 1.5x interquartile range. (H) The magnitude of CCL20 expression in the TIM population (y axis) is strongly correlated with the overall abundance of TIM cells, measured as a proportion of all myeloid cells (x axis).