An Endothelial Gene Signature Score Predicts Poor Outcome in Patients with Endocrine-Treated, Low Genomic Grade Breast Tumors

Breast cancer remains the most common malignancy and a leading cause of cancer-related death in women (1). Decreasing mortality rates in recent decades have come at a cost of both more extensively applied toxic adjuvant therapies, stressing the importance of finding reliable prognostic and treatment predictive markers. The success of administering trastuzumab to patients with overexpression of the HER2 protein serves as an excellent example (2,3), but equally useful predictive markers for treatment response are largely lacking for other targeted therapies, including antiangiogenic therapies like bevacizumab and sunitinib.

The last 15 years have seen the application of microarray technology to tumor samples with the aim of finding better prognostic and treatment predictive strategies for breast cancer patients. This research has resulted in a plethora of genomic classifiers ranging from binary good/poor prognosis signatures (4–6) to multilevel classifiers capable of diving breast cancer into prognostically relevant molecular subgroups (7). Further studies have tested the predictive capacity of gene signatures for both chemotherapy and tamoxifen (8–10) with favorable findings, although given the central role of proliferation-related genes in many classifiers (11,12) their value over traditional immunohistochemical markers such as Ki67 remains unclear.

Recent clinical trials into antiangiogenic therapies (13–15) have served to once again highlight the belief that breast cancer progression can be impeded through targeting of tumor angiogenesis. Despite this continued interest in these therapies, and question marks over their ability to prolong overall patient survival, the clinical relevance of angiogenesis-related transcription as a treatment predictive factor in breast cancer has remained largely unexplored. Here, we aim to examine a previously published set of 57 gene transcripts (ref. 16; representative of a normal endothelium) first through comparison with the traditional microvessel density (MVD) by correlation and gene ontology (GO) analysis and finally through assessment of the signature in six different breast cancer cohorts, with particular focus on its performance in subgroups of endocrine and chemotherapy-treated patients.

We have previously described both the Uppsala (N = 253) and Stockholm (N = 159) datasets (17–19) with an extensive overview for both cohorts found here (12). Both microarray studies were approved by the ethics committees at Karolinska Institutet and Karolinska University Hospital (Stockholm, Sweden), respectively and are publically available at NCBI GEO under accession numbers GSE4922 and GSE1456, respectively.

Data from the Netherlands Cancer Institute (NKI; N = 295; ref. 20), Erasmus Medical Center (Rotterdam; N = 286; ref. 21), and the John Radcliffe and Guys hospitals (Oxford and London, United Kingdom; N = 99 and 87, respectively; ref. 22) were used, the Oxford and Guys data for further analysis of the findings in relation to endocrine therapy. The NKI dataset is publicly available as a part of the breastCancerNKI R package and the Rotterdam and Oxford/Guys datasets are available under accession numbers GSE6532 and GSE5432. Clinical information for 14 patients treated with sunitinib plus docetaxel or docetaxel alone has been previously described (ref. 23; ClinicalTrials.gov identifier NCT00393939) and expression data are publicly available at NCBI GEO under accession number GSE54323. For the comparison of MV score to lymphovascular invasion (LVI), we used previously published data (24), retrievable under the accession number GSE5420.

RNA was extracted from homogenized tumor material with RNeasy spin column kits (Qiagen) and quality was assessed with an Agilent 2100 Bioanalyzer. Two to 5 μg of RNA was used to produce biotinylated cRNA. Hybridization to HGU133 A and B microarrays (Affymetrix) and scanning was performed according to Affymetrix protocols. A comprehensive account can be found in Pawitan and colleagues (19).

For immunohistochemical analysis, formalin-fixed paraffin-embedded sections (4 μm) were deparaffinized in xylene and rehydrated in graded concentrations of ethanol to TBS. Antigen retrieval by microwave treatment was performed for 20 minutes in Tris-EDTA Buffer (10 mmol/L Tris Base, 1 mmol/L EDTA Solution, 0.05% Tween 20, pH 9.17). CD31 antibody (clone JC70A; DAKO) was diluted 1:50, and staining was carried out in a Tech Mate Autostainer (DAKO). Slides were counterstained with hematoxylin and dehydrated. MVD was determined using previously described methods (25). Briefly, tumor sections were examined at low power to determine the areas containing the greatest numbers of microvessels. Individual microvessels were then counted in these areas at ×400 magnification (three fields per tumor section). Mean counts per high field were then calculated and the resulting value was normalized to yield an MVD score between 0 and 100 (for ease of comparison to microvascular gene expression score) using the rescale command of the “scales” package in R (26).

