A Transcriptome Atlas of a Cross Section of a Multifocal Prostate Cancer
Emelie Berglund1,§, Niklas Schultz2,§, Maja Marklund1,¤, Ludvig Bergensthråhle1,¤, Joseph Bergenstråhle1, Reza Mirzazadeh1, Ludvig Larsson1, Firaz Tarish2, Anna Tanoglidi3, Jonas Maaskola1, Thomas Helleday2 and Joakim Lundeberg1,*
1Department of Gene Technology, School of Engineering Sciences in Chemistry, Biotechnology and Health, Royal
Institute of Technology (KTH), Science for Life Laboratory, Solna, Sweden.
2 Division of Translational Medicine & Chemical Biology, Karolinska Institutet, Science for Life Laboratory,
Solna, Sweden.
3Department of Clinical Pathology, University Uppsala Hospital, Uppsala, Sweden. §These authors contributed equally to this work.
¤These authors contributed equally to this work.
*Correspondence should be addressed to J.L. (joakim.lundeberg@scilifelab.se)
SUMMARY
Understanding the heterogeneous molecular landscape of prostate cancer is fundamental to improve treatment of the disease. Here, we aim to provide a broad molecular view of a cross-section of a prostate organ. This study manifests several unique tumor gene expression subtypes, with adjacent stroma, occupying distinct and different anatomical regions. The interplay between multiple molecularly defined tumor gene expression factors and/or Gleason scoring and the corresponding tumor microenvironment was hereby studied in detail. We perform a molecular sub-categorization of the tumor microenvironment throughout the whole prostate, using three different molecular principles; (i) AR staining, since loss of stromal AR is directly proportional to the degree of differentiation (Gleason score), (ii) Masson staining for its role in marking reactive stroma, and (iii) spatial transcriptomics analysis. In particular, we performed a detailed analysis of spatial distribution in histologic sections at the invasive border of a tumor foci contrasting it to the tumor core. Here, we show spatially confined DEPDC1, CHN1 and CRISP3 upregulation at the invasive border with implications of signal transduction pathways such as; PTEN, TGF-beta Receptor Complex, NOTCH4, ERBB2, and VEGF signaling, some of which are druggable. Overall, our results show an unprecedented view of the molecular heterogeneity of a prostate cancer not evident by other means. Here, we reveal patient-specific gene expression hallmarks from low to aggressive cancer within the same cross section of a prostate.
INTRODUCTION
Prostate cancer (PCa) is the most commonly diagnosed non-skin cancer (Grönberg 2003, Siegel et al. 2019). If the tumor is organ-confined at the time of diagnosis, prognosis is often favorable (Chi et al., 2009). The serum biomarker, Prostate Specific Antigen (PSA), is often the first and only sign that cancer is present. The current approach for validation and detection of PCa in situ involves sampling 12 needle biopsies using transrectal ultrasound (TRUS) guidance. Clinical risk assessment is primarily determined by pathologic grade based on the Gleason histological
grading system, initially proposed in the 1960s and modified multiple times thereafter. The grade as reported by the Gleason system is based on the combination of the most dominant and second most dominant histological patterns observed (Gleason et al., 1966, Gordetsky and Epstein, 2016). Since 2014, a new contemporary PCa grading system stratifies PCa into five prognostic Grade Groups based on the modified Gleason score groups: Grade Group 1 = Gleason score ≤6, Grade Group 2 = Gleason score 3 + 4 = 7, Grade Group 3 = Gleason score 4 + 3 = 7, Grade Group 4 = Gleason score 8, Grade Group 5 = Gleason scores 9 and 10 (Gordetsky and Epstein, 2016). Still, the Gleason scoring system is known to suffer from high variability across pathologists, and even within the same pathologist (Egevad et al., 2013). Moreover, it has been estimated that the Gleason grade assigned to TRUS biopsies differs compared to the whole radical prostatectomy specimen in around 50% of cases (Kvåle et al., 2009, Humphrey 2004). Although not fully understood, evidence shows that PCa exhibits intra-tumor heterogeneity (ITH) (VanderWeele et al., 2018, Berglund et al., 2018, Yadav et al., 2018), and it is not unusual to find three or four different grades present within the same prostate (multifocal) (Andreoiu et al., 2010, Wise et al., 2002).
The prostate consists of four zones: anterior fibromuscular stroma, central zone (CZ), peripheral zone (PZ), and transition zone (TZ) (McNeal 1981). Studies have characterized the prevalence of cancer arising in these zones, and the majority of prostate cancers (70%–75%) originate from the PZ with the remaining 20%–30% originating from the TZ. Cancers in the CZ are uncommon (McNeal et al., 1988, Patel et al., 2011). Generally, the TZ is regarded as the site of benign prostatic hyperplasia (BPH), i.e. low Gleason grades. While PCa arising in the TZ are frequently found incidentally in TURP specimens, the detection of TZ and anterior stromal tumors pose a clinical challenge (McNeal et al., 1988, Noguchi et al., 2000, Augustin et al., 2003). Indeed, no
report has, to our knowledge, surveyed the gene expression differences from the PZ to CZ region in PCa to better define co-existing tumor entities.
Most studies of the genome landscape of tumors focus on aberrations in epithelial cancer cells. However, emerging evidence puts emphasis on the crucial role of the tumor microenvironment (TME) and the immune system, which can both promote and suppress cancer progression (Martin et al., 2016, Hanahan and Weinberg 2011, Chen and Mellman, 2017). The current histological grading system of PCa does not capture relevant information from the TME. Ayala and coworkers (2003) were the first to develop a grading system for reactive stroma. They showed that higher levels of reactive stroma correlate with tumor grade and prognosis, which highlights an active role of tumor-stroma in PCa progression (Ayala et al., 2003). Several studies have subsequently confirmed the association of the amount of reactive stroma with PCa
aggressiveness (Yanagisawa et al., 2007, McKenney et al., 2016). However, we lack a more comprehensive transcriptome analysis of the stroma and tumor interaction that would potentially advance our understanding of tumorigenesis in PCa. A molecular map of the prostate may also provide us a useful tool to complement existing tools to subcategorize the reactive stroma in the context of different tumor classifications.
An increasing number of high-throughput RNA-seq based methods that assess molecular
conditions with spatial resolution are emerging (Junker et al., 2014, Crosetto et al., 2015, Chen et al., 2015, Moffitt et al., 2016, Ståhl et al., 2016, Lein et al., 2017, Rodriques et al., 2019, Eng et al., 2019). Here, we analyze a cross section of a whole prostate using the Spatial Transcriptomics (ST) method. In this RNA sequencing based technique, fixed tissue sections are placed onto barcoded reverse transcription primers providing both a bright field image of the tissue as well a spatially resolved transcriptome. We carried out this study to capture gradual and spatial
differences in gene expression between pathological and normal tissue compartments of a prostate. Altogether, this work represents a resource and foundation of knowledge to study PCa in the context of its tissue environment at organ-scale.
RESULTS
Spatial Transcriptomics establishes a detailed atlas of cellular heterogeneity in the prostate cancer organ
Intra-tumor heterogeneity (ITH) and its consequences over time challenges diagnosis and
effective treatment of PCa (VanderWeele et al., 2018, Berglund et al., 2018, Yadav et al., 2018). Therefore, we aimed to present an initial effort to structure and molecularly define a cross section of a prostate organ containing cancer. We took advantage of recent advancements in the field of transcriptomics that are beginning to provide spatial resolution ranging from a few cells to subcellular resolution (Vickovic et al., 2019, Rodriques et al., 2019). First, we generated high-quality Spatial Transcriptomic libraries (Salmén et al., 2018) from barcoded regions (spots) to provide an overview of a cross section of the prostate and its multifocal tumors (Figure 1). The prostate was resected by open radical prostatectomy and horizontally cut into two pieces. The upper part was used in this study and subdivided, in two steps, into 21 cubes (Figure 1, Figure S1) and then snapfrozen.
We cryosectioned 10 µm sections from the bottom of each cube (two sections per cube) for ST analysis. The sections (areas) were mounted onto spatially barcoded microarray slides. The sections were fixed, hematoxylin-eosin (HE) stained and imaged by bright field microscopy. Sections were later annotated by a trained pathologist (Figure S2). The annotations were used for comparison with results from data-driven ST gene expression analyses. The slides were then processed with the ST protocol (Berglund et al 2018) and sequencing data were analyzed using the ST pipeline (Navarro et al., 2017). Transcriptional profiling was acquired from more than 23 000 barcoded regions from the entire cross section of the prostate, and approximately 4000 UMIs and 2100 genes were identified on average per spot (Figure S3). In this study, we define low-grade cancers as Gleason score ≤ 7 and high-low-grade cancer as Gleason score ≥8.
Next, we developed a method for analyzing the barcoded gene expression data using factorized negative binomial regression (Methods). Our method yields a set of gene expression factors (GEFs), each encoding the spatial distribution of a distinct gene expression program. We first included the entire set of 21 cubes in our analysis to achieve an overview of the cross section and to compare the inferred expression factors with the pathologist's annotations. Then, we repeated
the analysis separately for each cube annotated with cancer foci or identified as cancerous based on the inferred GEFs (see next paragraph).
An initial histological overview of the prostate identified a large amount of intra-tumor
heterogeneity (ITH) within the prostate (Figure 1B, Figure S2) with up to eight cancerous tissue regions with Gleason scores ranging from (3+3) to (4+5). The amount of acute and chronic inflammation and normal stroma (smooth muscle cells) were significantly higher within the TZ, as expected. Normal epithelium was found throughout the whole prostate. High‐grade prostatic intraepithelial neoplasia (HGPIN), generally believed to be a precursor of PCa, was identified in both the TZ and PZ, but was most prevalent in the TZ. The pathologist identified one small tumor in TZ with Gleason score (3+3).
The initial data-driven analysis of the 21 cubes confirmed the pathologist annotation by comparing with the implicated function of the genes within the data-driven factors
(Supplementary Data files 1 and 2). For example, in Figure 1B and C, visualizing the inferred spatial distribution of GEF, we observe inflammation as olive-green, stroma as violet and magenta, and epithelia as light pink. The high‐grade prostatic intraepithelial neoplasia (HGPIN) was also identified (pink) as well as the small low-grade cancer in TZ (orange). In a normal prostate, neuroendocrine (NE) cells are rare and interspersed among the epithelium. In the setting of PCa, NE cells can also stimulate surrounding prostate adenocarcinoma cell growth. We
identified NE features in the prostate, most of them located in cubes within normal epithelium, and not within cancer regions (Figure 1B and C). Furthermore, they were found at multiple locations in all zones. Gene expression analysis also revealed, concordant with previous results (Henry et al., 2018), that the CHGA gene was highly upregulated in the NE features.
The first overview of the data analysis revealed seven major tumor GEFs (Figure 1B and C, Figure S4) overlapping with the eight regions annotated as cancer. Here, a tumor GEF is defined by the gene expression profile of the factor and does not imply the concept of genetic clones. As described previously (Berglund et al., 2018), there was no evident correlation between the Gleason score and the tumor GEF. In fact, areas annotated to one Gleason score were shown to contain multiple tumor GEFs. In order to study the relationships between the identified tumor GEFs, hierarchical clustering of extracted spots (at least 10 spots per clone) was performed (Figure S4). This analysis revealed two major groups—one branch comprises tumor GEFs 1, 2,
and 6, which are all spatially confined, while the other cancer clones are more spread through the prostate (Figure S4).
Next, we performed hierarchical clustering of stromal GEFs (Figure S4) comprising stroma far from tumor, stroma in proximity of tumor, and/or stroma with immune infiltrates. The analysis of all spots revealed 11 stroma GEFs, overlapping with histological annotations (Figure S4). We observed two major groups of stroma with immune cells. The first group represents an acute inflammatory response and includes markers for macrophages and neutrophils (upregulation of e.g. OLFM4, the CXCL family, MMPs, SAA1, and SAA2). The identified chemokines of the (C-X-C motif) ligand (CXCL) family are reported to be involved in various carcinogenic and anti-carcinogenic signaling pathways. In particular, CXCL1 is known to enhance angiogenesis and inhibits extracellular matrix (ECM) synthesis in PCa cells (Kuo et al., 2012). SAA, a major acute-phase protein, leads to induced expression of M2 macrophage markers including IL-10. The second stroma group involves chronic inflammation (e.g. B-cells) (Figure S4). Interestingly, the first group, acute inflammation, is spatially located within an area of high-grade cancer. Notably, as previously described (Berglund et al., 2018), pro-tumor inflammation signatures appeared to be in proximity to reactive stroma and CAFs.
Furthermore, we identified two distinct reactive stroma GEFs. Reactive stroma has previously been shown to be associated with poor outcome in clinically localized PCa (Ayala et al., 2003). The first GEF had high expression of ACKR1, RERGL, CLSTN2, DEPP1, VIM, ADGRF5, and
IGFBP7. This profile was spatially located in both TZ and PZ and, interestingly, found in both
tumorous and non-tumorous tissue. The second reactive stroma GEF had high expression of genes such as FOS, EGR1, NR4A1, ATF3, JUNB, IER2, and CTGF, many of which we have previously observed in stroma close to inflammation and cancer (Berglund et al., 2018).
Figure 1. Intra-tumoral heterogeneity and genetic characteristics of prostate cancer A) Schematic overview of data collection and analysis. B) UMAP summary of GEFs. C) UMAP visualization of the inferred spatial distribution of tumor GEFs shows cellular heterogeneity at the global level. Seven GEFs were identified.
Detailed analysis of the tumor microenvironment in individual tissue sections
To increase the detail of our data-driven annotation of the TME landscape, we repeated the factor analysis for selected tissue sections (Figure 2A and B, Supplementary Data files 3 and 4) to identify more subtle events. Here, we chose to focus and combine the analysis on the seven tissue sections (areas) harboring the largest GEFs and accompanying stroma, identified in the initial overview analysis (Figure 1). Overall, results were consistent with the previous analysis and in agreement with pathologist annotations. Yet, in some cases, we could also achieve increased granularity when matching GEFs to morphology (Figure 2B). For example, we could identify ten tumor GEFs (an increase from seven). GEFs 1--7 closely match the corresponding GEFs in the previous analysis. Tumor GEF 4 and Tumor GEF 9 represent de novo identification of areas not identified by the routine pathology. Indeed, Tumor GEF 4 is only expressed in one of the more distal parts of the prostate (area 1), and several factors are primarily present in only one tissue area (Tumor GEFs 5, 6 and 9), which is consistent with the overall heterogeneity in this PCa. Tumor GEF 10 is mainly found in two neighbouring areas (area 5 and 7), and this is also the case for Tumor GEF 2. Indeed, Tumor GEF 1 is found to be widely expressed in most cancer tissue areas although with variation in intensity, except for the region annotated as intraductal prostatic carcinoma (Tumor GEF 7). Overall, the increase of the number of cancer gene expression factors (from 7 to 10) was validated by AMACR/p63 basal immunostaining as a positive control of PCa (Figure S5).
The new analysis of seven tissue sections also facilitated a discrimination between three GEFs of stroma with immune infiltrates (hereinafter denoted tumor-immune) and three other GEFs of stroma proximal to tumor areas (hereinafter denoted tumor-stroma) (Figure 2D). Overall, the tumor-stroma patterns were observed in all tissue areas, as expected, with limited spatial overlap between the individual GEFs. The three tumor-immune GEFs had much more distinct spatial patterns with no spatial overlap between themselves. The tumor-immune GEF 3 was only found in the most distant part of the prostate along with Tumor GEF 4.
Numerous studies have shown that inflammation can promote tumor growth and metastasis by changing the biological behavior of tumor cells and activating the stromal cells in the TME, such as vascular endothelial cells, tumor activated macrophages (TAMs), and cancer associated fibroblasts (CAFs) (Coussens & Werb., 2002, Diakos et al., 2014). As mentioned, we observed three distinct tumor-immune profiles. Two serum amyloid protein coding genes (SAA1 and
SAA2) were highly upregulated in tumor-immune GEF 1 (Figure 2D). Serum amyloid A proteins
are important clinical markers for acute-phase response or inflammation and might play a role in the local inflammation of cancer (Uhlar et al., 1994, Kovacevic et al., 2008). Moreover, they are associated with tumor size and lower survival rates in several cancers (Kovacevic et al., 2008). TAMs (expression of CXCL1) and tumor associated neutrophils (TANs) (expression of LCN2) were also observed in tumor-immune GEF 1. Importantly, the spatial distribution of CXCL1 and
LCN2 overlapped with each other. CD177 was highly upregulated in tumor-immune GEF 2, and
CD177 is known to be expressed on the cell surface of neutrophils during neutrophil infiltration in tissues (acute inflammation). A number of studies have shown that neutrophils can synthesize and store molecules with known angiogenic activity, including vascular endothelial growth factor (VEGFA) (Scapini et al., 1999, Tecchio et al., 2014). This is also observed in this study, as the tumor-immune GEF 2 exhibits high levels of both CD177 and VEGFA expression. Although significant progress has been made toward understanding the role of inflammatory mediators during the innate immune response in PCa, interestingly this study shows that tumor-immune GEF 2 is located within the tumor core (holding neutrophil infiltration), whereas tumor-immune GEF 1 was spatially distributed in proximity to the tumor core (carrying TAMs and TAN). The common disorders of the prostate include benign prostatic hyperplasia (BPH), prostatitis (chronic inflammation), and PCa. The first two develop in the transition zone while the latter in the peripheral part of the gland (McNeal 1981, 1988). Studies of the transition zone have
documented presence of chronic inflammatory infiltrates in and around nodules characteristic of histologic BPH, and the presence of chronic inflammation has been implicated as a risk factor for PCa (De Marzo et al., 2007). In this prostate, chronic inflammation (tumor-immune GEF 3) was as expected more evident in area 1, which is spatially located within the transition zone (Figure 2B and 2D). Among the top genes, MMP9 expression was observed in the tumor-immune GEF 3 (Figure 2D). Matrix metalloproteinases (MMPs), are essential for both tissue remodeling and angiogenesis, which have been viewed as the key modulators of tumor invasion (Nelson et al.,
2000). MMP9 has been found to be elevated in many cancer types, such as prostate cancer (Aalinkeel et al., 2001).
The importance of cross-talk between tumor cells and neighboring stroma in the TME has been demonstrated as a key component of malignant progression and the development of therapeutic resistance (Josson et al., 2010, Palumbo et al., 2015). The reciprocal complex interactions between the stroma and epithelium result in changes in the extracellular matrix (ECM) and the tumor stroma. The TME comprises a mixture of cells such as smooth muscle, immune, and vascular cells but frequently also fibroblasts called cancer associated fibroblast (CAFs) and tumor associated macrophages (TAMs). Indeed, the TME is frequently a site of chronic
inflammation and is often known as the reactive stroma, and the presence of reactive stroma has previously been linked with more aggressive tumors, angiogenesis, ECM remodeling, and
inflammation (Tuxhorn et al., 2001, Tuxhorn et al., 2002, Palumbo et al., 2015). Prostate reactive stroma is known to have a significant decrease in smooth muscle cells and increased amount of CAFs. The the role and function of CAF in PCa is not fully understood, but the prevailing view is that this cell population has an increased expression of FAP and COL1A1, and a decreased expression of ACTA2, DES, and CNN1 (Tuxhorn et al., 2002). Furthermore, increased TGF-ß levels in stroma stimulates stromal expansion and results in facilitation of metastatic spread (Tuxhorn et al., 2002, Tu et a., 2003).
As expected, non-reactive stroma adjacent to cancer (Tumor-stroma GEF 1) was associated with increased levels of ACTA2, DES, ACTG2, and CNN1 (Figure 2D) in fairly distinct patterns. On the other hand, reactive stroma (Tumor-stroma GEF 2) with increased expression of COL1A1,
COL3A1, CTGF, and RGS5 (Figure 2D) were more widely spread and in line with reported
fibroblasts reacting to stress by increasing collagen production.
Tumor-stroma GEF 3, observed both within and close to the cancer, showed high expression levels of PAGE4, SSTR2, ATXN8OS, and DCN (Figure 2C). PAGE4 is highly expressed in proliferative inflammatory atrophy (PIA) lesions and in HGPIN, both of which are thought to be precursors of invasive PCa (Nelson et al., 2003).
We identified a distinct profile that was shared among all cancer samples (Tumor GEF 1) (Figure 2D). For instance, SLC4A4, PPFIA2, GALNT3, ARMCX5, TTC32, ENOPH1, and HIF1A were among the top genes in this GEF. Recently, gene expression profiling was used to identify
PPF1A2 as an early predictive biomarker for high grade PCa (Gleason score ≥8) (Leyten et al., 2015). Interestingly, the hypoxia marker gene HIF1A was detected in all cancer samples,
indicating that hypoxia is one of the most prominent features of the regions annotated as cancer. Hypoxia is a well-known negative prognostic factor for survival and has been associated with a shift to a more glycolytic metabolism, and increased apoptosis in cancer cells and drug resistance (Harris et al., 2002, Bristow et al., 2002).
Notably, Tumor GEF 2 was spatially localized to the right domain of the prostate and absent in the left side (Figure 2B and 2C) and overlapped with histological Gleason grading (Gs 4+4, Grade Group 4). KLK11 was shown to be elevated in this GEF. Kallikrein 11 (KLK11) is a secreted type of serine protease and has been found to be highly expressed in PCa (Stavropoulou et al., 2005). Of note, the lower grade cancers (Gleason score ≤ 7) (Tumor GEFs 3 and 8) had high expression levels of PCAT14, HMGCS2, and VGLL3, consistent with prior reports (Shukla et al., 2016) (Figure 2B and 2D). Interestingly, Tumor GEF 3 was found both within HGPIN and in a high-grade region (Gs 4+4) (area 2 and area 3) within the same area, but its expression patterns (including HMGCS2, a Myc target) indicate a HGPIN type (Figure 2B and 2D). On the other hand, we examined the expression pattern of p63 in area 2 and 3 by
immunohistochemistry. We found that both areas lacked p63 staining, which indicates presence of malignant glands (Figure 2C). Despite this, Tumor GEF 3 (high expression levels of
HMGSC2) is located at the left side of the prostate and Tumor GEF 2 at the right side (high expression levels of KLK11). Hence, GEF 3 could be linked to location rather than Gleason score (Figure 2E). Area 6 was annotated to contain one large tumor annotated with Gleason score (4+3). Within this large tumor, we found three distinct tumor GEFs (5,6 and 9) positioned
adjacent to each other and with extensive differences in gene expression (Figure 2D). TZ tumors are estimated to account for approximately 20% of PCa cases and pose a clinical challenge because of the difficulty in detecting them. We identified one tumor profile located within the TZ (Tumor GEF 4). Some of the upregulated genes in this GEF are directly related to PCa, including one known PCa diagnostic biomarker, PCA3, and one lncRNA, PCGEM1, that have been implicated in PCa progression (Lee et al., 2011, Petrovics et al., 2004) (Figure 2B and D).
The cell layers lining the border between tumor and stroma are of particular interest in terms of tumor growth as well as spread of tumor (metastasis). We found that Tumor GEF 5 appears to be
spatially located at the border in area 6. The metastatic gene PI15 was pronounced in this particular GEF (Figure 2B and 2D). Interestingly, PI15 was previously found to be a potential novel prognostic marker capable of distinguishing PCa patients with metastatic-lethal tumors from non-recurrent tumors (Zhao et al., 2017). Moreover, in the proximal area 7 we observed tumor GEF 10 with higher expression of FOS, JUNB, EGR1, SPINK1, and NR4A1. Several reports showed that c-Jun or c-Fos overexpression increased cancer invasiveness (Ouyang et al., 2008). P63 staining was positive, whereas AMACR staining was negative in this region,
suggesting that the epithelial cells are still considered normal, even though a reactive stroma and immune landscape is growing to promote invasion (Berglund et al., 2018). Taken together, different gene programs linked to invasiveness appear to occur within short distances within the tumor.
Finally, one histopathologically unfavorable prognostic marker is the presence of intraductal carcinoma of the prostate (IDC-P), which occurs in 20% to 30% of all prostate cancers. Although IDC-P is formally not included in the Gleason score, numerous studies have linked intraductal carcinoma to more aggressive disease progression (Guo and Epstein., 2006, Kimura et al., 2014, Magers et al., 2015, Wobker et al., 2016). The genetics for this aggressive subtype remain unclear. We did, however, identify the molecular profile of this subtype, including elevated expression of ACSM3, PCSK1N, PI15, TMEFF2, ARG1, and ADAM9 (Figure 2D, GEF 7).
Fig 2. Tumor-immune-microenvironment organization in prostate cancer. A) Location of cancer samples within the whole organ. B) UMAP visualization of the inferred spatial
distribution of gene expression factors for each cancer sample with corresponding annotated brightfield image of the H&E-stained tissue section. C) Heatmap of p63 staining visualized on top of the tissue (area 2 and 3). Green color indicates p63 positive staining. D) Inferred spatial distributions of all stromal, immune, and tumor factors. ND corresponds to not determined, ASAP to atypical small acinar proliferation, IDC-P to prostatic intraductal carcinoma and HGPIN to high-grade prostatic intraepithelial neoplasia. E) Spatial distribution of KLK11 and
HMGCS2.
Characterization of reactive stroma in prostate cancer
From the previous analysis, we could identify COL1A1, COL3A1, IGFBP7, and RGS5 to be expressed within the reactive stroma of high-grade PCa (Tumor-stroma GEF 2). In order to investigate gene-expression changes in reactive stroma compared to normal stroma, we performed differential gene expression analysis between the two groups in area 5 (Figure 2B), which showed a distinct separation between tumor (Gleason score 4+4) and stroma by both histology and GEFs (Methods). Histologically, normal stroma further away from the tumor exhibited notably different morphology compared to the stroma within the cancer (Figure 3). Among the 40 most significantly differentially expressed genes, we found several previously described genes for normal (e.g. ACTG2, MYH11, DES, TALGN, and MYL9) and reactive stroma (e.g. COL1A1) (Tuxhorn et al., 2002) (Figure 3C). Interestingly, several genes that we found upregulated in reactive stroma have been implicated in bone remodeling processes. For example,
COL1A1 is an osteoblastic differentiation marker (Nakajima et al., 2014). As PCa most
commonly metastasizes to bone, and Gleason 8 tumors are more likely to metastasize than lower grade tumors, overexpression of bone remodeling genes in high-grade stroma is particularly interesting. Similar findings were observed for IGFBP7, which was significantly higher expressed in reactive stroma compared to normal stroma. IGFBP7 could be a novel tumor stromal marker and promotes EMT processes in epithelial cells (Rupp et al., 2015).
Next, we investigated significantly enriched biological processes in the reactive stroma (p<0.05, fold change>1). Interestingly, processes linked to extracellular matrix organization were
significantly upregulated in reactive stroma (ECM proteoglycans). One of the main components of ECM is proteoglycans (PGs), which are expressed on cell surfaces and in the ECM. In line
with previous reports, these results suggest that changes in stromal cell-surface proteoglycans occur within the microenvironment during prostate cancer progression (Edwards et al., 2012). Hence, these molecules might be of interest for their roles in ECM structural signalling to permit cancer growth, angiogenesis, and metastasis (Theocharis et al., 2014).
Fig 3. Gene expression analysis of reactive stroma in high-grade prostate cancer. A) Spatial gene expression of selected marker genes for reactive stroma B) Enlarged H&E regions of reactive and normal stroma C) Differential gene expression heatmap comparing normal stroma to reactive stroma. Columns correspond to spots and rows correspond to the 40 genes with the highest absolute z-score (FDR < 0.01).
Extracellular matrix composition in prostate tissues
Ayala et al. (2003) showed that the degree of reactive stroma, stained with Masson's trichrome, can predict prognosis or tumor recurrence. They reported that tumors with either a low (<5%) or high (>50%) fraction of reactive stroma behaved aggressively compared to patients with a moderate (6–50%) fraction of reactive stroma fraction. Currently, no stromal features are included in the Gleason grading system, although several stromal features have been reported to predict poor prognosis (Tuxhorn et al. 2002, Ayala et al. 2003, Yanagisawa et al. 2008). To validate our suggested marker for reactive stroma, we took advantage of Masson's trichrome staining. First, we stained Masson's trichrome and used a software‐based method to obtain reproducible measures of the reactive stroma fraction in most tumor areas of the prostate, thus avoiding intra‐observer variability (Figure 4A and B). The staining enabled quantification of collagen in the ECM through Masson's trichrome stain (Figure 4A). The software-based method was able to detect regions with reactive stroma (Figure 4B).
Subsequently, the tissue was partitioned into four grades of reactive stroma (RE) (Figure 4C), and genes upregulated within each grade were calculated (Figure 4D). The highest grade of reactive stroma (RE3) displayed genes such as IGFBP7, COL1A1, COL1A2, DUSP1, MGP,
A2M, BGN, and VIM (Figure 4D). We observed that β2-Microglobulin (B2M) transcript is
highly expressed in the RE2 grade (Figure 4D). Serum B2M levels have been shown to be elevated in patients with metastatic, androgen-independent prostate cancer, and B2M expressing cells have been suggested to serve as a novel drug target for the treatment of bone metastasis (Gross et al., 2007, Huang et al., 2006). Once again, ACTG2, DES, and MYL9 were almost exclusively seen in regions with low amounts of reactive stroma (RE0-1). Notably, genes like
IGKC, IGHAI, and AZGP1 were upregulated in RE0 and RE1 (Figure 4D). Among these,
AZGP1 has, for example, been shown to exhibit tumor-suppressive properties in studies
involving breast, prostate, and other malignant tumors (Huang et al., 2013), which is in line with our results. Furthermore, loss of AR in stroma, is a significant predictor of prostate cancer recurrence, independent of Gleason grade and PSA levels (Henshall et al., 2001). Here we have shown a significant association between low AR levels in regions with high grade of reactive stroma (Figure 4E). We observed that the nuclei of the stroma in and near the tumor border have negligible AR. Additionally, stroma further away from the tumor border display a lot of AR. Interestingly, the main tumor itself had very little AR within the nuclei (Figure S6).
Fig 4. A) Masson's trichrome staining of area 5. Blue color indicates ECM proteins (collagens). B) Heatmap displaying degree of reactive stroma (dark blue = 0% reactive stroma, orange = 100% reactive stroma). C) Enlarged regions showing four different grades of reactive stroma (RE) D) Gene expression profiles for each grade of RE (0-3). E) Heatmap displaying median nuclear staining intensity of the AR in stroma cells in local stroma regions of 100X100 micrometers. AU corresponds to arbitrary units.
A higher resolution transcriptional map of prostate cancer
Next, we revisited tissue areas 1 to 7, containing the major tumor foci in the prostate, to both validate and improve our understanding of the evident complex molecular landscape in a prostate. For these purposes, we sectioned new samples, close to (but not adjacent) to the previously described and processed the new samples with five times higher density of spatially barcoded spots (Visium), enabling us to study more subtle gene expression differences in situ. In total, we obtained data from 23 282 spots from the seven tissue sections. The data was analysed using the same probabilistic approach as before and resulted in 25 spatially meaningful GEFs (Supplementary Data files 5 and 6) with an expected increase in resolution, as fewer cells per spot are captured and analysed. The top expressed genes for each factor were used to assign its cell-type identity. The spatial distribution of the factors display once again ten unique tumor GEFs, representing the different histological Gleason scores previously identified (Figure 5A and C, Supplementary Data files 5 and 6). Additionally, GEFs for normal stroma, reactive stroma, normal epithelia, and immune cells with accompanying marker genes are reproduced similar to the previous analysis (Figure 5A, B and C). One of the new genes identified was
microseminoprotein-beta (MSMB), one of the most highly secreted proteins by the prostate gland. As its expression is lost in tumorigenesis, it has been identified as a potential biomarker of PCa risk and prognosis (Dahlman et al., 2011, Grönberg et al., 2015, Sjöblom et al., 2016). We observed a decrease in expression of MSMB in tumor compared to normal tissue (Figure 5C, Figure S7).
Interestingly, a monocyte enriched GEF was readily identified with the increased resolution. Spatially, these cells were found in the periphery of the tumors. Genes such as MMP7, PIGR,
CP, SAA1, CXCL2, NNMT, NCOA7, and LCN2, characterized the monocytes (Figure 5C).
MMP7 expression has been suggested to facilitate cancer invasion and angiogenesis by degrading extracellular matrix macromolecules and confers a selective advantage for the
colonization of bones (Lynch et al., 2005). This factor is found adjacent to the tumor core (Figure 5B and C), and, generally, infiltration into the tumor indicates a poor clinical prognosis in cancers (Lin et al., 2014; Shen et al., 2014).
Another tumor immune GEF was found both nearby tumor cells and within the core of the tumor. Notably, CD74 was upregulated in that region (Figure 5B and C). Apart from its role in inflammation, many studies have shown that macrophage migration inhibitory factor (MIF), after binding CD74, has an important role in promoting tumorigenesis (Meyer-Siegler et al., 1998). This activates the PI3K/Akt pathway, and, recently, overexpression has also been shown to be associated with epithelial-to-mesenchymal transition (EMT) (Funamizu et al., 2013).
Genes in the reactive stroma were similar to the lower-resolution analysis (e.g. IGFBP7, RGS5, and VIM), but genes such as ACKR1, COMP, and BGN were now also included in the
corresponding GEF (of which all were spatially confirmed to confer the same spatial expression as shown by the GEF) (Figure 5B and C, yellow color), indicating growth of blood vessels in these cancer regions (Huss et al., 2001). The switch to an angiogenic phenotype, i.e., the ability of a tumor to recruit new vasculature, is considered to be a fundamental determinant of
neoplastic progression and a requisite for growth beyond a small nodule (Folkman 1971). Among them, BGN, a proteoglycan of the extracellular matrix, has been associated with poor prognosis in tumors harboring PTEN deletions (Jacobsen et al., 2017).
Similarly, new genes were identified for the IDC-P GEF, which now includes classic
neuroendocrine markers and genes associated with advanced stages, such as SCGN, F5 and GAL (Figure 5C, red color). As shown in Figure 5C, a Gs 3+4 cluster was identified (light purple color), expressing high levels of NPR3 and HGPD (not described previously). Strikingly, the GS 4+4 cluster in the right corner of the prostate was not shared between area 5 and 6 previously but was now identified as a unique and shared cluster with the new data (Figure 5C, green color). The reason for that could be explained by the increased resolution and sensitivity with Visium-barcoded arrays, i.e. more detected genes per spot. FABP5, PCA3, PPFIA2, TFF3, SPON2 and AMACR were discovered in this tumor type. FABP5 is not expressed in the normal prostate but becomes highly upregulated in advanced metastatic PCa (Fujita et al., 2017), in concordance with our results (Figure 5C, Figure S8).
The transition zone usually harbors benign hyperplasia and in the re-analysis of area 1 with the additional tissue section further into the ‘cube,’ we identify a clearer tumor GEF with a Gleason
score of 3+3 (Figure 5C, light green color). This GEF suggests that SRARP could be a marker of the transition tumors, as this gene has previously been described as a tumor suppressor (Naderi, 2018) (Figure 5C, green color, Figure S9). However, this finding needs to be validated further. For the analysis of tumors in the peripheral zone, we could confirm the previous findings of the poor correlation between histological scoring and GEFs. Indeed, comparing tissue area 2 and 3, although somewhat further into the tissue cubes, we have three different histological scores, (Gs 4+5, 4+4, and 3+4), and our data-driven analyses indicate a shared GEF for the same region (Figure 5C, brown color). Importantly, we also observe that the same GEF is identified in tissue area 6 at a substantial distance, suggesting a shared genetic origin. Thus, genes such as
HMGCS2, H2AFJ and VGLL3 seem to be upregulated in the tumor regions, irrespective of
cancer grade (Figure 5C, brown color). A novel spatial pattern was also observed for area 5, in the data from the higher resolution arrays. A tumor area annotated as GS 4+4 was now shown two have two tumor GEFs, with one GEF being more spread and surrounding the more nodular tumor. This particular dendritic-looking tumor cluster expressed high levels of DEPDC1 and
CRISP3 (Figure 5C, dark green color), suggesting that these genes are important during cancer
Fig 5. Spatially resolved analysis of anatomical locations of epithelial and stromal cell types in PCa. A) UMAP projections of Visium data of the gene expression factors for each cancer sample B) Regulators of the tumor immune system. C) UMAP projections overlaid on tissue sections. D) Hematoxylin & eosin stains display variations in cell shape and Gleason score throughout the prostate.
Transcriptome-wide differences in the core and the invasive tumor border
With the increased resolution, we decided to revisit and investigate area 5 having a single
Gleason grade (4+4) tumor. For this purpose, we took advantage of a recently developed R-based analysis package (STUtility; github.com/jbergenstrahle/STUtility). We performed a combined analysis of all spots from area 5 and performed non-negative matrix (NMP) factorization and then visualized the results using UMAP. The identities of the factors were annotated, as before, by cell-type markers, spatial location within the tissue, and the pathologist annotations.
In total, 17 clusters (Figures 6A and 6B) were judged to be spatially meaningful. We could confirm the previous finding of a separation between normal stroma and reactive stroma with an upregulation of MGP, COL1A1, COL3A1, COL1A2, CTGF, SPARC, and IGFBP7 in the latter. We also confirmed that the monocyte-enriched region consists of MMP7, LCN2, SAA1, and
CFB.
Interestingly, we observe a zonation of tumor cells with three distinct clusters of tumor cells located within the core and along the tumor border. We explored the genes that are significantly upregulated within each region of the cancer (Figure 6C). The tumor factor within the core of the tumor demonstrated increased expression levels of ALB, S1PR3, PLAT, CSKMT, and AMACR (cluster 10), while the dense tumor at the border had high expression of NTS, DDS, TRPM8,
MESP1, AGR2, and BANK1 (cluster 2). Interestingly, among these genes, NTS has been
suggested to contribute to the progression of several tumor types, such as PCa, and has been implicated to have an important role in the development of castration resistant prostate cancer (CRPC) (Zhu et al., 2019). Particularly, we find the less dense appearing tumor cells located outside the denser area intriguing. Our data-driven analysis suggests a local, stepwise change in gene expression from the dense core to more spread tumor cells in the same spatial context (Figure 6C), as if tumor cells were ‘shed-off’ from the tumor core. Interestingly, the more spread tumor cluster includes genes such as DEPDC1, CHN1, CRISP3, and PCAT1 (cluster 9).
signaling pathways to accelerate cell cycle progression (Huang et al., 2017). NF-kB activity not only promotes tumor cell proliferation, suppresses apoptosis, and attracts angiogenesis, it also induces epithelial mesenchymal transition, which facilitates distant metastasis (Taniguchi et al., 2018). Next, further analysis identified enriched pathways within the edge of the tumor (p<0.05, fold change>1), related with immune system (neutrophil degranulation, gene and protein
expression by JAK-STAT signaling after Interleukin-12 stimulation), stress (cellular response to stress), gene expression (transcriptional regulation by TP53), axon guidance and cell migration (Signaling by ROBO receptors, Regulation of expression of SLITs and ROBOs), signal
transduction (PTEN regulation, signaling by TGF-beta Receptor Complex, signaling by
NOTCH4, signaling by ERBB2, RHO GTPases activate PAKs, signaling by VEGF, signaling by FGFR2, MAPK6/MAPK4 signaling, signaling by WNT, mTOR signaling), metabolism
(unfolded protein response), vesicle-mediated transport (clathrin-mediated endocytosis), and DNA repair (global genome nucleotide excision repair (GG−NER)). Interestingly, the Slit/Robo pathway is a key regulator of many oncogenic pathways (e.g. HER2, EGFR, VEGFR, CD20, mTOR), regulating axon guidance, cell proliferation, cell motility, and angiogenesis, many of which could be promising anticancer drugs. Thus, the high-resolution data provides us with some insight into a gradual change in the tumor edge that nicely delineates spatially defined areas of an individual tumor.
Fig 6. Organoid tumor with invasive peripheral border and surrounding stroma. A) UMAP projections of Visium data for prostate area 5. B) UMAP projections overlaid on tissue section C) Gene expression markers of the tumor core, tumor border, and tumor cells close to the border. D) Spatial view of selected gene expression markers from (C).
DISCUSSION
PCa holds distinct tumor phenotypes ranging from indolent to aggressive disease. Due to early diagnosis by PSA screening, most men are currently identified with localized PCa, suitable for curatively planned treatment. Although PSA screening reduces the prevalence of advanced disease and mortality, trade-offs include overdiagnosis and resultant overtreatment. A better defined panorama of molecular markers reflecting different stages of tumor development might improve prognostication and reduce the risk of overtreatment of men with localized PCa. As the disease progresses, additional aspects of tumor heterogeneity evolve with multiple tumor foci and variable stromal components that also need to be considered.
The application of new methods for profiling the transcriptome of the tumor enables extraction of important cancer hallmark characteristics. These new techniques (e.g. spatial transcriptomics) are of particular relevance as they enable comprehensive characterization of the cellular
composition of the TME. We sought to establish a baseline understanding of the transcriptional landscape in a PCa from normal epithelia to high-grade cancer in a single patient. Hereby, we could both achieve a gene expression dictionary for distinct morphological entities as well as identify gradients of gene expression to distinguish, for example, reactive stroma from normal stroma and to understand the tumor edge effects in contrast to the tumor core. Notably, we could identify gene expression profiles associated with aggressive PCa, marked by expression in particular of FABP5 and AMACR. Moreover, we could translate the gene expression signatures to explore the underlying biological processes in the context of a spatial landscape.
The Gleason grading system, developed in the 1960s, remains the strongest prognostic factor for PCa. Improvements to the Gleason grading system have aimed to provide more accurate
diagnostics, especially for the Gleason score 7 group (scores 3+4 and 4+3), which exhibits great variability in the cancer outcome (Amin et al. 2011, Pierorazio et al., 2013, Epstein et al., 2016). However, despite the upgraded classification system, it relies on a visual assessment, which remains a limitation. Combining different molecular markers into a panel has been suggested to enhance diagnostic performance, and efforts are ongoing to validate and use such strategies at a large scale (Grönberg et al., 2015). In this study, we highlighted the challenges associated with the diagnosis of a multifocal tumor unless the individual tumor foci are not investigated
differences that potentially could be important for diagnosis. In fact, we demonstrated that a single patient had more than 10 different cancer GEFs. Some of these overlapped with multiple and different Gleason scores assessed by trained pathologists, and, in addition, some of these spatially distinct areas were not identified in the first round of histological examination, thus underpinning the need for new strategies. Indeed, our data implies the need to perform an even broader spatial analysis before developing the next generation of stratified biomarker repertoires. Modelling of single tissue sections and their spatial expression at various stages will be
important to achieve a temporal understanding of how the tumors evolve over time.
Growing evidence indicates that we need and should consider the stroma in a more diagnostic manner, yet not the field is fully explored. The paradigm of targeting the tumor cell is now being expanded beyond classic chemotherapy, with the recognition that cancer cells exist in complex microenvironments that offer therapeutic targets, including the immune system and the
supporting stromal cells.
This study demonstrates that we can use data-driven approaches to identify areas of reactive stroma in proximity of the tumor, which we confirm by Masson's trichrome staining. The identified GEFs may provide new therapeutic avenues to treat the tumor by addressing its microenvironment. Of particular interest would be to monitor the immune cells and their distribution. For example, we could in this study identify a monocyte-enriched cluster, with the higher resolution arrays, demonstrating that we will need larger data sets to gain further insight of the role of particular immune infiltrates.
It is clear that publicly available RNA sequencing data from prostate tumors do not fully capture the full extent of the heterogeneity that may be present in a PCa tumor. Thus, there is a need to revisit existing tumor banks with spatial technology to achieve a draft of the molecular anatomy of PCa that considers both tumor and stroma. We have used principles of AI to identify shared gene expression profiles from the spatial RNA sequencing data, and the same approach can be used to expand the dictionary. In the future, super-resolution ST could provide us with an even more granular view of disease progression in PCa (Bergenstråhle et al., 2020). Moreover, our transcriptional atlas of the prostate can provide a resource for the single cell community to map their annotated cell clusters (i.e. cell types) to achieve further understanding of the spatial orientation of cell types in normal as well pathological conditions (Andersson et al., 2019). At
present high-quality single cell resources for prostate and PCa is limited but are likely to soon evolve from Human Cell Atlas initiative (www.humancellatlas.com).
In this study, we have presented an unprecedented view of cancer in a cross section of the prostate. We observe significant molecular heterogeneity that goes beyond the histology. We employ data-driven approaches to find gene expression commonalities, and the spatial
coordinates give us the chance to structure the data into GEFs representing different phases of normal and pathological processes, albeit as snapshots in time. The role of many of the genes in the development or progression of PCa remains to be determined. Nevertheless, this study suggests that the GEFs may represent molecular signatures of biological processes relevant to tumor progression that are not appreciable by pathological analysis. The richness of the data is immense and can initially be used to gain a biological understanding of the tumor and its microenvironment, but it also gives us new molecular tools to better diagnose PCa in a more unbiased manner, which is greatly needed.
METHODS
Prostate cancer material
A prostate was resected by open radical prostatectomy. The prostate was cut in two halves by a horizontal cut, the upper part (closest to the patient's head) was used and cut on a 5 mm mold to obtain a 5 mm high cylinder. Next, stripes were cut out from the cylinder, and each stripe was cut into smaller cubes (21 in total). All sections used for ST experiments were taken from the bottom of the cubes from the cylinder.
We obtained 21 fresh-frozen tissue cubes, representing the whole organ. Collected samples were fresh-frozen in liquid nitrogen and stored at −80 °C until embedding for cryosectioning. An expert pathologist evaluated all hematoxylin and eosin-stained sections covering all cubes.
Spatial Transcriptomics
ST was developed based on the capture of mRNA from a tissue on glass slides coated with oligonucleotides associated with spatial barcodes followed by library preparation and
sequencing. ST combines high-throughput accurate quantitative data, spatial information, and high-resolution imaging. More specifically, every capture probe has a molecular id (barcode) that is unique to the spot. This barcode is present in the reads generated by the sequencer and can therefore be used to place the reads (gene expression values) in a spatial context. The ST slides are coated with oligo-dT probes, whereas library preparation subarrays contain 1,007 spots each, with ~200 million oligo-dT probes characterized by a spot-specific barcode (Ståhl et al., 2016, Salmén et al., 2018).
We followed a protocol described in Ståhl et al. and Salmén et al. The protocol starts with cryosectioning followed by tissue fixation to preserve tissue morphology and RNA quality. The fixation step on the slide surface was performed using 4% formaldehyde for 10 min at room temperature. The sections were stained using standard H&E staining, as previously described (Ståhl et al., 2016). Bright-field imaging is subsequently performed to identify and record the tissue histology. Tissue sections were permeabilized using exonuclease I buffer for 30 min at 37 °C and 0.1X pepsin (pH 1) for 10 min at 37 °C. Polyadenylated transcripts are captured by the surface probes, which function as primers in the overnight reverse transcription reaction (42 °C). The release of the cDNA–mRNA hybrids and subsequent steps are not included below, as they are described in the accompanying protocol by Ståhl et al and Salmén et al. The material was
processed into libraries as described in Jemt et al. and sequenced on an Illumina Novaseq using paired-end 300-bp reads. FASTQ files were processed using the ST Pipeline v.1.5.1 software (Navarro et al., 2017). R1 contained only a spatial barcode followed by a UMI. R2 transcripts were mapped with STAR (Dobin et al., 2013) to the GRCh38.79 human reference genomes. Mapped reads were counted using the HTseq count tool (Anders et al., 2015). Spatial barcodes were demultiplexed using an implementation of TagGD (Costea et al., 2013). UMI filtering is carried out to remove duplicated reads. The output is a matrix with gene counts for each spatial barcode.
Factorized negative binomial regression
To deconvolve the expression data into GEFs, we developed a method for factorized negative binomial regression. Briefly, we assume that the observed number of transcripts for a given spot 𝑠 and gene 𝑔 is a sum of hidden counts 𝑦 from 𝑇 GEFs,
We assume that the hidden count 𝑦%& follows a negative binomial distribution with a type-specific rate. Model parameters are found by maximum a posteriori estimation and used to infer 𝑦%&.
To visualize the GEFs spatially, we first reduce the relative hidden counts in each spot to three dimensions using UMAP. The resulting data is scaled into the unit cube using a channel-wise affine transformation and used as RGB color coordinates. See Supplementary Note 1 for a complete description of our method and Code Availability for the code that we used to run the analyses.
In all analyses, we factorized the data into 𝑇 = 25 GEFs and ran the optimization for 5000 iterations. Spots were annotated based on their section to control for sample-wise batch effects.
Immunohistochemistry
Tissue sections of 10 micrometer on superfrost slides, stored in −80 °C were thawed in RT to be fixed with 3% freshly made paraformaldehyde in TBS for 10 min in RT. Tissues were then permeabilized for 10 min in TBS+ 0.1%Triton-X100, rinsed three times in TBS for 5 min, and rinsed and blocked with 2% bovine serum albumin in TBS for 2 h. Incubation with the primary antibodies (AMACR: (13H4), Life technologies, 1:150, P63: (4A4), Abcam, 1:150, AR: (N-20), SCBT, 1:400, KI67: (MIB-1), Dako, 1:300) was done overnight at 4 °C. Thereafter, tissue were rinsed 3 × 5 min with TBS before incubation with the secondary antibodies donkey anti-mouse immunoglobulin G (IgG)-AlexaFluor 568 (1:500 Molecular Probes) and donkey anti-rabbit IgGAlexaFluor 647 (1:500 Molecular Probes) for 1 h at RT in darkness. DNA was
counterstained with DAPI (Molecular Probes) and slides were mounted with Prolong Gold (Molecular probes). Fluorescence images were obtained with a Zeiss LSM 780 inverted confocal microscope, using a Plan Apochromat 20 × /NA 0.7 objective. Tiled images were acquired from optical sections of 5 micrometer.
Hierarchical clustering
For hierarchical clustering, we used the ward.D aggregation method of the hclust function from the hclust package in R, with the Euclidean distance as the dissimilarity measure.
Differential gene-expression analysis of the stroma
We extracted 13 spots from reactive stroma located within the tumor and 12 spots covering normal stroma (located outside the tumor area). Only genes with more than one count in both groups were considered for further analysis. The samples were transformed using the rlog function from the DESeq2 package in R (Love et al., 2014). Genes with significant differential expression (FDR < 0.01) were considered, and z-scores were calculated for each gene and group. Final results were plotted using the heatmap.2 function from the gplots package in R.
Masson's trichrome staining
The staining was performed with the Trichrome Stain Kit from Abcam. The protocol was modified to suit snap frozen tissue. The modification mainly followed the protocol in SOP (ID) Number MDC1A_M.1.2.003. Frozen tissue sections were thawed in RT, fixated with 4% freshly made paraformaldehyde in TBS for 1h in RT and subsequently re-fixed in Bouin's solution overnight. After one min rinse in tap water, followed by a two second rinse in ddH2O, the tissue
was incubated at 5 min in freshly made Weigert's solution. Tissues were then rinsed in 40 °C tap water for 7 min followed by a one min rinse in ddH2O to be incubated 5 min in Biebrich Scarlet acid fuchsin solution. After a 3 × 1 min rinse in ddH2O the tissue were differentiated in
phosphomolybdic/phosphotungstic acid solution for 15 min followed by a 7 min incubation with aniline blue stain and finally rinsed 2X1 min with ddH2O. Slides were then dehydrated in 70%, 90% and 100% etOH for 3 min each followed by 5 min in xylol before draining and mounting with Depex.
Differential gene-expression analysis of the different RE grades
We extracted a group of spots (at least 6 spots) from each reactive stroma grade. Only genes with more than one count in both groups were considered for further analysis. The samples were transformed using the rlog function from the DESeq2 package in R (Love et al., 2014). Genes with significant differential expression (FDR < 0.01) were considered, and z-scores were
calculated for each gene and group. Final results were plotted using the heatmap.2 function from the gplots package in R.
Clustering and annotation of spots in Fig. 6
The analysis was performed in R version 3.6.2. Count matrices were loaded into R and explored using the STUtility package (github.com/jbergenstrahle/STUtility) and the Seurat package (Butler et al., 2018). Spots with less than 500 genes and genes with less than a total read count below 100 were filtered away. Mitochondrial genes were removed. Normalization and batch correction were achieved using the variance-stabilizing transformation implemented in the SCTransform function from the Seurat package, with the tissue section specified as the batch variable in the “vars.to.regress” option. Subsequent dimensionality reduction was done using non-negative matrix factorization (NMF), and the resulting factors were embedded in a two-dimensional space using UMAP for visualizing spot similarities. Clustering of the data was performed using the FindNeighbors and FindClusters function from the Seurat package with the resolution parameter set to 1.5. Differentially expressed genes for each cluster were determined by FindMarkers with parameters min.pct=0.25, logfc.threshold=0.25, and test.use=”wilcox”. Differential gene expression analysis between normal stroma and stroma close to the border of the tumor was done using the DESeq2 package in R. Genes that met the above criteria with adjusted P-value < 0.05 (Benjamini-Hochberg correction) and fold change > 1 were used for
pathway analysis. Pathway analysis of the differentially expressed genes was accomplished using the ReactomePA package in R.
ACKNOWLEDGMENTS
This study was supported by AstraZeneca and Science for Life Laboratory. We would like to thank the National Genomics Infrastructure (NGI), Sweden for providing infrastructure support. We thank Alma Andersson for helpful discussions.
AUTHOR CONTRIBUTIONS
E.B, N.S, M.M and R.M performed the experiments. E.B, N.S, M.M, J.B, L.L analyzed the data. E.B, N.S, M.M, L.B and J.M designed the research. T.H and J.L supervised the research. F.T provided the prostate material and A.T (expert pathologist) annotated all sections. E.B, N.S, M.M, L.B and J.L wrote the manuscript.
DECLARATION OF INTERESTS
REFERENCES
Aalinkeel R, Nair BB, Reynolds JL, et al. Overexpression of MMP-9 contributes to invasiveness of prostate cancer cell line LNCaP. Immunol Invest. 2011;40(5):447–464. doi:10.3109/08820139.2011.557795
Amin A, Partin A, Epstein JI. Gleason score 7 prostate cancer on needle biopsy: relation of primary pattern 3 or 4 to pathological stage and progression after radical prostatectomy. J Urol. 2011;186(4):1286–1290.
doi:10.1016/j.juro.2011.05.075
Alma Andersson, Joseph Bergenstråhle, Michaela Asp, Ludvig Bergenstråhle, Aleksandra Jurek, José Fernández Navarro, Joakim Lundeberg. Spatial mapping of cell types by integration of transcriptomics data. bioRxiv 2019.12.13.874495; doi: https://doi.org/10.1101/2019.12.13.874495
Andreoiu M., Cheng L. Multifocal prostate cancer: Biologic, prognostic, and therapeutic implications. Hum. Pathol. 2010;41:781–793. doi: 10.1016/j.humpath.2010.02.011.
Augustin H, Erbersdobler A, Graefen M, et al. Differences in biopsy features between prostate cancers located in the transition and peripheral zone. BJU Int. 2003;91(6):477–481. doi:10.1046/j.1464-410x.2003.04140.x
Ayala G., Tuxhorn J. A., Wheeler T. M., Frolov A., Scardino P. T., Ohori M., Wheeler M., Spitler J., and Rowley D. R. (2003) Reactive stroma as a predictor of biochemical-free recurrence in prostate cancer. Clin. Cancer Res. 9, 4792–4801
Bergenstråhle L, He B, Bergenstråhle J, Andersson A, Lundeberg J, Zou J, et al. Super-resolved spatial transcriptomics by deep data fusion. BioRxiv 2020:2020.02.28.963413. doi:10.1101/2020.02.28.963413. Bristow RG, Hill RP. Hypoxia and metabolism. Hypoxia, DNA repair and genetic instability. Nat Rev Cancer. 2008;8:180–192. doi: 10.1038/nrc2344.
Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017;541:321–330. doi: 10.1038/nature21349.
Chen KH, Boettiger AN, Moffitt JR, Wang S, Zhuang X. RNA imaging. Spatially resolved, highly multiplexed RNA profiling in single cells. Science. 2015;348(6233):aaa6090. doi:10.1126/science.aaa6090
Chi KN, Bjartell A, Dearnaley D, Saad F, Schroder FH, Sternberg C, Tombal B, Visakorpi T. Castration-resistant prostate cancer: from new pathophysiology to new treatment targets. Eur Urol. 2009;56(4):594–605.
Coussens L. M., & Werb Z. (2002). Inflammation and cancer. Nature, 420, 860–867. 10.1038/nature0132
Crosetto N., Bienko M., van Oudenaarden A., Spatially resolved transcriptomics and beyond. Nat. Rev. Genet. 16, 57–66 (2015).
Dahlman A, Rexhepaj E, Brennan DJ, et al. Evaluation of the prognostic significance of MSMB and CRISP3 in prostate cancer using automated image analysis. Mod Pathol. 2011;24(5):708–719.
doi:10.1038/modpathol.2010.238
De Marzo AM, Platz EA, Sutcliffe S, Xu J, Grönberg H, Drake CG, Nakai Y, Isaacs WB, Nelson WG. Inflammation in prostate carcinogenesis. Nat Rev Cancer. 2007;7:256–269.
Diakos C. I., Charles K. A., McMillan D. C., & Clarke S. J. (2014). Cancer‐related inflammation and treatment effectiveness. The Lancet Oncology, 15, e493–e503. 10.1016/S1470-2045(14)70263-3
Edwards IJ. Proteoglycans in prostate cancer. Nat Rev Urol. 2012;9(4):196–206. Published 2012 Feb 21. doi:10.1038/nrurol.2012.19
Egevad L, Ahmad AS, Algaba F, et al. Standardization of Gleason grading among 337 European pathologists. Histopathology. 2013;62:247–256. doi: 10.1111/his.12008.
Epstein JI, Zelefsky MJ, Sjoberg DD, et al. A Contemporary Prostate Cancer Grading System: A Validated Alternative to the Gleason Score. Eur Urol. 2016;69(3):428–435. doi:10.1016/j.eururo.2015.06.046 Folkman J. Tumor angiogenesis theraperutic implications. N Engl J Med. 1971;285:1182–6.
Fujita K, Kume H, Matsuzaki K, et al. Proteomic analysis of urinary extracellular vesicles from high Gleason score prostate cancer. Sci Rep. 2017;7:42961.
Funamizu N, Hu C, Lacy C, et al. Macrophage migration inhibitory factor induces epithelial to mesenchymal transition, enhances tumor aggressiveness and predicts clinical outcome in resected pancreatic ductal adenocarcinoma. Int J Cancer. 2013;132(4):785–794. doi:10.1002/ijc.27736
Gleason DF. Classification of prostatic carcinomas. Cancer Chemother Rep. 1966;50:125–128.
Gordetsky J, Epstein J. Grading of prostatic adenocarcinoma: current state and prognostic implications. Diagn Pathol. (2016) 11:25. 10.1186/s13000-016-0478-2
Gross M, Top I, Laux I, et al. Beta-2-microglobulin is an androgen-regulated secreted protein elevated in serum of patients with advanced prostate cancer. Clin Cancer Res. 2007;13(7):1979–1986. doi:10.1158/1078-0432.CCR-06-1156
Grönberg H, Adolfsson J, Aly M, et al. Prostate cancer screening in men aged 50-69 years (STHLM3): a prospective population-based diagnostic study. Lancet Oncol. 2015;16(16):1667–1676. doi:10.1016/S1470-2045(15)00361-7. Grönberg H. Prostate cancer epidemiology. Lancet. 2003;361(9360):859‐864. doi:10.1016/S0140-6736(03)12713-4 Guo CC, Epstein JI. Intraductal carcinoma of the prostate on needle biopsy: histologic features and clinical
significance. Mod Pathol. 2006;19(12):1528–1535.
Han S, Cao C, Tang T, Lu C, Xu J, Wang S, Xue L, Zhang X, Li M. ROBO3 promotes growth and metastasis of pancreatic carcinoma. Cancer letters. 2015;366:61–70.
Harris AL. Hypoxia - a key regulatory factor in tumour growth. Nat Rev Cancer. 2002;2:38–47. doi: 10.1038/nrc704 Henry GH, Malewska A, Joseph DB, et al. A Cellular Anatomy of the Normal Adult Human Prostate and Prostatic Urethra. Cell Rep. 2018;25(12):3530–3542.e5. doi:10.1016/j.celrep.2018.11.086
Henshall SM, Quinn DI, Lee CS, et al. Altered expression of androgen receptor in the malignant epithelium and adjacent stroma is associated with early relapse in prostate cancer. Cancer Res. 2001;61(2):423–427
Huang L, Chen K, Cai ZP, et al. DEPDC1 promotes cell proliferation and tumor growth via activation of E2F signaling in prostate cancer. Biochem Biophys Res Commun. 2017;490(3):707–712. doi:10.1016/j.bbrc.2017.06.105 Hill R, Song Y, Cardiff RD, Van Dyke T. Selective evolution of stromal mesenchyme with p53 loss in response to epithelial tumorigenesis. Cell. 2005;123: 1001–11. 10.1016/j.cell.2005.09.030
Huang WC, Wu D, Xie Z, et al. beta2-microglobulin is a signaling and growth-promoting factor for human prostate cancer bone metastasis. Cancer Res. 2006;66(18):9108–9116. doi:10.1158/0008-5472.CAN-06-1996
Huang CY, Zhao JJ, Lv L, et al. Decreased expression of AZGP1 is associated with poor prognosis in primary gastric cancer. PLoS One. 2013;8(7):e69155. Published 2013 Jul 23. doi:10.1371/journal.pone.0069155
Humphrey PA. Gleason grading and prognostic factors in carcinoma of the prostate. Mod Pathol. 2004;17(3):292– 306.
Huss WJ, Hanrahan CF, Barrios RJ, Simons JW, Greenberg NM. Angiogenesis and prostate cancer: identification of a molecular progression switch. Cancer Res. 2001;61(6):2736–2743.
Jacobsen F, Kraft J, Schroeder C, et al. Up-regulation of Biglycan is Associated with Poor Prognosis and PTEN Deletion in Patients with Prostate Cancer. Neoplasia. 2017;19(9):707–715. doi:10.1016/j.neo.2017.06.003 Josson S., Matsuoka Y., Chung L. W., Zhau H. E., Wang R. (2010). Tumor-stroma co-evolution in prostate cancer progression and metastasis. Semin. Cell Dev. Biol. 21, 26–32. 10.1016/j.semcdb.2009.11.016
Kimura K, Tsuzuki T, Kato M, et al. Prognostic value of intraductal carcinoma of the prostate in radical prostatectomy specimens. Prostate. 2014;74:680–7. doi: 10.1002/pros.22786.
Kovacevic A, et al. Expression of serum amyloid A transcripts in human bone tissues, differentiated osteoblast-like stem cells and human osteosarcoma cell lines. J Cell Biochem. 2008;103(3):994–1004.
Kulkarni P., Jolly M.K., Jia D., Mooney S.M., Bhargava A., Kagohara L.T., Chen Y., Hao P., He Y., Veltri R.W., et al. Phosphorylation-induced conformational dynamics in an intrinsically disordered protein and potential role in phenotypic heterogeneity. Proc. Natl. Acad. Sci. USA. 2017;114:E2644–E2653. doi: 10.1073/pnas.1700082114. Kuo PL, Shen KH, Hung SH, Hsu YL. CXCL1/GROα increases cell migration and invasion of prostate cancer by decreasing fibulin-1 expression through NF-κB/HDAC1 epigenetic regulation. Carcinogenesis. 2012;33:2477–2487. doi: 10.1093/carcin/bgs299
Kvåle R, Møller B, Wahlqvist R, et al. Concordance between Gleason scores of needle biopsies and radical prostatectomy specimens: a population-based study. BJU Int. 2009;103(12):1647–1654.
Lee C. H., Akin-Olugbade O., Kirschenbaum A. Overview of prostate anatomy, histology, and pathology.
Endocrinology and Metabolism Clinics of North America. 2011;40(3):565–575. doi: 10.1016/j.ecl.2011.05.012.
Lee GL, Dobi A, Srivastava S. Prostate cancer: diagnostic performance of the PCA3 urine test. Nat Rev Urol. 2011;8(3):123–124. doi:10.1038/nrurol.2011.10
Lein E., Borm L. E., Linnarsson S., The promise of spatial transcriptomics for neuroscience in the era of molecular cell typing. Science 358, 64–69 (2017).
Leng X, Wu Y, Arlinghaus RB. Relationships of lipocalin 2 with breast tumorigenesis and metastasis. J Cell
Physiol. 2011;226(2):309–314. doi:10.1002/jcp.22403
Leyten GHJM, Hessels D, Smit FP, Jannink SA, de Jong H, Melchers WJG, et al. Identification of a candidate gene panel for the early diagnosis of prostate Cancer. Clin Cancer Res. 2015;21:3061–3070. doi:
10.1158/1078-0432.CCR-14-3334.
Lin GN, Peng JW, Xiao JJ, Liu DY, Xia ZJ. Prognostic impact of circulating monocytes and lymphocyte-to-monocyte ratio on previously untreated metastatic non-small cell lung cancer patients receiving platinum-based doublet. Med Oncol. 2014;31(7):70. doi:10.1007/s12032-014-0070-0
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
Genome Biol. 2014;15(12):550. doi:10.1186/s13059-014-0550-8
Lynch, CC, Hikosaka, A, Acuff, HB, Martin, MD, Kawai, N, Singh, RK, Vargo-Gogola, TC, Begtrup, JL, Peterson, TE, Fingleton, B, Shirai, T, Matrisian, LM & Futakuchi, M 2005, 'MMP-7 promotes prostate cancer-induced osteolysis via the solubilization of RANKL', Cancer Cell, vol. 7, no. 5, pp. 485-496.
