Categories
Prostate cancer

Spatio-temporal analysis of prostate tumors in situ suggests pre-existence of treatment-resistant clones

Heterogeneous response of prostate tumor cells to castration

We have previously demonstrated that ADT overall reduces the expression of non-homologous end joining repair proteins, such as Ku70 and P-DNA-PKcs, and thereby hampers the DNA repair capacity of PCa cells, which agrees with increased response to radiation therapy54. However, we noticed an extensive variability in the levels of these repair proteins in the biopsies, both within and between patients.

To further investigate this, we analyzed a set of formalin-fixed needle biopsies pre- and post-ADT from five patients with Gleason scores ranging from six to nine (GG 1–5) and analyzed repair protein expression patterns along with nuclear AR, using immunohistochemistry (IHC). Epithelial nuclei were distinguished by an in-house-developed software (Fig. 1a), and epithelial nuclear AR intensities were measured from biopsies pre- and post-ADT, showing the presence of nuclear AR also post-ADT in a fraction of the epithelial glands, in all five patients, which indicates treatment unresponsiveness (Fig. 1b, c).

Fig. 1: IHC images for the AR activity in epithelial nuclei performed on needle biopsies pre- and post-ADT.
figure 1

a Prostate tissue from a paraffin-embedded biopsy stained for the AR (red) and DNA (blue). AR levels in epithelial cell nuclei pre- and post-ADT were extracted, encircled in white (right panel; for details see Online methods). Scale bars 100 µm. Insets 1–3 show examples of occasional misclassifications that may occur using the algorithm. Inset 1 shows failure of separating two nuclei. Inset 2 shows failure of nucleus detection due to low DNA-staining intensity. Inset 3 shows misclassification of stromal cell nuclei as epithelial, due to their proximity to epithelial cells. b Biopsies pre- and post-ADT from patient B. Small clusters of cells with nuclear AR still exist after ADT. Scale bar in whole figure, 1 mm, scale bars in close up, 100 µm. c Nuclear AR intensity from five patients with Gleason scores ranging from 6 to 9 (GG 1–5). The x-axis indicates the positional rank of each cell in the biopsy’s length direction while the y-axis shows the mean intensity of the androgen receptor in each cell’s nucleus. Source data is provided as Source Data file. AR androgen receptor, ADT androgen deprivation therapy, P X patient X, GG grade group, GS Gleason score.

We found evidence that high levels of nuclear AR predict high levels of said repair proteins and a correlation between levels of nuclear AR and serum PSA post-ADT (Supplementary Fig. 1; Supplementary Table 1). Non-responsiveness, was thus strongly linked to the expression of nuclear AR post-ADT. ADT resistance may also be the result of the activation of other genes or oncogenic pathways. To gain a broader molecular understanding of implicated genes and pathways contributing to ADT resistance in situ we continued with an exploratory transcriptome analysis.

Design of spatial gene expression experiments

To chart the molecular landscape of non-responsive cell clusters, fresh core needle biopsies pre- and post-ADT were collected from the prostates of three patients. The biopsies were snap frozen to facilitate gene expression analysis. We investigated the spatial gene expression profiles of each biopsy using ST50, a transcriptome-wide methodology capturing mRNA from tissue sections using spatially barcoded spots on microscopic slides, followed by a data driven analysis that identifies gene expression patterns across the tissue sections (overviewed in Supplementary Fig. 2).

In total eight core needle biopsies per patient were analyzed; four biopsies pre-ADT and four biopsies eight weeks after gonadotropin-releasing hormone (GnRH)-analogue treatment (post-ADT). Clinical data during this period were collected (Supplementary Table 2). Sections from each biopsy were hematoxylin and eosin (H&E) stained, scanned, independently annotated by two pathologists, and analyzed with ST. Adjacent tissue sections were immunostained for AR to enable a comparison to in situ gene expression patterns. Spatially variable high levels of nuclear AR were observed in the new set of fresh frozen biopsies, validating previous observations (exemplified in Fig. 2a).

Fig. 2: Overview of AR activity and bulk analysis results pre- versus post-ADT.
figure 2