Human homologs for the 57 of 58 previously published mouse microvasculature transcripts (16) were extracted from the HomoloGene database at the National Center for Biotechnology Information (NCBI, Supplementary Table S1). Datasets were RMA normalized and median centered before the expression levels of signature genes were added per tumor, and the resulting sum was scaled to yield a microvasculature signature score (hereafter called MV score) of between 0 and 100 within each dataset using the rescale function of the R genefu package. The exception to this was the NKI dataset; here, we used the original normalized and median-centered data from the breastCancerNKI R package for scaling as described above. Of note, mean centering has been demonstrated as sufficient to remove much of the dataset bias associated with gene expression data from different cohorts, allowing for meaningful prognostic comparisons to be made (27). Fifty-seven of the genes (corresponding to 115 probes) were present on the HGU133 A and B platforms, and could be used for MV score determination in the Uppsala and Stockholm data. In the external datasets, identical methodology was used to determine MV scores and for the NKI, Rotterdam, Oxford, Guys, and metastatic tumor data, 49, 46, 57, 57, and 57 signature genes were present on the respective platforms corresponding to 56, 79, 115, 115, and 115 probes, respectively. In the case of multiple probes mapping to the same gene, an average expression of probes was taken. When data were pooled for combined analysis, the signature score was calculated and scaled in individual datasets before pooling.

Molecular subtyping according to the PAM50 signature was performed as outlined in the original publication (7), using the code provided by Parker and colleagues on the UNC Microarray Database website as a data supplement to the original article. Of note, we have previously published our code and the PAM50 subtypes for the Uppsala and Stockholm cohorts (12).

Genomic grade index (GGI) was calculated as described in the original publication (6) and we have previously published our code for this signature along with the GGI calls for the Uppsala and Stockholm datasets (12).

To assess the similarity between the MV and MVD scores, we performed a Pearson correlation comparing (i) both scores directly and (ii) each of the MV score genes to MVD score in 182 patients of the Uppsala cohort. P values were adjusted for multiple testing using the FDR method as part of the p.adjust command of the R stats package. To test for enrichment of biologic themes, the significance analysis of microarrays (SAM; ref. 28; quantitative response) was used to rank genes for association with the MV score (gene expression) and MVD scores (vessel count), respectively. Over-representation of gene ontology categories (GO) reflecting biologic processes was determined with conditional hypergeometric tests using the R GOstats package. As many of the MV score genes are annotated as vascular or angiogenesis related, and by necessity will correlate to the signature, the 57 signature genes were removed from the data prior to GO testing. Estimation of false discovery rate control was performed according to Storey and Tibshirani (29) using the R package qvalue 2.0.0. All occurring categories in respective gene lists were tested.

To test for differences in mean expression of the MV score in relation to clinicopathologic parameters, Student t test (for comparing means between two groups) or ANOVA (for comparisons between more than two groups) with post hoc Tukey were used as indicated in table legends. Survival analysis was performed using Kaplan–Meier and multivariable proportional hazards (Cox), where only variables with demonstrated prognostic significance in univariable testing were included in the final multivariable model. To ensure consistent survival endpoints across all datasets, distant metastasis-free survival (DMFS) was chosen. Hazard ratio (HR) is reported per 25 increments of MV score and MVD score, to simplify interpretation. ROCs were used to determine best MV score cutoff (with DMFS at 5 years as endpoint) in testing datasets before being applied to validation datasets and correlograms were constructed using the Spearman rank correlation metric. Concordance index (c-index) calculation was performed in R using the “concordance” output from a Cox regression analysis for the GGI and PAM50 signatures alone or in combination with the MV score. All analyses were performed using the R Statistical language (26), with the GOstats, ROCR, survival, and samr packages.

