Identification of Interacting Stromal Axes in Triple-Negative Breast Cancer

In clinical practice, breast cancer is stratified into subtypes defined by differential protein expression of the estrogen receptor (ER) and progesterone receptor (PR), as well as expression and/or genomic amplification of human epidermal growth factor receptor 2 (HER2). Together with other clinical variables including patient age, tumor size, grade, and lymph node status, these factors form the basis for clinical decision-making with respect to treatment (1). The 10%–15% of breast carcinomas that lack expression/amplification of ER, PR, or HER2 form a subtype that is termed triple-negative breast cancer (TNBC; ref. 1). Tumors of the TNBC subtype are associated with earlier age of onset, higher grade at presentation, and overall poorer patient prognosis (2). Despite extensive efforts, this subtype still lacks targeted therapies and these tumors are generally treated with untargeted chemotherapy and radiation (3, 4).

Previous studies, including several high-throughput profiling efforts, have indicated that the TNBC subtype has higher levels of intertumoral (patient-to-patient) heterogeneity when compared with other subtypes with respect to both gene expression (5), and somatic genomic aberrations (6, 7). This heterogeneity may at least partially underlie why TNBC is a poor outcome subtype (8–11). Several efforts have investigated whether there are subtypes within TNBC with distinct cellular processes and responses (12–14). These studies, however, have used gene expression profiling of bulk samples enriched for epithelial cells of the tumor proper. Lehmann and colleagues (13) identified six subtypes within TNBC patients: BL1 (basal-like 1), BL2 (basal-like 2), M (mesenchymal), MSL (mesenchymal stem-like), IM (immunomodulatory), and LAR (luminal androgen receptor) – and demonstrated that these “TNBCType” subtypes are associated with differential responses to neoadjuvant chemotherapy (15). This scheme has since been translated into an assay for clinical implementation (INSIGHT TNBCTYPE; ref. 16).

It is well-established that the microenvironment of the tumor coevolves with the tumor proper, and that this coevolution is central in tumor progression, and ultimately patient outcome (17). To date the few studies that have investigated stromal gene expression profiles have utilized a pan-breast cancer cohort (18–20). Our effort here is the first to investigate stromal heterogeneity specifically within TNBC using microdissected stromal compartments of both the tumor proper and adjacent morphologically normal cells from 57 TNBC clinical samples. An unbiased class discovery approach identified four stromal axes with distinct molecular and clinical qualities associated with the TNBC microenvironment. Combining the stromal axis scores to generate a novel subtyping approach revealed interactions between the stromal axes with distinct prognostic associations.

Samples were collected from patients undergoing breast surgeries at the McGill University Health Centre (MUHC) between 1999 and 2012 who provided written, informed consent (MUHC REB protocols SDR-99-780 and SDR-00-966). All tissues were snap-frozen in O.C.T. Tissue-Tek Compound within 30 minutes of removal. Information regarding clinical variables was obtained through review of medical records. Samples used for this cohort (n = 57) were reported as negative for ER (based on a 10% threshold) and HER2 via IHC (ER and HER2) and/or FISH (HER2). All patients were PR-negative (IHC), with the exception of one case with weak expression. Hematoxylin and eosin-stained sections from each sample were evaluated by an attending clinical pathologist with expertise in breast tissue to identify representative areas of tumor and tumor-associated stroma, as well as histologically normal breast epithelium and stroma. Laser capture microdissected RNA was isolated from samples and hybridized to Agilent Technologies SurePrint G3 Human GE 8 × 60K Microarrays (See Supplementary Information).

R/Bioconductor (vers 3.20; Bioconductor 3.1; ref. 21) was used for most analyses. Normalization was performed using the limma package (22) where loess was applied for dye bias correction, and quantile normalization was used across arrays. Replicates of noncontrol probes were aggregated by taking their mean value. To investigate technical error introduced in the LCM/microarray procedure, the stromal and matched adjacent normal sample from a single patient were repeated and found to be highly concordant (>0.8). Replicate expression profiles were then averaged for the remainder of the analysis. The most variable probe was chosen when there were multiple probes for the same transcript. Clear distinctions between normal and tumor stromal samples were also observed (See Supplementary Information; Supplementary Fig. S1). Raw and normalized microarray data have been deposited in the Gene Expression Omnibus database under accession number GSE90505.