a Nuclear AR activity pre- and post-ADT identifies responding and resistant prostate cancer cells post-ADT (green and yellow arrows, respectively). The scale bars in the upper panels are 1 mm, in lower panels 200 µm. b Principal component analysis of pseudo-bulked spots per needle biopsy pre- and post-ADT, after rlog transformation. Source data is provided as Source Data file. c DGE analysis (Welch’s t-test) of epithelial spots post- versus pre-ADT in patient 1. AR regulated genes are significantly down-regulated post-ADT (magenta). Genes shaded gray are relevant to PCa or genes of interest. Source data is provided as Source Data file. d Activated pathways for genes upregulated post- and pre-ADT (q < 0.05), color scale represents significance. ADT androgen deprivation therapy, AR androgen receptor, DGE differential gene expression.

The histological annotations of the tissue sections were conducted using the Gleason grading system. Information regarding immune response/inflammation and high-grade prostatic intraepithelial neoplasia (HGPIN) was also notated (Supplementary Fig. 3). In addition, the cells covering the spatially barcoded spots were further categorized into four tissue types: (i) stroma, (ii) 1–10% epithelium, (iii) 11–50% epithelium, (iv) 51–100% epithelium. An overview of the collected data, annotations, and performed analyses is illustrated in Supplementary Fig. 4.

Each patient was assigned to a representative clinical response group regarding ADT-treatment: responder (patient 1), moderate responder (patient 2), and non-responder (patient 3), based on the clinical data (Table 1, Supplementary Table 2), such as PSA nadir, PSA progress, and metastasis status.

Table 1 Clinical data for the three patients used in the study

Gene expression analysis in situ before and after ADT treatment

Gene expression analysis in tissue sections was achieved using barcoded arrays with spots on the surface, each with a known x- and y-coordinate. Each spot had a diameter of 100 micrometers, capturing the transcriptome from around 10–50 cells. The sample handling protocol for prostate tissue55 was adjusted to be compatible with the limited amount of tissue provided by core needle biopsies (Supplementary Figs. 56). On average, we detected approximately 4000 expressed genes per spot for all biopsies (Supplementary Fig. 6c).

First, to obtain an overview of patients and their biopsies, we bulked all spots per biopsy into individual pseudo-bulk samples and performed a principal component analysis (PCA). Most of the biopsies separated patient-wise (Fig. 2b) as has previously been observed in patients with PCa56. For patient 1 (responder), biopsies separated pre- and post-ADT. Patient 3 (non-responder) exhibited only a small separation, while patient 2 (moderate responder) showed more extensive spread between the biopsies.

To further assess the overall gene expression differences, we next conducted differential gene expression (DGE) analysis between histological areas comparing pre- and post-ADT samples. By merging gene expression data from epithelial and stromal spots, respectively, followed by DGE analysis, we could compare the temporal changes between the two histological entities. AR-regulated genes, such as KLK3, KLK2, and NKX3-157, were downregulated post-ADT in epithelial spots, for the three patients (Fig. 2c, Supplementary Fig. 7a–d). In line with clinical responsiveness, patient 1 had more differentially expressed AR-regulated genes compared to the other patients, indicative of successful and persistent downregulation of AR by ADT.

To identify biological pathways associated with differentially expressed genes (DEGs, q < 0.01), we performed functional enrichment analysis, querying against the KEGG database58 (Fig. 2d, Supplementary Fig. 7e–g). Among the three patients, pre-ADT, pathways related to AR (e.g., PPAR signaling, steroid hormone biosynthesis, sphingolipid signaling) and protein processing (e.g., protein processing in endoplasmic reticulum, lysosome, phagosome) were activated. Pathways related to cell migration (regulation of actin cytoskeleton, focal adhesion) were activated post-ADT for all patients. Subsequently, the same analysis procedure was performed for stroma demonstration and an upregulation of immune response genes post-ADT in all patients was observed (Supplementary Fig. 8).

Spatially resolved transcriptomes of patient biopsies

Next, we investigated the spatially resolved transcriptomes for the study cases. In total, we generated spatial and transcriptome-wide data for more than 4000 barcoded spots from 48 core needle biopsy sections including two consecutive tissue sections per biopsy. ST data from patients were analyzed individually by applying Spatial Transcriptome Decomposition (STD)52 (Supplementary data files 16; schematically overviewed in Supplementary Fig. 2). STD is a probabilistic model that factorizes the observed transcript data into latent gene expression factors (Methods). The factors characterize distinct metagenes, groups of genes that are likely to be co-expressed, and their spatial expression patterns.