To understand if any similarities exist between genes associated with microvasculature gene expression and traditional MVD score, we stained and scored the tumors of 182 patients from our previously published Uppsala cohort (12,17–19) using the endothelial/MVD marker CD31. Next, we derived a normal MV score using gene expression data from the same tumors by adding the mRNA expression levels for 57 of 58 previously identified (16) endothelial-specific gene transcripts (Supplementary Table S1). Both the MV and MVD scores were then normalized and scaled (for ease of comparison) to yield values between 0 and 100.

A simple correlation analysis did not reveal any clear overall similarity between the MV and MVD scores (Pearson correlation = 0.096, P = 0.20, data not shown) and in individual signature gene analysis, only a weak correlation was found between five MV score genes and the MVD score (Supplementary Table S1). Next, we compiled two lists of the top 200 genes most associated with the MV score (having first removed the 57 signature genes) and MVD scores and assessed whether any biologic processes were over-represented in these lists through GO analysis. The top 10 GO terms for both lists are displayed in Table 1, where an enrichment of terms related to angiogenesis and cardiovascular/blood vessel/endothelial development is notable among MV score–correlated genes (Table 1, top list). Conversely, no angiogenic or vascular development terms are found within MVD (CD31) correlated genes, which rather display an enrichment of terms associated with immune response (Table 1, bottom list). First, these results show that the MV score strongly reflects angiogenic/endothelial processes on a transcriptional level and second, based on the CD31 GO terms, we may be unaware of the extent to which the immune response is involved in blood vessel formation and maintenance.

With the aim of determining the prognostic and treatment predictive capacity of the MV score, we extended our analysis to three additional gene expression breast cancer datasets (12,19–21). Again, MV scores were normalized and scaled within each dataset and the resulting score was assessed in relation to the clinicopathologic parameters shown in Table 2. 

We noted a trend toward lower MV scores in grade 3 tumors across all datasets (vs. grade 1 tumors; Table 2) and similarly, a statistically significant low signature score in the estrogen receptor–positive Luminal B subtype (vs. Luminal A; Table 2, P ≤ 0.001, 0.008, 0.016, and <0.001, for the Uppsala, Stockholm, NKI, and Rotterdam datasets, respectively). A trend toward a higher MV score was found in tumors of a Normal-like subtype, which reached statistical significance in two of four datasets (vs. Luminal A subtype; Table 2, P = 0.003 and 0.047, for the Stockholm and Rotterdam datasets, respectively). Of interest, as information on lymphovascular invasion (LVI) was not available in these cohorts, we analyzed an additional publicly available dataset (ref. 24; n = 128) and could not demonstrate a difference in MV score in the absence or presence of LVI (MV score mean ± SEM = 45.71 ± 2.23 and 49.54 ± 2.15, for tumors with and without LVI, respectively. P = 0.22, Welch t-test, data not shown).

The prognostic capacity of many first-generation gene expression signatures tends to be limited to ER-positive breast tumors (30). As such, we assessed the MV score in univariable analysis across all, ER-positive, and ER-negative patients in each dataset. No consistent relationship to distant DMFS was found in any of these groupings; however, in all patients, two of four datasets demonstrated a lower HR with increasing MV score (Supplementary Table S2, see “All patients” Uppsala and Stockholm cohorts, HR 0.69; 95% CI, 0.49–0.96 and 0.46; 95% CI, 0.24–0.87, respectively).

