{'fraction_non_alphanumeric': 0.036418816388467376, 'fraction_numerical': 0.013895087625338452, 'mean_word_length': 5.419022154316272, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neurogenesis and angiogenesis are associated features of aggressive breast cancer

The tumor microenvironment (TME) is important for cancer growth and progress. Whereas angiogenesis is an established hallmark of cancer, the field of cancer neuroscience is now emerging. Here, we asked whether neurogenesis and angiogenesis are associated features in breast cancer and related to aggressive behavior. Single cell-based spatial analysis by imaging mass cytometry (IMC) indicated significant colocalization of neural and vascular structures, suggesting the presence of neuro-vascular niches within the tumor tissue. By protein and gene expression data from human breast cancer cohorts, we found that neurogenesis and angiogenesis, together with a consolidated neuro-angiogenic signature score, were linked to high-grade breast cancer characteristics and reduced survival, also within the luminal tumor subgroup. Cases with high neuro-angiogenic score were split in two subgroups, one hormone receptor negative and one receptor positive. Our findings might be relevant for improved patient stratification and further exploration of novel treatment strategies.


Introduction
Tumor-stroma communication is important for cancer progression and is suggested as a therapeutic target. Programs of inflammation, angiogenesis, and fibroblast activation are among the tumor microenvironment (TME) characteristics of aggressive cancer (1,2). Recently, evidence has indicated links between tumor cells, nerve elements and progressive disease (3,4), and experimental studies suggest that cancer cells may secrete certain factors that simultaneously promote the recruitment of both sprouting axons and endothelial cells (5)(6)(7)(8)(9). Studies of neural tissue in breast cancer stroma are limited (10)(11)(12)(13)(14)(15)(16), and the potential coexistence and coregulation of microaxons and microvessels have, to our knowledge, not been examined in human tumor tissues.
Here, we investigated whether neurogenesis (i.e., axonogenesis) and angiogenesis are associated features in breast cancer TME, and the relationship between these structures and aggressive cancer phenotypes were examined. Using single cell-based spatial analysis by imaging mass cytometry (IMC), together with protein and gene expression data from human breast cancer cohorts, our tissue-based findings indicate that neurogenesis and angiogenesis are coexisting and appear to have an impact on disease progression. Our findings provide a basis for better tumor classification and improved patient management.

Markers of neurogenesis and angiogenesis associate with aggressive breast cancer
Our purpose was to study the presence and significance of neurogenesis (i.e., presence of nerve fibers) and angiogenesis in human breast cancer, their association, and how they are linked to clinico-pathologic phenotypes and patient prognosis. We 4 also wanted to examine the spatial relationship between these components in the TME. We first applied immunohistochemistry (IHC) on tissues from 534 primary breast cancers and 95 matched lymph node metastases (Bergen Breast Cancer Cohort).
No associations were observed between NBD and primary tumor features.

Markers of neurogenesis and angiogenesis are associated traits in breast cancer tissues
By case-based correlation analysis, neurogenesis (by MAD count) and angiogenesis  (Figure 1g-k). In contrast, no associations were observed between nerve bundles and vascular markers. Taken together, our IHC studies revealed positive correlations between neurogenesis and angiogenesis in primary breast cancer tissues.

Spatial analysis of neural and vascular structures indicates co-localization in breast cancer tissues
To investigate the single-cell spatial distribution and potential co-localization of neural and vascular structures in breast cancer, we designed an imaging mass cytometry 6 (IMC) panel of 38 markers; 12 of these markers were included in the current analyses (Supplementary Table 5) to profile 10 samples from luminal-like (n=5) and basal-like (n=5) tumors. IMC uses isotope-labelled antibodies combined with laser ablation and mass-spectrometry to produce high dimensional images which are further segmented into single-cell information to investigate co-expression and spatial distribution of various markers. Altogether, we identified 46,574 single-cells or equally denoted cell units (i.e., 25,532 from basal-like and 21,042 from luminal-like tissues), quantifying the marker expression together with spatial information of each cell unit. An unsupervised clustering approach was used to identify cell phenotypic clusters of different epithelial phenotypes, T-cells, B-cells, macrophages, endothelial cells, and stromal cells (17).

Vascular proliferation was determined by dividing endothelial cells on the basis of Ki67
and CD31 co-expression, using density estimates from the flowDensity software (18), resulting in the following four subsets (denoted as vascular phenotypes): Ki67low/CD31-high, Ki67-high/CD31-high, Ki67-high/CD31-low and Ki67-low/CD31-low (Supplementary Figure 3a). Quantification of cell units in each subset (Figure 2a) indicated that higher number of proliferating vessels was present in basal-like compared to luminal-like tumors (69.03% vs. 30.97%). Our analysis showed that putative pre-existing and not proliferating vessels were present in both subgroups, with luminal-like cancers having a higher relative percentage.
To elucidate whether markers of neural fibers and vasculature are co-expressed within the same cell units (or within overlapping and densely packed units), we used density estimation of the Neurofilament marker that resulted in groups of low and high expression (Supplementary Figure 3b). Our analyses revealed that cell units of higher vascularization (CD31-high) exhibited higher levels of neural fiber elements based on Neurofilament expression (P<0.001; Figure 2b-c). Thus, using median expression values of CD31 and Neurofilament across the studied cell units, our analysis indicated consolidated neurogenic and angiogenic features, i.e. the existence of a "neuro-angiogenic profile" (Figure 2d). We then asked whether the consolidated neuro-angiogenic profile could be used to discriminate between basal-like and luminallike cases at the sample level. To do so, we used as input the consolidated neuroangiogenic profile and projected basal-like and luminal-like samples into 2D space, using Multidimensional Scaling (MDS). We found that the basal-like samples were projected in closer proximity to each other, and they were clearly separated from the luminal-like samples (Figure 2e).
The identified neuro-angiogenic profiles or units were also detected visually, when cell units were mapped to the 2D space using t-distributed stochastic neighbor embedding (TSNE) (Figure 3a). Further, we used IMC-based images to detect areas with spatial co-localization of neural and vascular structures. For each endothelial cell unit, we estimated the number of neighboring cells using a k-nearest approach, and we visualized the average distance from all cell units with high Neurofilament marker expression. Our results indicated that proliferating vessels (Ki67-high/CD31-high) were on average in closer proximity to cell units of high neural participation (by Neurofilament expression; Kolmogorov-Smirnov statistic; Figure 3b), when compared with the rest. Our results were validated using pseudo-colored IMC images with representative markers indicating proximity between neural and vascular structures (Figure 3c-d).
Taken together, our single-cell spatial analysis revealed that markers reflecting neurogenesis (Neurofilament) and angiogenesis (CD31), are co-expressed and/or colocalized, supporting relationships which might be functionally important. Thus, our 8 single cell-based IMC data are in concordance with observations from our case-based IHC analyses.

Proteomic analysis indicates differential expression of neurogenic and angiogenic proteins in breast cancer tissues
We then used mass spectrometry profiling of 24 laser microdissected tumor epithelial samples from basal-like (n=12) and luminal-like (n=12) breast cancer tissues to understand whether neurogenic and angiogenic proteins are differentially expressed in luminal and basal cell types, Gene Ontology (GO) annotations were used to identify gene sets related to axonogenic (GO:0007409), neurogenic (GO:0022008) and angiogenic (GO:0001525) processes. Due to the high degree of overlap between axonogenic and neurogenic annotations from GO, these two sets were merged and are referred to as "neurogenic" in the following. After duplicate removal, we compiled a set of 473 genes characterizing axon morphogenesis, axon growth, neural process generation, nervous system cell generation, neural cell differentiation, and another set of 494 genes characterizing blood vessel formation from pre-existing vessels.

Gene expression analysis support neurogenesis and angiogenesis as features of aggressive breast cancer
Following up on our results from IHC, IMC and proteomics analyses, we wanted to elucidate how gene expression patterns representing sprouting axons and vascular proliferation relate to breast cancer phenotypes, by exploring gene expression programs reflecting tumor progression. By analyses of global gene expression data from human breast cancer in METABRIC (n=1710) and TCGA (n=505) databases, and breast cancer cell line data from CCLE (n= 59) and in-house cell lines (n=12), we mapped gene expression signatures reflecting stroke-associated sprouting axons (19) and vascular proliferation (20), and we then estimated the sprouting axons score (SAS) and vascular proliferation score (VPS; See Methods).
High SAS and VPS were both associated with ER and PR negative tumors (P<0.001), and their distribution across breast cancer subtypes revealed consistently higher scores in basal-like and HER2 enriched categories (Figure 5a- Next, we derived a consolidated score reflecting both sprouting axons and vascular proliferation, termed a neuro-angiogenic score, NAS (see Methods). We found a distinct separation between molecular subtypes, as also found in breast cancer cell lines; the highest scores were seen in basal-like and HER2 enriched subtypes (

High neuro-angiogenic score associates with signatures reflecting hypoxia responses.
To further pursue our GSEA-based results, we found strong positive correlations between NAS and independent hypoxia signature scores (28,29) Phenotypic heterogeneity among tumors with increased neurogenesis and angiogenesis scores. When clustering dichotomized signature scores reflecting EMT, stemness, and hypoxia, along with phenotypic clinico-pathologic markers in accordance with high or low NAS, phenotypic subgroups among tumors were evident (four clusters, C1-C4; Figure 6a-b). One subgroup of NAS-high tumors (C4) was ER negative and basal-like with high scores for EMT, TGF-β activation, stemness and hypoxia; the other subset of NAS-high cases (C3) was ER positive and non-basal but still associated with high EMT and stemness scores. No survival difference was observed between the two NAS-high subclusters (C3; C4).
The survival pattern for these clusters demonstrated best prognosis in cases of C1 and poorer survival in C2-C4 (Figure 6c), with no differences between C2-C4 (P>0.15). To summarize, we demonstrate phenotypic heterogeneity in the subset of tumors with NAS-high tumors, where ER status and hypoxia split the high-score cases. 13 As we demonstrated correlations between markers of sprouting axons and angiogenesis in breast cancer tissues and cell lines (protein and mRNA levels), and associations of these with aggressive breast cancer phenotypes, we next assessed how measures of neurogenesis and angiogenesis relate to breast cancer survival.

The neuro-angiogenic tumor profile associates with reduced patient survival
In our IHC-based tissue studies using Neurofilament positive microaxon densities (MAD) and nerve bundles (NBD) in primary tumors and lymph node metastases, limited prognostic information by survival analyses was found (Figure 7a-b). Still, presence of microaxons in lymph node metastases was associated with significantly shorter survival (Figure 7c). Angiogenesis markers (VPI, pMVD) in primary tumors were associated with shorter survival (Figure 7d), whereas MVD was not significant.
Vascular proliferation in lymph node metastases indicated reduced survival (p=0.05;  Table 7). A separate multivariate model for the basal-like subgroup showed a trend of prognostic value for VPI (P=0.07).
Using METABRIC data, both high SAS and high VPS (dichotomized) were associated with reduced survival in univariate analyses (Figure 7g-h). The consolidated neuroangiogenic mRNA score (NAS) was associated with shorter survival, both in univariate and multivariate analysis, adjusting for tumor diameter, histologic grade, and lymph node status (Figure 7i-j) 14 with a similar trend in non-luminal tumors (P=0.09), when adjusting for tumor diameter, histologic grade, and lymph node status.

Discussion
The influence of the microenvironment on tumor progress has been increasingly recognized (30)(31)(32). Since the present topic has not been previously studied in human cancer, we asked whether neurogenesis and angiogenesis are associated traits in breast cancer and related to more aggressive behavior. Tissue-based protein markers and gene expression data from patient cohorts indicated that both neurogenesis and angiogenesis were consistently linked to high-grade tumors and reduced patient survival. By using spatial mapping at the single-cell level, neural and vascular structures were significantly colocalized. Our novel findings indicate that neurogenesis and angiogenesis might be coordinated and closely related in aggressive breast cancer.
In this integrated study, tumor-associated neurogenesis was strongly associated with markers of high-grade tumors, such as negative hormone receptors, tumor cell plasticity, stemness, and hypoxia responses. The highest level of microaxon density was found in triple negative and HER2 positive cancers, while the lowest density was observed in the luminal A subgroup. The findings are supported by recent data (15).
Nerve bundles, assumed to be pre-existing structures, were observed in many primary tumors but did not associate with any high-grade traits. Notably, nerve bundles were not observed in lymph node metastases, where microaxons might possibly develop de novo from tissue or circulating cells (33). Patient survival was reduced in cases with presence of microaxons in lymph node metastases.
We then asked whether neurogenesis and angiogenesis are associated features in human breast cancer tissues. This was found to be the case, and a combined neuroangiogenic mRNA-based score (NAS) was strongly linked to high-grade breast cancer and independently related to reduced patient survival by multivariate analysis. Notably, the neuro-angiogenic signature was even stronger than the tumor cell-based classification. Taken together, our data support a coordinated activation of neurogenesis and angiogenesis in human breast cancer.
Spatial mapping of single cells and exploration of potential neuro-vascular niches indicated significant proximity of neural and vascular structures, and this was more pronounced in basal-like as compared to luminal-like tumors. This colocalization, by single cell-based analysis using sensitive imaging mass cytometry, support our findings using case-based protein data and transcription profiles. Our results lend support to a functional relationship between neural and vascular components, as reflected in human breast cancer samples. Experimental studies have indicated commonalities for neuronal guidance and vascular sprouting, supporting that the two processes might be co-regulated (7,8), but this has not been previously addressed in human cancer tissues.
The consolidated neuro-angiogenic score was associated with markers of epithelialmesenchymal plasticity, tumor stemness, and tumor hypoxia. These findings indicate that the connectivity of aggressive breast cancer features, based on tumor cell markers with characteristics of more invasive and less differentiated cells, are supplemented by varying degrees of "aggressive stroma". Our data suggest an existence of neurovascular interactions in the TME which could potentially be targeted for more efficient cancer management. 16 Notably, when using multiple breast cancer cell lines (34), higher levels of the neuroangiogenic score were observed in the most aggressive cancer subtypes. This indicates that neuro-angiogenic programs appear to be part of intrinsic tumor cell properties. The finding that our neuro-angiogenic score was prognostic within the luminal-like breast cancers, and the fact that some of the cases with high neuroangiogenic score were hormone receptor positive, might indicate discordant evolution of tumor cells and TME in subsets of breast cancer.
Tumor angiogenesis is an established hallmark of cancer (1) and is well described in breast tumors (35,36). Here, we found consistent and strong associations between activated angiogenesis and aggressive breast cancer. Vascular proliferation was increased in basal-like cancers, as we have reported (36,37). Our previous findings were confirmed by the present single cell analysis using imaging mass cytometry. As a novel finding, low vascular proliferation within the subgroup of aggressive basal-like cancers indicated a better prognosis compared to the others. Also, elevated vascular proliferation in lymph node metastases was associated with reduced survival, which has not been previously reported.
In summary, our data support that neurogenesis and angiogenesis are associated features of aggressive breast cancer. High neuro-angiogenic score was consistently linked to high-grade breast cancer phenotypes and reduced patient survival. Notably, among some hormone receptor positive cases, a high-grade subgroup could be defined by increased neuro-angiogenic score. Single cell tissue mapping indicated colocalization of neural and vascular structures, suggesting the presence of neurovascular niches of possible functional importance. Taken together, our findings might be relevant for improved patient stratification and exploration of novel therapy.

Patient series
The study cohort included women diagnosed with primary invasive breast cancer as 12 cases due to lack of available tissue, 534 cases were finally included as previously described (38). The patients received treatment according to the national guidelines at the time (38). Detailed information about this study population and the specimen characteristics is presented in Supplementary Material.

Immunohistochemistry (IHC)
Sections from primary tumors and the matched axillary node with the largest metastatic focus (≥ 2 mm) were used for assessment of nerves (Neurofilament antibody). For Neurofilament staining (pan-nerve marker), the monoclonal mouse antibody Neurofilament (DAKO M0762) was used.
Sections from the same block of the primary tumor and the axillary lymph node with metastases that were used for assessment of tumor-related nerves, were applied for evaluation of angiogenesis markers; the polyclonal rabbit antibody (DAKO A0082) for Factor-VIII, and a monoclonal mouse antibody for Ki67 (DAKO 7240), were used. The blood vessels in the adjacent normal mammary or fat tissue stained with Factor-VIII and the proliferating tumor cells stained with Ki67 were used as internal controls for 18 these markers. In lymph nodes, the adjacent capillaries and the proliferating lymphoid cells were used as internal control cells for Factor-VIII and Ki67, respectively. The comparison with HE slides was always done to guide the location of primary tumors and metastatic foci in the breast and affected lymph nodes. For further details from the immunohistochemistry, see Supplementary Material.

Evaluation of tissue-based IHC for neurogenesis and angiogenesis markers
Microaxon and nerve bundle density. Microaxon and nerve bundle densities for each case were assessed separately in the whole tumor area using light microscopy (MAD=microaxon density; NBD=nerve bundle density). Microaxons were observed as single, small, and thin structures (Figure 1a), while a nerve bundle was identified as a group of axons (Figure 1b). The number of microaxons or bundles per mm 2 were determined.
Microvessel density (MVD). As a measure of angiogenesis, the count of all vessels in ten high power visual fields (HPF x400) was performed as described (35), and reported as counts/mm2. Microvessels included both vessel-like structures with a visible lumen and single endothelial cells or cell clusters (Figure 1c), as defined by Factor-VIII positivity according to Weidner (39).
Proliferating microvessel density (pMVD). Using the same ten microscopic fields (x400) that were used for microvessel density counts, vessels containing Ki67 positive endothelial cells were counted. Thus, pMVD represents a separate count of vessels with proliferating endothelial cells and was reported as count/mm² (Figure 1d).

Vascular proliferation index (VPI).
For estimation of activated angiogenesis, the proportion (%) of vessels with proliferating endothelial cells of the total number of microvessels counted in ten fields (Vascular proliferation index; VPI) was calculated.

IMC staining protocol
Antibody hybridization was performed as described in the "Imaging Mass Cytometry Staining Protocol for FFPE Sectioned Tissue" (Fluidigm) with a few modifications. In short, the freshly cut TMA slide underwent dewaxing, rehydration and antigen retrieval in a Ventana Ultra Discovery Autostainer (Roche Diagnostics GmbH). The antibody mix was applied to the slide, which was then stored overnight at 4°C in a hydration chamber. After antibody incubation, the slide was first washed with PBS with and without detergent, before incubation first in 0.3 uM Iridium (Ir)-intercalator (Fluidigm), washed, and then in 0.0005% Ruthenium (RuO4)/PBS (Electron Microscopy Sciences) (40). Finally, the slides were washed in MaxPar H2O, air-dried, and stored at 4 °C.

Antibody validation
All antibodies were validated by IHC using a test-TMA with positive control tissues for all antibodies included in the panel: tonsil, placenta, hippocampus, cerebellum, autonomic ganglion and peripheral nerve tissue, normal breast tissue, and selected breast carcinomas (ER+/PR+, HER2+, ER-/PR-/HER2-and basal-like breast carcinomas). In some cases, our IHC staining was compared to staining performed by 20 the Human Protein Atlas (proteinatlas.org). A pilot-TMA was also established and stained (five basal-like; five luminal-like).

IMC analysis and pre-processing
Data from the 10 pilot-TMA breast cancer cases were acquired by a Helios time-offlight mass cytometer (CyTOF) coupled to a Hyperion Imaging System (Fluidigm) and administered using the CyTOF Software (v7.0.8493; Fluidigm). The square inscribed in the circular TMA cores were laser ablated at 200 Hz at a resolution of approximately 1µm 2 .
The highly multiplexed data was processed using the ImcSegmentationPipeline(41) to segment data into spatial regions corresponding to single cells. First, the IMC output (.txt files) were converted to tiff-images (ImcTools v2.1.7), which were used for pixel classification (Ilastik v1.3.3post2) to generate probability maps of nuclei, cytoplasm/membrane, and background. Each probability map, corresponding to one TMA core, were then segmented into single cells, with output as a cell mask image, using Cell Profiler (v4.0.7). The pre-processed data were imported into Histocat (v1.7.6; (42)) and immediately exported as csv-files to be used in down-stream analyses. We note that even with high-quality segmentation, the imaging of tissue segments cannot rule out cases of "overlapping cell units" that do not capture the nucleus of an individual cell. Therefore, in a few cases nuclei-mismatched signal can be assigned to neighboring cells, especially in (cellular) dense areas.

IMC data normalization and analysis workflow
For each marker, the total intensity per cell unit was computed. Values were divided by the cell size and were arcsinh-transformed. The single-cell data per channel were censored at the 99 th percentile to remove outliers. Single-cell annotation was 21 performed in a hierarchical scheme using unsupervised clustering/meta-clustering and prior-knowledge of cell type defining markers of our antibody panel, similar to the methodology presented in (17). The procedure operates in two steps: initially, single cells are categorized as Immune and Non-Immune. FlowSOM (43) is used to cluster the data into 120 clusters, that further merged based on cluster profile similarity. After this step, Phenograph (44) is used to refine the annotation and split Immune cells to three groups (B cells, T cells, macrophages), and non-Immune cells to three groups (epithelial cells, endothelial cells, stromal cells). For the clustering approach, data were standardized across markers using standard normalization.
To characterize vascular proliferation, we focused on the endothelial phenotype. All endothelial cells were gated into four subgroups (Ki67-low/CD31-high, Ki67high/CD31-high, Ki67-high/CD31-low and Ki67-low/CD31-low) on the basis of CD31 and Ki-67 marker expression, using the flowDensity software (18). flowDensity is a supervised clustering algorithm based on density estimations that is designed to overcome the tedious work of setting static gates. flowDensity was also used to stratify single cells to Neurofilament high and Neurofilament low expression groups, given the density estimation of the Neurofilament marker, see Supplementary Material for further details on the evaluation of the tissue-based neurogenesis and angiogenesis markers.
To estimate the distance of endothelial cells, to Neurofilament-high cells, we used a knearest neighbor approach. For every cell of the Ki67-low/CD31-high, Ki67-high/CD31high, Ki67-high/CD31-low and Ki67-low/CD31-low phenotype, we first selected the kclosest cells. From the set of k-closest cells we filtered the ones with Neurofilament high expression, and we estimated the Euclidean distance from the corresponding endothelial cells of interest. The average log2-transformed distance was used for visualizations, and the procedure was repeated for different k values starting from 3 to 10 to assess reproducibility.

Proteomic analyses of microdissected breast cancer epithelial cells
Proteins were extracted using the filter-aided sample preparation (FFPE-FASP) protocol (47) and analyzed by liquid-chromatography tandem mass spectrometry (LC-MS/MS). Raw data were analyzed using the freely available Maxquant software (v1.6.0.16) with recommended settings for label-free quantification (48) (see Supplemental Material for an extended description of methods). Gene ontology datasets were collected from the AmiGO 2 resource (version 2.5.12). Protein intensities were compared between groups (basal-like and luminal-like) using the Students t-test. Correction for multiple hypothesis testing was done with the Benjamini-Hochberg method (49). Differential analysis of protein expression was performed using the DEP R package (50). Volcano plots and other visualizations of the proteomics results were performed in R.

Gene expression (mRNA) analyses of breast cancer tissues and cell lines
In  Other gene expression signatures. For correlation analyses, we mapped two hypoxia signatures (28,29), an EMT signature (21), a TGF-β signature (22), and two signatures reflecting stemness features (25,26), to the mRNA data we analyzed. For further detailed information about gene expression analyses, see Supplementary Material.

Statistics and survival analyses
Data were analyzed using SPSS (Statistical Package of Social Scienc1es), Version 25.0 (Armonk, NY, USA; IBMM, Corp) and R-package via Rstudio interface version1.3.9. A two-sided P-value less than 0.05 was considered statistically significant. A P-value of 0.05-0.10 was considered to represent a statistical trend.
Categories were compared using Pearson`s chi-square or Fisher`s exact tests when appropriate. Non-parametric correlations were tested by Spearman's rank correlation, while Mann-Whitney U and Kruskal-Wallis tests were used for comparing continuous variables in groups. Wilcoxon Signed Rank test was used to compare the differences between two continuous variables. The Kolmogorov-Smirnov statistic was used to quantify a distance between the between the empirical distribution functions of two samples. Odds ratios (OR) and their 95% confidence intervals were calculated by the Mantel-Haenszel method. Kappa statistics were used to test inter-and intra-observer agreement of categorical data. For survival analyses, the endpoint was breast cancer specific survival, defined as the time in months from the date of histologic diagnosis to the date of death of breast cancer. Univariate survival analysis (Kaplan-Meier method) was performed using the log-rank test to compare differences in survival time between categories. Patients who died of other causes or who were alive at last date of followup were censored in the analyses. The influence of co-variates on breast cancer specific survival was analyzed by Cox`s proportional hazards method and tested by the backward stepwise likelihood ratio (lratio) test. All variables were tested by logminus-log plots to determine their ability to be incorporated in multivariate modelling.
When categorizing continuous variables, cut-off points were based on median or quartile values, also considering the distribution profile, the size of subgroups, and number of events in survival analyses.

Study approval
This study was approved by the Western Regional Committee for Medical and Health Research Ethics, REC West (REK 2014/1984). Written informed consent was not obtained from the patients, but in accordance with national ethics guidelines and procedures for retrospective studies, all participants were contacted with written information on the study and asked to respond if they objected. In total, 9 patients (1.7%) did not approve participation.

Author contribution
LAA conceived and designed the study and contributed to histologic analyses (IHC, IMC), data interpretation and in writing the major parts of the manuscript. EW contributed to the study design, analyzing the transcriptomic data, general data 26 analyses and interpretation, and in writing the major parts of the manuscript. SA contributed to the study design, performed the tissue-based work and IHC, performed statistical analysis, data interpretation and contributed to writing the manuscript. DK contributed to study design and analyzing the mRNA and IMC data, and contributed to writing the manuscript. HV and KF contributed to analyses of proteomics and IMC data.
CA contributed to histologic analyses and work on IMC data, and to collection of clinicopathologic data. IW contributed to IMC data. GK participated in data collection and interpretation. KK, AAS, EB, SK, and IMS participated in data collection and interpretation. TA and BD participated in collection of clinical data. All authors commented on and approved the final manuscript. EW and SA made the initial core findings, whereas EW played a larger role in manuscript completion; thus, EW is listed first in the authorship. All authors have read and approved the submitted paper.     Clustering was conducted on binary variables reflecting either high (red color) / low (gray color) signature scores, positivity (red color) / negativity (gray color) for ER and HER2, or presence (red color) / absence (gray color) of histologic grade 3, and basal-like phenotype from PAM50. Four cluster groups were evident visually, and the tree was cut horizontally to further explore these clusters marked as C1, C2, C3, C4; c) Kaplan-Meier curve illustrating the difference in survival between patients in the cluster groups C1-C4 using data from the METABRIC discovery cohort from (a). Differences in breast cancer specific survival between groups were tested using the log-rank test.
Neuro−angiogenic score EMT score (Jechlinger 2003 C-4 C-3 C-2 C-1 C-4 C-3 C-1 C-2 Figure 7. Markers of angiogenesis, neurogenesis, and neuro-angiogenesis and their association with survival in breast cancer. a-b) Microaxon density (MAD), and nerve bundle density (NBD) do not associate with survival in primary breast cancers; c) Presence of microaxons in lymph node (LN) metastases associates with reduced breast cancer specific survival; d-f) High vascular proliferation index (VPI) in primary breast cancer, and LN metastases associates with reduced cancer specific survival, also in the subset of primary basal-like breast cancers; g-h) High Sprouting axons score (SAS), Vascular proliferation score (VPS), and Neuro-angiogenic score (NAS) associate with reduced breast cancer specific survival in patients from the METABRIC cohort. In all Kaplan-Meier analyses the log-rank test was used to quantify differences. Next to every survival curve we report the number of breast cancer deaths, out of the total number of cases in the category; j) Cox' proportional hazard modelling shows that the Neuro-angiogenic score demonstrates independent prognostic value in multivariate survival analysis, adjusting for tumor size, histologic grade, and LN status, using data from the combined METABRIC cohort.