The patient-specific factor analysis provided 13, 16, and 10 gene expression factors for patient 1, 2, and 3, respectively (Methods), and a control was made to ensure the independence between the factors (Supplementary Fig. 9). The ranked list of marker genes within each factor was used for molecular annotation. We broadly categorized factors into three entities; stroma, immune enriched stroma, and tumor activity (Supplementary Figs. 1012). Factors representing normal epithelium were not identified. Overall, the molecular annotation of the factors overlapped with the histological annotations, but the factor analysis provided an increased resolution in the analysis. For example, multiple tumor factors were identified enabling the temporal subdivision into responding or non-responding factors. These factors were then compared with AR staining in adjacent sections.

Hierarchical clustering of the factors was performed for the patients (Supplementary Figs. 1316) providing a summary of the factors, the corresponding histological entities and the molecular annotations for the tumor factors. The overview also depicts the annotation of temporal differences in AR staining.

The UMAP visualization of factor activities and its spatial position in the core needle biopsies is shown in Fig. 3. The two-dimensional embedding (UMAP) presents the distinct factors that correspond to histological features such as stroma, and tumor. In addition to the temporal differences observed for the tumor factors, we identified stroma factors that also displayed temporal differences.

Fig. 3: UMAP visualization of gene expression factors (left) and the spatial position of gene expression factors (right).
figure 3

a Patient 1- clinical responder. b Patient 2 – moderate responder. c Patient 3 – non-responder. Clusters of spots from factors of the same type are encircled. UMAP Uniform Manifold Approximation and Projection.

The molecular heterogeneity as determined by the number of tumor factors and spread in UMAP space (left panels) appeared to be most prominent for patient 3. In patient 3, at least eight distinct tumor factors could be identified while patient 1 displayed five tumor factors. The UMAP results for patient 2 points to a broader representation, represented by eleven tumor factors, in line with the initial PCA data on pseudo-bulk data on tissue sections, which indicate a higher degree of tissue heterogeneity.

The spatial position of the tumor factors along the needle biopsies outlines the temporal differences after eight weeks of ADT, raising certain observations from our analysis; Patient 1, the treatment responder, has a reduction of tumor factors post-ADT, where they largely are replaced by stroma factors (right panel). Interestingly, even though this patient responded clinically, as determined by a persistent low PSA-value, we can still observe some tumor factors post-ADT. We cannot determine whether the post-ADT tumor cells that express these factors are in a transient mode of disappearing, still capable of having AR-activity, or represent pre-existing resistant cells.

In patient 2, the moderate responder with high initial and low PSA levels at week eight, we observe four responding and seven non-responding factors. However, we observe a major change in tissue makeup post-ADT, similar to Patient 1. Most of the tissue after treatment is composed of stroma and immune response genes.

Patient 3, the non-responder with moderate PSA levels before and after eight weeks, had no responding factor but eight non-responding factors. The spatial distribution and amount of the tumor factors were not changed for Patient 3, over the therapy period. Although patient 3 only displayed non-responding factors, we detected a down-regulation of AR-regulated genes when comparing epithelial spots pre- and post-ADT (Supplementary Fig. 7c). This suggests that spots designated to non-responding factors also contain a fraction of cells that do respond to ADT.

Overall, the ratio of responding and non-responding tumor factors after eight weeks of treatment agrees with the clinical outcome of ADT (Supplementary Fig. 13). In all three patients, the tissue area corresponding to identified tumor factors post-ADT, and thus, the number of cancer cells, decreased substantially due to the apoptotic effect of ADT59. Responding factors were overall present in larger groups of spatially confined spots pre-ADT than non-responding factors (pre-ADT) (Supplementary Figs. 1012). This is also true for the non-responding tumor factors, which are generally present in larger areas pre-ADT compared to a more dispersed pattern post-ADT, except for patient 3 (Fig. 3, Supplementary Fig. 17).