To determine whether the type of adjuvant therapy received influenced these results, we subdivided each cohort into patients who did not receive systemic treatment, who received endocrine treatment, or those who received chemotherapy. The MV score did not provide consistent statistically significant information regarding DMFS in untreated, untreated ER-positive, or chemotherapy-treated patients (Fig. 1, see “Untreated,” “Untreated (ER positive),” and “Chemotherapy,” All). For patients receiving endocrine therapy, a trend toward decreased risk of distant metastasis was observed in two of three datasets (Fig. 1, see “Endocrine treatment,” Uppsala and Stockholm, HR 0.69; 95% CI, 0.42–1.12 and 0.45; 95% CI, 0.20–1.03, respectively); however, this trend did not reach overall statistical significance (Fig. 1, see “Endocrine treatment,” All, HR 0.82; 95% CI, 0.57–1.18). As the HR for the NKI dataset in endocrine-treated patients (Fig. 1, HR 1.56; 95% CI, 0.36–6.6) was in the opposite direction to that of the Uppsala and Stockholm datasets, and as 50% of the endocrine patients in the NKI datasets also received chemotherapy, we further examined the MV score in an independent dataset of 186 patients collected at the John Radcliffe and Guy's Hospitals (Oxford/Guys dataset; ref. 22). All patients had ER-positive tumors and received adjuvant tamoxifen monotherapy. Here, a higher MV score was associated with a reduced risk of DMFS in endocrine-treated patients (Fig. 1, see “Additional datasets,” Both, HR 0.58; 95% CI, 0.38–0.89), consistent with the trend found in Uppsala and Stockholm endocrine-treated patients. This statistical significance remained in multivariable analysis combining the Uppsala, Stockholm, Oxford, and Guys data, when considering standard prognostic markers (Table 3, n = 380, see “Histologic Grade,” MV score HR 0.70; 95% CI, 0.50–0.97). Taken together, these data suggest an endocrine therapy predictive capacity for the MV score.

Although numerous gene expression signatures exist for endocrine-treated patients (5,9,31), their prognostic capacity relies on proliferation-related genes (12,32). To address a potential relationship to proliferation for the MV score, we classified the Uppsala, Stockholm, and Oxford/Guys datasets according to the GGI, a strongly proliferation-related gene signature (12,22). When histopathologic grade was replaced by genomic grade in multivariable analysis, the prognostic ability of the MV score was lost (Table 3, see “Genomic Grade,” MV score HR 0.89; 95% CI, 0.65–1.21). Similarly, neither MV score nor the PAM50 subtypes were prognostic in a multivariable analysis containing both variables (Supplementary Table S3, MV score HR 0.84; 95% CI, 0.58–1.21). To further examine the reason for this loss of prognostic power, we first identified a cutoff (of 39) for the MV score in our training datasets (Uppsala and Stockholm cohorts) that would perform best for prediction of DMFS in endocrine-treated patients (Supplementary Fig. S1A, iv). Second, we used this cutoff to produce Kaplan–Meier curves for both our training and validation datasets (Oxford and Guys cohorts), split according to genomic grade (GG1 or GG3). This analysis showed that the strength of the MV score resides in the lowly proliferative genomic grade 1 group of tumors (Fig. 2A, GG1, training group, MV score high/low n = 66/48; Fig. 2B, GG1 validation group, MV score high/low n = 84/26; and Fig. 2C, GG1 both groups, MV score high/low n = 150/74; P = 0.31, 0.012, and 0.007, respectively. Compare to Fig. 2D–F, GG3, training, validation, and both groups together, respectively.). In this subgroup, as in the full set of endocrine-treated tumors, the MV score had independent prognostic capacity over standard prognosticators (Supplementary Table S4, HR 0.47; 95% CI, 0.26–0.84). It is pertinent to highlight here that although the separation of GG1 curves does not reach formal statistical significance in the training dataset (Fig. 2A), there is a crossing of curves after approximately 4 years, potentially rendering the log-rank test underpowered. Moreover, this cutoff displays strong prognostic capacity in other treatment subgroups of the training dataset (Supplementary Fig. S1B, i–iii). Next, we hypothesized that a combination of both signatures (MV and GGI) would provide more prognostic information than either one alone. To test this, we calculated the c-index for the GGI, MV, and PAM50 signatures alone and in combination in all patients of our original four datasets. In line with our hypothesis, we found that, in general, the addition of the MV score to the PAM50 and GGI gene signatures may provide more prognostic information all datasets than either signature alone (Supplementary Fig. S2, compare green vs. blue bars—PAM50 vs. PAM50 + MV score and red vs. yellow bars—GGI vs. GGI + MV score, in all datasets).