Probes with an interquartile range > 2.0 across all samples were used as features in hierarchical clustering (Ward's algorithm, Pearson correlation distance). Samples were mean-centered and scaled for each transcript across all patients prior to clustering. Pvclust (version 1.3-2; ref. 23), was used to measure cluster stability with 100,000 iterations and we selected clusters containing at least 12 genes with an Approximately Unbiased (AU) value of >85%. Gene clusters deemed unstable by Pvclust were not selected for further analysis, and no further characterization was performed before selection.

This unbiased approach, which is described in ref. 24, ranks samples based on their expression of a specific gene set. In our case, we use the characteristic genes for each stromal and Lehmann and colleagues axis. This method estimates each patient sample as either low, intermediate, or high using a random resampling technique with 1,000,000 iterations.

To identify differentially expressed (DE) genes for each stromal axis, we fitted a linear model comparing the levels for each stromal axis using the R package limma (22) and corrected with Benjamini–Hochberg (P < 0.05). DE gene lists were examined using QIAGEN's Ingenuity Pathway Analysis (IPA, QIAGEN; https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/) and compared against the Molecular Signatures Database (MSigDB v5.2; http://software.broadinstitute.org/gsea/msigdb) for pathway analysis.

We used our compendium of 5,901 bulk expression profiles of invasive breast cancer samples from 13 nonoverlapping datasets generated on different technologies (5). We define poor outcome as an observed distant metastasis within 5 years of diagnosis (where available) and used ER and HER2 status as reported for each dataset where available. As many datasets lacked information on PR status, we used ER and HER2 negativity to define TNBC patients. All patients within the compendium were also labeled with intrinsic subtype values via PAM50 (25), and TNBC patients were labeled by TNBCType via the web-based tool (26).

To assign scores along the stromal axes to TNBC patients across the compendium of bulk expression profiles, we performed ROI₉₅ using the lists of differentially expressed genes for each stromal axis. Our software to generate stromal axis scores in a sample is available as a Bioconductor package entitled STROMA4 (http://bioconductor.org/packages/devel/bioc/html/STROMA4.html). STROMA4 was applied to each (of the 13) datasets of our TNBC compendium independently and the three classes (low, intermediate, high) were combined across all of these datasets. The Lehmann axes were similarly assigned using the characteristic genelist for each Lehmann subtype from the original publication. Assignments by ROI₉₅ were compared to assignments by the TNBCType web-based tool (26). Only four disagreements for “high” classification by both tools were observed (Supplementary Table S1) confirming the accuracy of the ROI₉₅ assignments.

Cohen kappa statistic was used to measure agreement between two ROI₉₅-based categorizations into low, intermediate, and high (fmsb package vers. 0.5.1). Enrichment analyses were performed via a one-sided Fisher exact test in R. The minimum P value between each one-sided test was used to determine whether there was significant enrichment or depletion. For variables with more than two levels, multiple tests were performed to determine enrichment for each level individually against all other levels.

For comparisons between a ternary variable and a binary variable (for example, lymph node status versus high, intermediate, low stromal axis score), we removed the intermediate category and then used Cohen kappa for the two binary variables. To determine association with distant metastasis free-survival, we used a Cox proportional hazards regression model via the coxph function in R (vers. 2.38) using the three (ordinal) levels estimated by the ROI₉₅.

To determine whether the fraction of observed combinations of stromal axis scores was higher than expected, we used a one-sided binomial test (stats package in R; ref. 21). The observed number of patients with the combination being tested was used for the number of successes, the total number of patients was used for the number of trials, and the hypothesized probability of success was determined as the product of the fractions for the individual axis scores being tested as observed in the TNBC patient compendium. A P value less than 0.05 was deemed to be significant.

To investigate stromal heterogeneity across TNBC tumors, 57 patient samples were selected based on negative ER, PR, and HER2 status according to clinical-pathological reports (Supplementary Table S2). Tumor nonepithelial (stromal) compartments were separately isolated by laser capture microdissection (LCM) and subjected to microarray-based gene expression profiling (see Supplementary Methods; refs. 18, 19, 27). Matched histologically normal stromal tissue was isolated for a subset of cases (n = 11). Following quality control and data normalization (see Materials and Methods), the success of the LCM procedure in isolating distinct tissue compartments was confirmed (see Supplementary Methods).

To establish whether differences are observed in TNBC stroma, the most variable genes in TNBC tumor stromal samples [interquartile range (IQR) > 2, n = 211 genes] were subjected to hierarchical clustering (Ward algorithm, Pearson correlation distance). Four distinct clusters were observed that contained a significant number of genes with strong pairwise gene-gene correlations of expression (Fig. 1A, colors along rows). These clusters were statistically stable and reproducible (pvclust, AU ≥ 85%). Genes within each cluster that exhibit strong coexpression across the patient cohort are termed the characteristic gene set.

To assign values along each stromal axis to TNBC tumors, patients were linearly ordered on the basis of the sum of ranks of the characteristic genes for each stromal axis signature independently (see Materials and Methods). A rank-based permutation test (ROI₉₅; ref. 24), was applied to each linear ordering to estimate boundaries of regions that delineate samples that are low, intermediate, or high for the characteristic gene set (Fig. 1B, black bars). Hence, each patient sample is independently scored for each of the four ternary axes (low, medium, high). This approach differs from traditional subtyping approaches that partition the patient cohort into distinct, nonoverlapping subtypes (Fig. 1C).

To characterize the molecular pathways and presence of specific cell types in each stromal axis, we identified differentially expressed genes between patients deemed low versus those deemed high for each stromal axis (LIMMA, FDR-adjusted P < 0.05 after ROI₉₅, Supplementary Fig. S2; Supplementary Table S3). For the first axis, genes differentially expressed include both general (CD2, CD3D, IL2Rα, IL2Rβ, IL2Rγ), and lineage-specific (CD4, CD8A, CD8B) T-cell–associated markers, as well as markers of a Th1-mediated antitumor response including IL15 (28), granzymes (GZMBA, GZMB, GZMK, GZMH; ref. 29), markers of an IFN response (IFI30, IFIT5; ref. 30), transcription factors involved in Th1 differentiation (STAT1, STAT4; refs. 31, 32), and TNFα-induced genes (TNFAIP2, TNFAIP8; refs. 33, 34). These genes had greatest expression in patients with high scores for this axis (purple, Fig. 1B and C). Pathway analysis for the first axis identified signatures linked to the proliferation of T lymphocytes and activation of cytotoxic T cells, confirming the associations with T cells (purple, Fig. 1B and C; see also Materials and Methods and Supplementary Table S4).

For the second axis, genes differentially expressed between low and high classes (magenta, Fig. 1B and C) include B-cell markers (CD19, CD79A, CD72), immunoglobulins (IGLL5, IGLL1, IGJ), and transcription factors associated with B-cell activation (POU2AF1, XBP1). These genes have greatest expression in samples scoring as high along this axis. Pathway analysis for the second axis identified signatures of B-cell proliferation (magenta, Fig. 1B and C; see also Materials and Methods and Supplementary Table S4).

For the third axis, differentially expressed genes (teal, Fig. 1B and C) include keratins (KRT6B and KRT23), and metallothioneins. These genes are expressed by tumor epithelial cells (35) and thus may represent invasive tumor cells that have retained some of their epithelial characteristics due to tumor plasticity (36).

For the fourth axis, differentially expressed genes (orange, Fig. 1B and C) include multiple collagens (collagens 1A1/2, 3A1, 5A1/2, 8A1/2, 10A1, 12A1, 16A1), PDGFRB, FAP, in addition to collagen stabilizing and modifying enzymes (P4HA2, MMP2, LOXL1). All of these are factors associated with a desmoplastic reaction (37, 38). Analysis of the fourth axis identified a signature of desmoid-type fibromatosis (hypergeometric test, P < 0.05; orange, Fig. 1B and C; see also Materials and Methods and Supplementary Table S4; ref. 39).

On the basis of these observations, we labeled the four stromal axes as T (T-cell, purple), B (B-cell, magenta), E (invasive epithelial cells, teal), and D (desmoplastic reaction, orange), respectively.

To test associations between the stromal axis scores and clinical variables (e.g., tumor size, grade, stage, outcome, etc.) requires a larger cohort of TNBC patient samples. Because of the unavailability of TNBC stromal datasets, we developed and tested a statistical method to estimate the status along each stromal axis in bulk expression data (Fig. 1D; Supplementary Table S5). This method, entitled STROMA4, was applied to a large cohort of TNBC patient samples (n = 1,098) selected from 13 individual nonoverlapping publicly available breast cancer datasets (see Materials and Methods; Supplementary Methods; ref. 5). Stromal axis assignments were computed independently per dataset, and pooled across the constituent datasets (Supplementary Table S6). This enabled us to test if the low, intermediate, and high partitions of each axis stratified patients by clinical outcome. While the D score (orange) did not demonstrate significant association with outcome, the T, B, and E scores (purple, magenta, teal) were significantly correlated with outcome [log-rank test, distant metastasis free survival (DMFS) at 5 years all P < 0.05; Fig. 1E]. This demonstrates that scores along the T, B, and E stromal axes inform on clinical outcome for TNBC patients.

To establish whether the stromal axis scores are associated with subtypes of TNBC derived from bulk tumor gene expression profiles (TNBCType; ref. 13), we first applied our approach to the TNBCType subtypes. Using the data from Lehmann and colleagues (13), the gene sets that underlie each of the six TNBC subtypes were subjected to our methodology, estimating their activation as either low, intermediate, or high across the TNBC compendium (Materials and Methods). This procedure generates six “Lehmann axes” in a format that allows comparison with the four stromal axes.

To examine if any of the Lehmann axes interact (e.g., whether samples high for one subtype are also high for the second), associations between all possible states of all possible pairs of Lehmann axis scores were determined. Subtyping schemes partition patient samples into disjoint groups with the assumption that if a patient belongs to one subtype (for which they have the appropriate molecular profile), they do not belong to another subtypes (their profile is sufficiently distinct). However, we observe strong statistical correlations between 11 out of 15 pairs of TNBCType axes (Fig. 2A, Kappa test, all P values <0.01). While the mesenchymal (M) and immunomodulatory (IM) scores display a near perfect (P < 1e−10) anticorrelation that is consistent with distinct subtypes, we also observe an equally strong (anti-)correlation between the basal-like-2 (BL2), luminal androgen receptor (LAR), and mesenchymal stem-like (MSL) scores. Hence, for example, even though by our method many patients are high for both BL2 and LAR, each will be assigned to only one of these subtypes under the TNBCType approach. This correlation does not support their identification as distinct subtypes.

To establish whether the four stromal axis scores show similar associations, we repeated this analysis with the stromal axes. Of the 6 pairs tested, only the T and B stromal axis scores were strongly correlated (P < 1e−10), while the D score showed a weak anticorrelation (P < 0.01) with the T and E scores (Fig. 2B).

To investigate associations between the Lehmann axis scores and our TNBC stromal axis scores, we performed a similar analysis comparing the Lehmann and stromal scores (Kappa test, Fig. 2C–F). Notably, the T stromal axis score (Fig. 2C), and to a lesser extent the B score (Fig. 2D), capture the inversely-correlated M and IM Lehmann axis scores (P < 1e−10). This observation is in line with the observation by Lehmann and colleagues that the IM score is correlated with lymphocytic infiltration (40). The stromal E score exhibits strong correlation with the BL1 score and anticorrelation with the LAR score (P < 1e−10; Fig. 2E). Patient samples estimated high for the D axis are almost always estimated high for the BL2, LAR, and MSL axes and low for the BL1 axis (P < 1e−10, Fig. 2F), which is in line with the observation that the MSL subtype is correlated to the presence of CAFs in the tumor (40). Together these observations highlight that the “Lehmann axes” are strongly associated with the stromal axes, and suggests that TNBC heterogeneity can be succinctly summarized by three distinct scores related to immune infiltration (B and T), androgen receptor signaling (inversely related to E), and a desmoplastic stroma (D).

High expression of the MSL and BL1 Lehmann axes was observed to be negatively and positively correlated with the proliferative index of the tumor, respectively (13), and in the preceding subsection we show that our D stromal axis is also strongly positively and negatively associated with the MSL and BL1 Lehmann axes, respectively. This suggests that the D axis score reflects the stromal response to the proliferation of neighboring tumor epithelial cells.

We investigated this relationship in three ways. First, samples with high scores along the D axis had a significantly lower percentage of Ki67-positive tumor cells than samples with low D scores (two-sided t test, P < 0.05, Fig. 3A). Second, using a gene signature of proliferation (41) and the ROI₉₅ to estimate proliferative states, we observed a strong statistical association between scoring along the D axis and expression of the proliferative signature (Fig. 3B and C; Kappa test, P < 0.01). Third, we observed that samples scoring low along the D axis are associated with higher grade (Fig. 3D). High-grade tumors tend to have a higher mitotic index (42). Together, these findings indicate that the D axis score is the stromal image of tumor proliferation from the neighboring epithelial cells.

The many alternative subtyping schemes for breast cancer propose varying numbers of patient partitions that range from just four subtypes via classic approaches based on ER, PR and HER2 status to 10 subtypes identified by IntClust through joint DNA and mRNA analysis (43). Similarly, there have been different numbers of “sub-sub”-types proposed for TNBC (13, 14). Our four ternary stromal axes suggests that as many as 3⁴(=81) such sub-subtypes could exist; each such subtype is described as a combination of the four axis scores (eg T-high, B-high, E-low, D-int). However, using our compendium of approximately 1,000 TNBC profiles, we observed that the number of patients assigned to each of the 81 possible combination subtypes varied significantly, with some subtypes populated by many samples and others with very few.

To measure this statistically, the enrichment or depletion of subtype populations were assessed relative to a background model (binomial-based test, see Materials and Methods). Only 15 of the 81 subtypes are significantly enriched beyond levels that would be observed solely by chance across our compendium (Supplementary Table S7). Many of the 15 enriched subtypes were either of the form T-high, B-high/B-intermediate, or of the form T-low, B-low/B-intermediate, suggesting a strong interaction between the T and B axis scores. Conversely, 17 of the 81 subtypes were significantly depleted, suggesting selection against some specific combinations of stromal axis scores (for example, the T-high, B-low/intermediate, D-high, E-high combination), in turn suggesting a complex interaction between these axes. The remaining 49 subtypes were populated as would be expected under a null binomial model.

Above, we have investigated the prognostic capacity (DMFS at 5 years) of the stromal axis scores individually, of which the B, T, and E axes were statistically significant. To establish whether the interactions between the stromal axes have prognostic capacity in the TNBC compendium, we focused on the 15 stromal subtypes that had larger than expected populations. We observe that all subtypes where the D score is high were not significantly different in terms of the prognostic capacity (Fig. 4A, log-rank test, P > 0.05). However, when the D score is low or intermediate, the spectrum of subtypes (varying states of B, T and E) do have significantly different survival characteristics (Fig. 4B; log-rank test, P < 0.05). Additional survival analysis across the full TNBC patient cohort supported these findings (Supplementary Methods; Supplementary Fig. S3).

To determine whether this effect was dependent on chemotherapy treatment, patients were stratified on the basis of whether they received chemotherapy and survival analysis was performed (Supplementary Figs. S4 and S5). As observed in the unstratified cohort, the D axis score acts as a master regulator of the other axes; none of the B, T, or E axes are prognostic in either chemotherapy-treated or -untreated D-high patients. However the B and T axes are prognostic among D-low/intermediate patients in chemotherapy-treated patients, while the B and E axes are associated with prognosis among D-low/intermediate patients not treated with chemotherapy. Thus, the D axis confounds the prognostic ability of the other stromal axes irrespective of adjuvant chemotherapy status.

The vast majority of previously reported prognostic gene signatures for breast cancer are derived from bulk expression data from the tumor proper. These signatures measure a broad range of tumoral hallmarks and cancer-related processes (e.g., proliferation, genomic instability, immune response). In a previous effort (5), we identified a subset of patients whose observed outcome was consistently mispredicted by almost all reported prognostic gene signatures. As our stromal D axis score appears to be a master controller of the prognostic capacity of the T, B, and E axes, we hypothesized that patient samples estimated high for the D axis might have higher inherent prognostic difficulty, that is, gene signatures that predict prognosis will almost always incorrectly predict D-high patients. If this is true, then knowledge of the state of the D axis might provide a significant breakthrough that would increase the accuracy of prognostic classifiers.

Using the inherent difficulty score from Tofigh and colleagues (5), we observed a significant difference in inherent difficulty of observed poor-outcome patients contingent on D-high versus D-low or D-intermediate status (two sample, two-sided t test; P < 0.05; Fig. 4C, DMFS at 5 years with blue and red representing good and poor outcome, respectively). Thus, those poor-outcome patients that are high for the D axis are systematically mispredicted as good-outcome.

Although the stromal axes were identified within a TNBC cohort, a natural question is to ask whether they have prognostic capacity in other breast cancer subtypes. It is well-established that signatures related to cell cycle and tumor proliferation provide prognostic information especially within ER-positive related cohorts (5, 10, 44), suggesting potential efficacy for the D stromal axis. The prognostic capacity of immune-related signatures are also well established, particularly within ER-negative cohorts, suggesting efficacy for the T and B axes (5, 10, 45, 46). We asked whether there is evidence that the microenvironmental states captured by our TNBC stromal axes can be used universally across the disease to predict disease progression.

Following the methodology of Tofigh and colleagues (5), we trained prognostic predictors for each of the stromal axes and compared these predictors with previously reported predictors (n ∼ 120), trained in the same manner for each breast cancer subtype, across a large collection (n ∼ 5,000) of invasive breast cancer profiles (Fig. 4D; Supplementary Fig. S6). More specifically, for the gene set of each stromal axis and each prognostic signature, a Naive Bayes' Classifier was generated under statistical cross-validation within datasets while reserving several complete datasets for additional independent validation (Supplementary Fig. S7). The use of such classifiers allows, for example, the learning procedure to “weight” specific genes within a gene signature as more or less important in predicting patient outcome. DMFS at 5 years was used as the clinical end point, and performance was measured by the product of accuracy, via standard survival analyses (log-rank test) and via a random sampling–based approach.

We observed that the stromal axis scores, despite being discovered among TNBC patients, were associated with patient prognosis in other patient cohorts. The predictors derived from the stromal axes are observed to cluster closely with other predictors polling similar processes (Fig. 4D; Supplementary Fig. S6). Specifically the D predictor clusters with a predictor derived from comparing normal and tumor stroma (47), the T predictor clusters with other signatures linked to the presence of immune cells in the tumor (48), and the B predictor clusters amongst other B-cell–related signatures (49). Interestingly, the E predictor clusters among predictors of metastatic potential [including metastasis to brain (50) and lung (51)] and among signatures related to immune suppression/protumor immune signaling (52, 53). This suggests that the E axis may represent an early signature in the microenvironment of invasion or metastases. Future work could determine whether the molecular events underlying the E axis could be therapeutically targeted to prevent metastases.

As perhaps expected, the T predictor was one of the best predictors in the TNBC cohort (Fig. 4D). However, it was also significantly associated with prognosis in all subtypes except the ER-positive/HER2-positive cohort. While the D predictor was also significant in these cohorts, it was among the top predictors of prognosis in unstratified pan-breast cancer analysis. Although the D predictor is significant in the TNBC subtype (P < 0.001), it was not among the highest-performing classifiers. It has been previously observed that signatures measuring proliferation have their best performance in unstratified analyses, likely because there are large concomitant differences in the proliferative indices and rate of poor outcome between ER⁺ and ER⁻ patients (5). Given that the D axis score correlates strongly with tumor proliferation, it is perhaps unsurprising that the D predictor is one of the best predictors in unstratified analysis. The B predictor is significantly associated with prognosis in the same cohorts as the T predictor, although it never appears among the highest performing signatures. Interestingly the E predictor showed a stronger association with prognosis in unstratified and ER-positive cohorts, than within the TNBC subtype. Together these results suggest that the properties underlying stromal axis scores are nearly universal to breast cancer, and have different prognostic capacities in different subtypes.

The poor-outcome TNBC subtype remains the focus of much research as we still lack effective classification markers, prognostic signatures, and targeted therapies. This study represents the first large-scale effort to investigate the tumor microenvironment in TNBC patients. With respect to the TNBC subtype, previous studies have focused on gene expression profiling (12–14) or DNA sequencing (54) of bulk material enriched for tumor cells. Efforts to study the tumor microenvironment, including our own (18–20, 27), have used LCM to isolate stromal elements in a pan-breast cancer fashion, not restricted to TNBC.

Recently the spatial and temporal heterogeneity of tumor stroma has received much focus. It has previously been reported that the stromal profile of a single section is sufficient to capture the spatial heterogeneity of the tumor stroma (55). In contrast, it has been observed that the response of the tumor stroma to neoadjuvant treatment is varied between patients, and that the stromal response to treatment is predictive of patient outcome (56). However, the heterogeneity of the stroma within a subtype has not been investigated.

Here, we identify four stromal axes in TNBC patients. A key distinction between this study and previous work is that each patient is scored along all axes, rather than being assigned to individual distinct subtypes. We observe that this method better captures the heterogeneity of the TNBC stroma. Despite being discovered in LCM-derived material, these stromal axes are shown to hold true even when applied to matching bulk expression sample profiles (Fig. 1D; Supplementary Table S5). This allows us to assign levels of the four axes to a large bulk-derived compendium of TNBC patients, revealing that the activation state of the B, T, and E axes were associated with patient outcome, and that the D axis score is associated with tumor proliferation.

Using the infrastructure and concepts from Tofigh and colleagues (5), we show evidence that the vast majority of existing gene signatures for predicting patient prognosis fail for many D-high patient samples. In particular, these signatures often incorrectly predict D-high poor-outcome patients to have good outcome. D-high patients are the least proliferative among the TNBC, although all of these tumors are highly proliferative in comparison with non-TNBC tumors. This suggests that the approximately 40% of TNBC tumors estimated to be D-high and the least proliferative have also been problematic for existing prognostic signatures. Therefore, when our novel D signature assigns a high score to a sample, current prognostic predictors should not be utilized as they will likely fail. Conversely, it is only when a patient sample is deemed low (or intermediate) for D that the T, B, and E axis scores provide additional prognostic information. These observations allow us to generate a subtyping schema that ablates the complexity of having many (81) potential subtypes (Fig. 4E). The variable survival rates observed among D-low patients also explains why the D score did not have prognostic capacity when studied in isolation (Fig. 1E).

When the D score is high, the immune response within the microenvironment (estimated via the T and B scores), and the E score are insufficient to predict outcome. The desmoplastic stroma as encompassed by the D signature may suppress the tumor-antagonistic effects of a stimulated immune response. This is consistent with the observation that the cohort of samples with high D scores have moderate to poor prognosis when compared to the entire TNBC cohort. Hence, the prior or concurrent targeting of desmoplastic stroma may enhance the therapeutic benefit of immunomodulatory therapy in this patient cohort.

Using bulk expression profiles, the TNBCType scheme estimates that six subtypes capture the heterogeneity of TNBC patients. By applying our methodology to the genes that define their subtypes, we present strong evidence that essentially every patient sample belongs to multiple Lehmann and colleagues subtypes. However, as their methodology assigns each patient to exactly one subtype (e.g., MSL and not LAR), these patients would be treated according to the standard of care for the MSL subtype, and potential antiandrogen therapies suitable for the LAR subtype would be ignored despite the evidence that the patient is also LAR positive. Therefore application of our approach may improve the already positive findings from Lehmann and colleagues (13, 15).

Ng and colleagues (10) suggest that TNBC may be characterized by three or possibly four signatures related to androgen signaling, immune infiltration, and a desmoplastic stroma, all of which are captured by our four stromal axes presented here. This effort is the first to offer a quantitative estimate of the number of distinct, populated TNBC subtypes. If there are four core axes of TNBC and each axis is scored as high, intermediate, or low, then there are a total of 3⁴(=81) possible subtypes. We found that 15 of these 81 subtypes had more patients than expected by chance, while 17 subtypes that had significantly fewer patients than expected by chance. The remaining 49 subtypes had populations within our datasets that are not significantly different than what we would expect by chance, and therefore larger TNBC cohorts are necessary to investigate their prognostic and predictive values. As the axes are associated with genomic subtypes, we note that development of an approach similar to AIMS (57) may be required in future.

Our approach has therefore identified four interacting stromal axes that can be combined to subtype TNBC patients. These axis scores have also displayed broader applicability among other cohorts of breast cancer patients, indicating that they play a key role in determining breast cancer outcome. The finding that immune-related (T and B) scores lose prognostic value in D-high cases is of special relevance in triple-negative breast cancer, where various immunotherapy approaches are entering clinical practice. The presence of a cohort comprising approximately 40% of TNBC cases, for which immune status is unrelated to prognosis, may underlie the lack of response to immune-directed agents observed in a majority of patients.

No potential conflicts of interest were disclosed.

Conception and design: A. Omeroglu, M.T. Hallett, M. Park

Development of methodology: S.M.I. Saleh, H. Zhao, A. Omeroglu, M.T. Hallett, M. Park

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): N.Bertos, M. Souleimanova, H. Zhao, A. Omeroglu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.M.I. Saleh, T. Gruosso, M. Gigoux, M.T. Hallett

Writing, review, and/or revision of the manuscript: S.M.I. Saleh, N. Bertos, T. Gruosso, M. Gigoux, A. Omeroglu, M.T. Hallett, M. Park

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): N. Bertos, M. Souleimanova, H. Zhao, A. Omeroglu, M.T. Hallett

Study supervision: M.T. Hallett, M. Park

Other (generated figures): S.M.I. Saleh

We thank Vanessa Dumeaux and Daniel Del Balso for their advice during the preparation of this manuscript and for the careful reading of this manuscript. We thank Dr. Marie-Christine Guiot for her assistance with anti-Ki67 immunohistochemistry.

This study was supported by funding from CQDM (Consortium québécois sur la découverte du medicament/Quebec Consortium for Drug Discovery), and the NIH to M. Park and M. Hallett. S.M.I. Saleh is supported by the CIHR Systems Biology Training program and a McGill Faculty of Medicine Internal Studentship (awarded as Gershman Memorial and Hugh E. Burke Fellowships). The breast tissue and data bank at McGill University is supported by funding from the Database and Tissue Bank Axis of the Réseau de Recherche en Cancer of the Fonds de Recherche du Québec-Santé and the Quebec Breast Cancer Foundation (M. Park). The bulk expression profile dataset was generated thanks to financial support provided by Genome Quebec (M. Hallett and M. Park).

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.