Further, we observed in the PCA analysis (Fig. 2b) that biopsy 1 (pre-ADT) for patient 2, overlaps with non-responding factors in patient 3, which might indicate presence of a resistant phenotype in patient 2.

Neuroendocrine differentiation (NED) can occur in prostate cancer. Prostatic adenocarcinomas that have undergone NED are resistant to ADT. To investigate if there was an enrichment of cells that had undergone NED in the areas expressing resistant factors, we stained all biopsies in patient 2, before and after ADT, with the neuroendocrine marker chromogranin A (CgA). Areas expressing a given factor were mapped against the corresponding area in the CgA stained section, and the ratio of CgA positive cells for the areas was quantified. The criteria for the area depicted was that at least five spots with a given factor should be clustered together in both ST-replicates. No difference in ratio of CgA positive cells was found in areas expressing resistant factors compared to areas expressing non-resistant factors. (Supplementary Figs. 1819, Supplementary Tables 34).

To interpret the spatial RNA data obtained from the ST method it is of importance to know how well it correlates with the corresponding protein levels in the tissue. In a large study, including expression data from 60 genes in several different tissues and cell-lines, the correlation between the number of RNA transcripts and number of protein molecules was good within each gene, but less accurate when comparing the number of transcripts with the number of protein molecules in-between different genes60. In accordance with this result, we showed in a previous article, that there is a good concordance between RNA expression detected with ST technology and protein levels detected with immunocytochemistry for all 7 proteins tested in prostate tissue applying similar ST protocol used herein55. Here we show that this relationship between RNA expression and protein levels also is valid for the androgen receptor (Supplementary Fig. 20).

Processes modified in non-responsive tumor factors

The factor annotation across the biopsy sections served as a tool to improve our investigation of the processes that underlie non-responsiveness. We, therefore, undertook a DGE analysis by first bulking the gene counts of the responding and non-responding factors, respectively, for patients 1 and 2. Here, spot selections were based on a stringent scheme: for example, spots with tumor factors above a defined factor intensity threshold were counted as non-responding spots if they consisted of histologically annotated tumor cells (see Methods for details, Supplementary Fig. 21a, b). Importantly, we only used information collected from spots pre-ADT to exclude ADT as a confounder in the DGE analysis.

The selection process resulted in an approximately equal number of responding and non-responding spots for patient 1 (71 and 90 spots, respectively, i.e., 56% non-responding spots from a total of 161 spots), while for patient 2 the fraction of non-responding spots was higher (64% non-responding spots from a total of 254 spots). The top differentially expressed genes (DEGs, q-value < 0.05 and logFC > 0.3 for patient 1 and > 0.5 for patient 2) between responder and non-responder spots were visualized as a gene expression heatmap (Fig. 4a, b). Only four DEGs were identified in patient 1, while patient 2 had more than 50 DEGs, potentially reflecting the clinical responsiveness of the patients.

Fig. 4: DGE analysis of responding versus non-responding prostate tumor factors.
figure 4

a Tissue regions in Patient 1 (q-value < 0.05, logFC > 0.3), pre-ADT and b tissue regions in Patient 2 (q-value < 0.05, logFC > 0.5). The top 50 DEGs are shown and an overlap with patient 1 of non-responding genes CD74 and HLA-DRA was observed. Scale-bars indicate normalized counts. Tissue regions were pre-ADT and a two-sided Wilcoxon rank-sum test adjusted for multiple comparisons was used. Source data are provided as Source Data files for a and b. c Gene expression of DHCR24 and SNHG25, upregulated in non-responding cells, plotted onto tissue sections pre- and post-ADT in patient 2. d Representative example of nuclear AR activity in an epithelial luminal gland (yellow arrow), post-ADT, at an area where the non-responding gene DHCR24 is highly expressed (marked with ‘1’ in the bottom row in c) in patient 2. This serves as a confirmation of the presence of non-responsive cells in these areas. Scale bar 50 µm. e Activated pathways for DEGs upregulated in non-responsive spots or upregulated in responsive spots (q < 0.05) in patient 2, color scale represents significance. Source data is provided as Source Data file. ADT androgen deprivation therapy, AR androgen receptor, pat patient, rep replicate, DEG differentially expressed genes, FC fold change.