To further explore the relevancy of our signature in a clinical setting, we calculated the change in MV score in a previously published cohort of 14 metastatic breast tumors (7 from the control arm and 7 from the treated arm) before and after treatment (14 days) with the angiogenic inhibitor sunitinib (23). These samples were taken as part of a substudy from a recent phase III clinical trial comparing the efficacy of sunitinib and docetaxel versus docetaxel alone in an advanced breast cancer setting (33). A heatmap displaying the intrapatient change in the expression of the MV score genes before and after treatment in both clinical trial arms is shown in Fig. 3A, where in general, a greater change in signature genes was found in the combination arm (Fig. 3A, red bar). Concomitantly, a significant increase in MV score was noted after treatment in the combination arm (Fig. 3B, right, P = 0.031 vs. baseline); however, low patient numbers prevent further analysis regarding survival. For the sake of completeness, we also show a table of the change in MV score (14 days, baseline) and RECIST response for all patients, no clear trend is observable (Supplementary Table S5). These results are in line with the concept of vascular normalization following treatment with an angiogenic inhibitor (34,35), but notably, are the first demonstration of this principle on a transcriptional level in human breast tumor samples.

In this study, we found that the abundance of normal microvascular transcripts was reproducibly related to both the Luminal B breast cancer subtype and the clinical endpoint DMFS. In 993 primary breast carcinomas, a simple summary signature was expressed at lower levels in Luminal B tumors and in endocrine-treated patients, high expression of the MV score displayed a trend toward a more favorable outcome in two of three datasets and a similar finding was noted in a set of 186 patients subjected to tamoxifen monotherapy. Multivariate and subgroup analysis revealed that this association was only present in a subgroup of tumors characterized by low genomic grade. In addition, a significant increase in signature score was found in seven metastatic breast tumors after 14 days treatment with sunitinib + docetaxel, an increase that was not present in tumors treated with docetaxel alone.

Although others have also produced microvascular gene expression signatures in a breast cancer setting, these signatures have generally been designed to capture the transcriptional differences between microdissected normal versus tumor microvasculature. In short, these signatures likely represent a tumor endothelium that is increasingly thought of as angiogenically active and chronically inflamed (36). This is in contrast to our signature score that is derived from a physiologically normal microvasculature and as such is highly expressed in low-risk tumors. These differences are further emphasized when comparing the overlap of our module genes with other published angiogenesis signatures: number of genes from our module present in Bhati and colleagues signature, 1/48; in Masiero and colleagues signature 13/43; Pepin and colleagues 0/494; and Mannelqvist and colleagues 0/18 (data not shown; refs. 37–40). Of note, two recently published endothelial metagenes do display a greater degree of overlap with our signature, those being “signature 4” and “signature 5” from Winslow and colleagues (ref. 41; overlap = 4/9 and 2/3, respectively). Interestingly, these metagenes were derived through correlation analysis to a core set of genes enriched in tumor stromal compartments. Related to this, we also characterized the MV score in the context of other published signatures and gene expression modules (42) and found that that the MV score (Supplementary Fig. S3A–S3D, MV_SIG, red arrows) is strongly inversely correlated to the AURKA proliferation–related gene module when considering all patients (Supplementary Fig. S3A–S3D, AURKA, blue arrows). Of note, although we also see moderate correlations to the GGI, PLAU, and Stoma1 modules in all patients, these become weakly correlated or not statistically significant in the endocrine-treated, low GG subgroup (data not shown). A second link between our signature and proliferation was in evidence when we examined the MV score within the PAM50 molecular subgroups. Here, the MV score was consistently lower in the Luminal B tumors of all four tested datasets relative to Luminal A tumors. It has previously been demonstrated that the one of the main factors distinguishing these two tumors groups is level of proliferation with higher levels found in Luminal B tumors (43). However, given that our signature retains prognostic significance in the lowly proliferative GG1 tumor subgroup (Supplementary Table S4), it is reasonable to state that the prognostic capacity of signature extends beyond that of a simple proliferative marker.