To check that the DE-results of patient 1, which are immune related genes, are belonging to the non-responding factors and not simply is the result of non-factor-related cell types (since we are doing the analysis on all transcripts in the selected spots), we plotted all factors in patient 1 for their log2 fold changes for these specific genes (Supplementary Fig. 22). The result shows that the non-responsive factor 11 has a FC between 2 and 3 for all these four genes.

For patient 2, spots from factor 5, 9, 10 and 14 were excluded from the DGE analysis, because of not fulfilling the specified criteria, such as presence of factor activity in more than one biopsy (see Methods for details). The DEGs in patient 2 could be separated into two groups – one group of genes that are predominantly expressed in the non-responding spots and another group of genes with less distinct separation between responding and non-responding spots. The non-responding group of genes (including DHCR24, SNHG25, TRPM8, IFI6, FKBP2, PPDPF, HLA-DRA, TAGLN, TIMP1, AEBP1, IGFBP7, MGP, A2M, and CD74) has previously been described to correlate with e.g., androgen-independence, cell migration, resistance, PCa, cancer, tumor stroma, immune response, and inflammation61,62,63,64,65,66,67,68,69,70,71,72,73,74 (see Supplementary Tables 56 for complete list including downregulated genes). To spatially analyze the observations, the gene expressions of DHCR24 and SNHG25 were plotted onto the biopsy sections (Fig. 4c). Both genes are present over investigated time but more spatially prevalent before treatment onset.

To evaluate the set of predominantly non-responding genes, we took advantage of patient 3, the non-responder. Indeed, we confirmed that many of the marker genes (DHCR24, TRPM8, IFI6, H2AFJ) in non-responding areas in patient 2 were also expressed in the non-responder spots of patient 3 (Supplementary Fig. 23a). To ascertain that predicted areas with non-responding gene profiles are indeed non-responding we took advantage of the AR staining in the adjacent tissue sections, demonstrating nuclear staining for tumor cells in all areas post-ADT that had been annotated as non-responding. Some of these areas were crowded with eAR(+) nuclei and some had only sparse nuclei of this type. Non-responsive prostate glands and scattered cells, with nuclear AR, were found in areas where DHCR24 and SNHG25 were expressed, and a representative image of corresponding AR-stained tissue region shows a small AR-active prostate gland (Fig. 4d). Representative areas of AR-staining across the tissue sections post-ADT are shown in Supplementary Fig. 23b.

We continued by performing pathway analysis on all significant DEGs found in non-responding versus responding spots in patient 2 (Fig. 4e). Analysis revealed pathways associated with migration75,76,77, survival75, metastasis75,76,78,79 and androgen-independence75, proliferation80, and angiogenesis80 (focal adhesion, regulation of actin cytoskeleton, cell adhesion molecules (CAMs), proteoglycans in cancer, and platelet activation). Furthermore, pathways such as ECM-receptor interaction, MAPK signaling, RAS signaling, Rap1 signaling, PI3K-Akt signaling, and immune-related leukocyte transendothelial migration were more prominent in the non-responding spots.

The underlying non-responding genes give further insight into molecular processes. For example it has been demonstrated that primary tumor cells can stimulate platelets to get activated81 which leads to the release of a wide range of growth factors and cytokines, such as TGFβ182 which has been shown to induce metastasis83,84,85,86.

Platelet activation can also activate intracellular signaling cascades, such as p42 MAPK, which can stimulate proliferation, survival, adhesion and chemotaxis of hematopoietic cells81. Further, MAPK- and the PI3K-Akt pathways play a key role in apoptosis and bone metastasis87. MAPK-mediated phosphorylation of the nuclear receptor co-activator 1 (NCoA1, or SRC-1) may increase the coactivators affinity for AR, contributing to disease recurrence and CRPC88.

The identified RAS/MAPK-pathway in the non-responding regions has been shown to contribute to PCa progression and metastasis89, and is activated in 43% of primary tumors while in 90% of metastatic tissues90. RAS signaling has, in cell lines, shown to decrease androgen dependence and promote metastasis91. Several studies suggest that the PI3K-Akt pathway is involved in androgen-independent growth of PCa92,93,94,95,96,97, and genetic alterations of PI3K-Akt pathway occur in 100% of metastatic PCa which suggest a key role in the progression to CRPC12. The PI3K-Akt signaling pathway is also known to induce stem-like properties, proliferation, migration, angiogenesis, regulation of cellular growth and survival98, as well as having a potential correlation with PTEN-loss99.

The enrichment of ECM-receptor interaction in non-responding areas has previously been shown to play a key role in metastasis since it needs to allow for a CAM-mediated coordinated balance between adhesion and detachments of tumor cells100,101,102. Akt signaling is a master regulator when it comes to inducing EMT and cancer stem cell phenotype by the ECM, and this can be mediated by various focal adhesion proteins and lead to activation of e.g., NF-κB103,104. Focal adhesion formations transduce ECM signaling into the tumor cells and activate the PI3K-Akt pathway105.

Furthermore, the pinpointing of Rap1 signaling in non-responding areas is interesting as this has shown to induce cancer cell proliferation and disease progression in several cancer types106,107,108, and particularly in PCa its activation affects integrins important in migration, invasion, and bone metastasis109,110. Increased Rap1 activity correlated with high metastatic potential in both PCa cell lines and in vivo, implicating Rap1 could be of therapeutic importance for curing PCa110. The leukocyte transendothelial migration is an important step in the initiation of an inflammatory immune response and chronic inflammation, suggested to serve as an anti-cancer therapy111. Further, for invasive cervical cancer, KEGG pathway enrichment analysis has revealed activated pathways such as ‘focal adhesion’, ‘ECM-receptor interaction’ and ‘platelet activation’112.

We also observe that these resistant tumor areas have a lower cell cycle activity as compared to non-resistant areas before treatment onset, when comparing to a 71 cell cycle gene signature113 (Methods, Supplementary Fig. 24).

The tumor microenvironment in prostate cancer

Previous studies have shown that lack or low levels of nuclear AR in stromal cells (sAR(-)), adjacent to tumor cells, are observed in the context of high Gleason scores, and metastasis114,115. In contrast, normal stroma expresses nuclear AR. Across our tissue sections, we observed multiple regions of sAR(-) cells (Fig. 5a) that was compared to our factor analysis. Hereby, we could investigate the stroma of responding and non-responding tumor factors, respectively.

Fig. 5: AR(+) stroma versus AR(−) stroma in biopsies pre-ADT.
figure 5

a Representative images of sAR(+) (top) and sAR(-) (bottom), in patient 2. Note the AR(+) luminal epithelial cells of the gland (bottom). Scale bars 50 µm. b Heatmap of the 26 DEGs (two-sided Wilcoxon rank-sum test adjusted for multiple comparisons, q-value < 0.05, logFC > 0.5) in the stroma with and without nuclear AR across spots in patient 2 and 3, respectively. Source data is provided as Source Data file. c AEBP1 and TIMP1 expression in a tissue section from patient 2 (pre-ADT, biopsy 2, left). HE-image with encircled spots annotated to GG5 and cell type composition of epithelial cells >50% (right). d Pathways activated for DEGs upregulated in sAR(−) areas or upregulated in sAR(+) areas (q < 0.05), color scale represents significance. Source data is provided as Source Data file. GG Grade group, GS Gleason score, ADT androgen deprivation therapy, GG Gleason grade group, GS Gleason score, AR androgen receptor, sAR(+) stromal nuclear AR positive, sAR(−) stromal nuclear AR negative, DEG differentially expressed genes, FC fold change.

The proportions of stromal AR-staining in stromal nuclei is illustrated as a pie chart per factor for each of the patient in Supplementary Figs. 1416. We hypothesized that a responding tumor factor would have a higher extent of surrounding stromal AR-positive (sAR(+)) cells, while non-responding tumor factors would be encircled by sAR(−) cells.

All stromal tissue spots were annotated by their AR-content using a defined threshold where tissue signals below the cutoff were treated as AR(−) stroma. In short, spots, by visual detection, composed solely of nuclei that lacked AR were counted as sAR(−), and spots with nuclei with at least 50% nuclear AR, were counted as sAR(+) (described in detail in Methods). Areas with a mix of sAR(−) and sAR(+) spots were annotated to sAR(mix) (Supplementary Fig. 21c).