Given the vast morphologic differences between the endothelium lining, the cardiovascular/lymphatic systems, and tumor endothelium, disparities in the quantity and type of genes expressed are to be expected. The tumor endothelium is characterized by atypical cell morphology, blood flow that can range from chaotic to nonexistent (44), and intracellular gaps that leak fluids and blood into the surrounding tissue (45). Taken together, these hallmarks of endothelial dysfunction not only influence gene expression patterns but also make pathologic assessments of microvascular density challenging. Indeed, spotty CD31 staining has been highlighted in the tumor endothelium in vivo, owing to a lack of expression in some cells and an absence of cells entirely in some areas of the vessel wall (46). This is likely one of the reasons as to why we found no overt similarity when comparing the MV score to microscopic assessment of MVD with CD31 staining. Staining issues notwithstanding, the value of MVD assessment as a prognostic marker has, on the whole, been called into question. In a systematic review of MVD and outcome, Uzzan and colleagues reported risk ratios in the range of 1.5 to 2, concluding that MVD has significant but weak prognostic capacity in breast cancer, and that standardization of MVD assessment is needed (47). Similar weak/negative findings were recently reported by Cheng and colleagues in a recent renal cell carcinoma meta-analysis of MVD (48).

This study had some limitations, the foremost of those being that this is a retrospective study performed in multiple patient cohorts (rather than a single, large, homogeneously-treated cohort) and that the patient numbers in the metastatic cohort are low (N = 14), as expected from a feasibility study. Furthermore, not all analyses were prespecified; our initial aim was to characterize the MV score in terms of its relationship to traditional MVD and to standard breast cancer clinicopathologic parameters in the context of treatment subgroups and DMFS. The subgroup analysis splitting patients into GG1 and GG3 and the examination of change in signature score before and after treatment in the metastatic cohort were exploratory in nature.

In summary, we report a MV score representative of a normal endothelium that reproducibly describes differential expression between Luminal A and B molecular subtypes, and identifies a subgroup of endocrine-treated patients with worse outcome. Moreover, we show the first evidence of normalization of tumor vasculature on a transcriptional level in response to an angiogenic inhibitor in human metastatic breast cancer samples. In light of these findings, evaluation of transcriptional changes in microvascular genes alongside assessment of MVD and angiogenic factors in clinical trials of antiangiogenic compounds appears warranted.

J. Bergh reports receiving research funding, through Karolinska University Hospital or Karolinska Institutet, from Astrazeneca, Amgen, Bayer, Merck, Pfizer, Roche, and Sanofi-Aventis. T. Foukakis reports receiving a commercial research grant, Karolinska Hospital, from Roche. No potential conflicts of interest were disclosed by the other authors.

Conception and design: N.P. Tobin, K. Wennmalm, L. Lindstrom, C. Betsholtz, J. Bergh

Development of methodology: N.P. Tobin, K. Wennmalm, T. Foukakis

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K. Wennmalm, T. Foukakis, G. Landberg, J. Bergh

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N.P. Tobin, K. Wennmalm, L. Lindstrom, T. Foukakis, L. He, G. Genove, A. Ostman, G. Landberg, C. Betsholtz

Writing, review, and/or revision of the manuscript: N.P. Tobin, K. Wennmalm, L. Lindstrom, T. Foukakis, G. Genove, A. Ostman, C. Betsholtz

Study supervision: K. Wennmalm, C. Betsholtz

The authors thank Dr. John Lövrot, Dr. Johanna Klinge, and Julie Lorent, for essential technical assistance and discussion.

This work was supported by grants from the Swedish Cancer Society (to J. Bergh, C. Betsholtz, G. Genové), the Stockholm Cancer Society (to J. Bergh), BRECT (to J. Bergh), the King Gustav V jubilee fund (to J. Bergh), the ALF Foundation (to J. Bergh), the Linné Foundation (to J. Bergh), the Swedish Research Council [to J. Bergh, C. Betsholtz, LSL (grant nos.: 524-2014-4483 and 521-2014-2057)], the Swedish Research Council for Health, Working life and Welfare, FORTE [LSL (grant no.: 2014-1962)], the Märit & Hans Rausing Foundation (to J. Bergh), the Knut and Alice Wallenberg's Foundation (to J. Bergh, C. Betsholtz), and AstraZeneca (to J. Bergh).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.