For patient 1, 72% of the 287 epithelial-containing spots with attributed non-responding tumor factors were associated with sAR(−), while 15% were associated with sAR(+) (the remaining spots contained a mixture of AR positive and AR negative cells). 27% of the 79 epithelial-containing spots with attributed responding tumor factors were associated with sAR(−), while 60% were associated with sAR(+) (the remaining spots contained a mixture).

For patient 2, >99% of the 424 epithelial-containing spots with attributed non-responding tumor factors were associated with sAR(−). 53% of the 61 epithelial spots with attributed responding tumor factors were associated with sAR(−), and 47% contained a mixture.

For patient 3, 81% of the 399 epithelial-containing spots with attributed non-responding tumor factors were associated with sAR(−), and 9% were associated with sAR(+). No responding factors were found in patient 3.

We continued to investigate sAR(−) areas in patient 2 and 3 in more detail. We sought to compare gene expression levels between the stromal tissue regions sAR(+) and sAR(−) located next to tumor areas, independent of responsiveness of potential nearby tumor factors. DGE analysis was performed on the expression data between spots belonging to each tissue type (Fig. 5b, Supplementary Fig. 21d). The fraction of sAR(−) spots was 42% of a total of 172 spots. 26 DEGs were identified of which 21 was upregulated in sAR(-) spots (Supplementary Table 7), for example the epithelial- mesenchymal transition (EMT)-associated genes COL1A1, COL1A2, COL3A1, BGN, POSTN, SPARC, and AEBP1116,117,118. BGN is also associated with poor prognosis and PTEN deletion119 and upregulation during tumor angiogenesis120, and SPARC promotes bone metastasis in prostate cancer121. Both SPARC and POSTN are glycoproteins important for the structural network in the ECM. POSTN has been shown, using in vitro models, to be upregulated in advanced stages of cancer stroma and in bone metastases, however not in advanced stages of tumor cells122, in line with our observations (Fig. 5b, Fig. 4b). Further, SFRP4 is a marker for aggressive PCa and also a post-surgery recurrent marker123, and TIMP1 expression has been shown to be elevated in PCa stroma, to stimulate cancer associated fibroblasts, and to promote tumor progression124. The genes AEBP1 and TIMP1 are plotted onto the tissue section of biopsy 2 from patient 2 (pre-ADT) (Fig. 5c), and a comparison with histology reveals that these stromal compartments are located adjacent to PCa cells annotated as GG5.

Finally, in sAR(-) regions, we found an upregulation of several interesting pathways (Fig. 5d), of which the four most prominent ones were also upregulated in the non-responding tumor factors (ECM-receptor interaction, focal adhesion, PI3K-AKT, and Proteoglycans in cancer; Fig. 4b), together with TGF-β. The ECM-receptor interaction pathway has previously been shown to correlate with high reactive stroma content in PCa125. Also, changes in proteoglycans in the tumor microenvironment occur during tumor progression and affect e.g. cell signaling, chemokines, growth factors, and apoptosis126. The fact that the TGF-β pathway was upregulated in the sAR(-) regions is of particular interest since this can be induced by platelet activation which was activated in the non-responding tumor areas. TGF-β signaling promotes tumor initiation, progression, metastasis, EMT, stroma-tumor crosstalk, inflammation, immune-response, and angiogenesis81,82,83,84,85,86,127,128. Also, upon dissemination to the bones, tumor cells activate osteoclasts to degrade the bone matrix and release the stored TGF-β, which in turn leads to enhanced tumor cell malignancy129.

In summary, albeit sparseness of the material in the needle biopsies, we observed that a majority of the stroma in proximity to non-responding cancer cells pre-ADT lacks AR to a higher extent, in all patients. We noted that patient 1, who clinically had the best response to the treatment, displayed more sAR(+) surrounding the cancer, independent of responsiveness, while patient 2 and 3, who developed CRPC, displayed a higher frequency of sAR(−) areas in proximity to the non-responding epithelial spots. Future validation of these findings is important to reveal biomarkers and drug targets connected to stromal changes during the development of CRPC.