Multi-Omics Profiling Reveals Distinct Microenvironment Characterization and Suggests Immune Escape Mechanisms of Triple-Negative Breast Cancer

Triple-negative breast cancer (TNBC) is the most aggressive breast cancer subtype, which is defined as no expression of estrogen receptor (ER) and progesterone receptor (PR) and no amplification or overexpression of HER2 (1, 2). Early relapse and lack of therapeutic targets are major problems in TNBC treatment (3, 4). Recent research has demonstrated that TNBC had higher immunogenicity than other subtypes do, suggesting immunotherapeutic strategies for patients with TNBC (5–9). However, recent clinical trials indicated that immunotherapy, such as immune checkpoint inhibitors (ICIs), showed low efficacy in the whole population of patients with TNBC (10–12). A lack of patient selection based on tumor microenvironment landscapes might be the major reason for these disappointing results. Meanwhile, the whole landscape of TNBC microenvironment phenotypes remains unknown. Previous studies have focused on only one or two microenvironment cell subtypes of TNBC (13–20), which might result in the biased understanding of TNBC microenvironment. As microenvironment cells have intensive cross-talk, it is more rational to consider them all as a whole.

The development of next-generation sequencing provides opportunities to systemically explore the tumor microenvironment. As RNA sequencing of tumor tissue usually contains microenvironment cells, researchers have developed some expression profile–based estimation of the abundance of microenvironment cells in tumor tissue (21–24). Furthermore, researchers have analyzed the connection between DNA-level alternations and immune filtration to explore genomic alterations that drive the low immune infiltration (25–27). However, owing to small sample size and the lack of multi-omics data, few studies conducted genomic analysis of TNBC from the perspective of immunology.

Collectively, we questioned whether TNBC has heterogeneous microenvironment phenotypes and what genomic events drive the formation of these phenotypes. With multi-omics data for the largest single-center TNBC cohort, we successfully classified 386 TNBC samples into three microenvironment clusters with distinct potential immune escape mechanisms and genomic drivers.

We retrospectively selected 386 consecutive patients with TNBC who underwent surgeries at the Department of Breast Surgery, Fudan University, Shanghai Cancer Center (FUSCC; Shanghai, China), from January 1, 2007 to December 31, 2014. Detailed inclusion criteria for the 386 samples were as follows: (i) female patients; (ii) unilateral invasive ductal carcinoma; (iii) pathologic examination of the ER, PR, and HER2 status performed by the Department of Pathology at FUSCC through immunochemical analysis and in situ hybridization (for HER2 status only); (iv) patients with no evidence of metastasis at the time of diagnosis; and (v) sufficient frozen tissue for further research. In summary, RNA-seq data (tumor tissues: n = 245; paired normal tissues: n = 90), HTA 2.0 microarray data (n = 141), whole-exome sequencing data (n = 268), OncoScan microarray copy-number data (n = 335), hematoxylin and eosin (H & E) sections data (n = 300), and tissue microarray data (n = 181) were obtained. In addition, we invited two pathologists to evaluate the stromal and intratumoral tumor-infiltrating lymphocytes (sTIL and iTIL, respectively) and fibrosis and necrosis in H & E pathologic sections on the basis of established guidelines (28, 29). The studies were conducted in accordance with the Declaration of Helsinki. All the tissue samples included in this study were obtained with approval from the independent ethics committee/Institutional Review Board at FUSCC Ethical Committee, and each patient provided written informed consent. The follow-up of our cohort was completed on June 30, 2017. The median follow-up length was 47 months (interquartile range, 29.3–73.5 months). The events included in the relapse-free survival (RFS) analysis were defined as the first recurrence of locally, regionally, or distantly invasive disease, a diagnosis of contralateral breast cancer, or death from any cause. Patients without events were censored at the last follow-up.

Detailed information regarding the sample processing and sequencing data generation is provided in the Supplementary Materials and Methods.

To construct a compendium of microenvironment genes related to specific microenvironment cell subsets, we considered three aspects in the gene selection: first, the gene is specifically expressed in one specific microenvironment cell subset; second, other normal tissues do not express the gene; and third, the compendium contains major cell types of the TNBC microenvironment. After researching papers, we modified two gene signatures, CIBERSORT (21) and MCP-Counter (23), to construct our compendium. First, we filtered the gene list of CIBERSORT (see their supplement 2). A 1.5-fold change cutoff from the highest expression cell type to the next highest expression cell type was considered as a criterion for cell type–specific expression. As a result, 324 genes were selected. As CIBERSORT do not contain signatures of fibroblasts and endothelial cells, which were important component of TNBC microenvironment, we added extra 40 genes for these cells (32 for endothelial cells and eight for fibroblasts) from MCP-Counter gene list to our compendium. In all, our TNBC compendium contained 364 genes representing 24 microenvironment cell types (Supplementary Fig. S1; Supplementary Table S1). Subsequently, we used single sample gene set enrichment analysis (ssGSEA, “GSVA” function in R) to calculate the abundance of each cell subset in each sample with expression data. Similarly, we also referred to another published article (22) to construct immune cell signatures for calculating the abundance of types 1, 2, and 17 T helper cells (Th1, Th2, and Th17) and myeloid-derived suppressor cells (MDSC).

To describe the constituent pattern (i.e., relative cell proportions) of microenvironment cell subsets within clusters, we first used the tumor purity, which was calculated by ASCAT (30), to adjust the enrichment scores of each microenvironment cell subset. The adjusted enrichment score was calculated as the enrichment score divided by (1 − tumor purity). We then illustrated the enrichment score distributions of each microenvironment cell subset to represent the relative cell proportions (“density” function in R).

We performed k-means (“kmeans” function in R) clustering and Nbclust testing (“NbClust” function in R, index = “all”) to determine the optimal number of stable TNBC microenvironment subtypes. To cluster samples based on the constituent pattern of each microenvironment cell type, we scaled each sample before clustering. Silhouette analysis was performed to confirm the stability of the clustering. For heatmap plotting (“pheatmap” function in R), we utilized the k-means clustering result to reorder the samples and scaled the original ssGSEA results before plotting.

We developed the univariate and multivariate Cox proportional hazards model to analyze the prognostic significance of microenvironment phenotypes. Age, tumor size, the number of positive lymph nodes, and PAM50 subtypes were first analyzed in univariate Cox proportional hazards model. All significant variables were then included as covariates in the multivariate Cox proportional hazards model. We further validated the prognostic value of microenvironment clusters by comparing the predictive efficacy between two models: one consisting of tumor size and the number of positive lymph nodes as covariates; and the other adding microenvironment phenotypes as the third covariates. The AUC of time-dependent ROC curves (“timeROC” package in R) was set as the indicator of prognostic efficacy. Statistical tests of the difference in AUC between the two models were conducted in the first 3 years of follow-up (“compare” function in R). We also evaluated the prognostic value of each cell subsets in the whole cohort and within each microenvironment cluster. In each analysis process, the abundance of cells was divided into high or low categories according to the optimal cutoff (“cutp” function in R). A univariate Cox proportional hazards model was then performed for the categorical variables of cell abundance.

The detailed calculation of several bioinformatic indicators, such as neoantigens, cancer testis antigens (CTAs), homologous recombination deficiency (HRD) scores, and intratumoral heterogeneity (ITH), are described in Supplementary Materials and Methods.

We referred to a published article (31) to establish a signature of 10 oncogenic pathways containing 331 genes. We then used the ssGSEA method on these gene sets to generate enrichment scores for each pathway in each sample. The enrichment scores were calculated as the activated score minus the repressed score. Subsequently, we compared the ssGSEA score of each pathway among the three clusters.

Genes with mutation frequencies greater than 2.5% were included in our comparison. To assess cluster-specific mutated genes, clusters were modeled (using logistic regression) as a function of the gene's mutational status (ignoring silent events) and mutation load. The latter one was included to diminish confounding effects. P value less than 0.05 after adjusting for mutation load was considered significant. Because logistic regression model only receipt dependent variable with binomial distribution, we performed comparison in every two clusters. We conducted permutation test (“glmperm” function in R program) in logistic analysis. Cluster-specific mutations were defined as significantly mutated genes from the comparison of the given cluster with one of the other two clusters (Supplementary Table S6, column 13).

The ASCAT algorithm was used to adjust the copy number of genes based on ploidy and purity. Samples were not included in further analysis when ASCAT-evaluated ACF score equaled 1 or if ASCAT failed to evaluate the purity. A regression approach similar to mutation analysis was then used to test for copy-number variation associations. To test a given gene, clusters were modeled as a function of the gene's copy number (categorized by the median value) and somatic copy-number variation (SCNV) load. To explicitly define amplification and deletion events, we applied this linear regression approach twice. The first run was amplification centric, and thus, the negative copy-number values were set to zero. The second run was deletion centric, and thus, the positive copy-number values were set to zero. The cluster-specific amplified “peak” was defined as a continuous stretch of genes (arranged in genomic order) with a nominal P value less than 0.05 and log₂ ratio greater than log₂(2.5/2). Similarly, the cluster-specific deleted “peak” was defined as a continuous stretch of genes (arranged in genomic order) with a nominal P value less than 0.05 and a log₂ ratio less than log₂(1.5/2). The comparison methods mainly referred to some previous immunologic articles (25, 26) and the cut-off values, log₂(2.5/2) and log₂(1.5/2), were also obtained from one published high-quality article of this field (32). Cluster-specific peaks were defined as significant peaks in each comparison between the given cluster and every one of the other two clusters.

Student t test, Wilcoxon test, and Kruskal–Wallis test were utilized to compare continuous variables and ordered categorical variables, such as mutation load, neoantigen load, HRD score, CTAs number, and ITH. Prior to the comparisons, the normality of the distributions was tested with Shapiro–Wilk test before comparison. Pearson χ² test or Fisher exact test were employed for the comparison of unordered categorical variables. Permutation test was conducted in the comparison of gene mutation frequencies among clusters. Correlation matrices were created with Pearson or Spearman correlation. Survival analysis was performed using the Kaplan–Meier method, and the survival of the clusters was compared using the log-rank test. All the tests were two sided, and P < 0.05 was regarded as indicating significance, unless otherwise stated. The FDR correction was utilized in multiple tests to decrease false positive rates. All of the analyses were performed with R software (version 3.4.2, http://www.R-project.org).

All data can be viewed in The National Omics Data Encyclopedia (http://www.biosino.org/node) by pasting the accession (OEP000155) into the text search box or through the URL: http://www.biosino.org/node/project/detail/OEP000155. The sequencing data is also available in GSE118527 (OncoScan), GSE76250 (HTA 2.0), and SRP157974 (WES and RNAseq).

We first established a reference microenvironment compendium that included 364 genes representing 24 microenvironment cell subsets, to systematically characterize the microenvironment phenotypes of TNBC (Materials and Methods, Supplementary Fig. S1; Supplementary Table S1). We then estimated the abundance of 24 microenvironment cell subsets in each sample and confirmed the accuracy of our results by comparing them with those of other microenvironment signatures (Supplementary Table S2). Subsequently, we performed k-means clustering of the TNBC microenvironment phenotypes. All 386 TNBC microenvironment phenotypes were classified into three heterogeneous clusters (Fig. 1A). Our analysis revealed that three was the optimal and stable clustering number (Supplementary Fig. S2). Cluster 1, the “immune-desert” cluster (type 1 “cold tumor”), was characterized by relatively low microenvironment cells infiltration, and thus, we selected a “red light” for this cluster, which represented low immune infiltration. Cluster 2, the “innate immune-inactivated” cluster (type 2 “cold tumor”), was characterized by the infiltration of inactivated innate immune cells, fibroblasts, and endothelial cells, and thus, we chose a “yellow light” to indicate moderate immune infiltration. Cluster 3, the “immune-inflamed” cluster (“hot tumor”), was characterized by relatively high innate and adaptive immune cells infiltration, and we therefore selected a “green light” for this cluster to denote abundant immune infiltration. Moreover, we described the distributions of the abundance of microenvironment cell subsets to analyze the relative proportions of cell subsets within clusters. Adaptive immune cells and inactivated innate immune cells represented the major proportions in the microenvironments of the three clusters. The relative weights of innate immune cells and nonimmune cells were increased in cluster 2, whereas the relative proportion of adaptive immune cells was increased in cluster 3 (Fig. 1B).

To validate the expression profile–based clustering, we first evaluated the sTILs and iTILs and CD8⁺ cells in pathologic sections among the clusters. Cluster 3 had significantly higher sTILs (mean, 19.86 vs. 18.06 vs. 32.65; P < 0.001), iTILs (mean, 5.27 vs. 4.58 vs. 11.02; P < 0.001), and CD8⁺ cells (mean 13.36 vs. 12.61 vs. 25.02; P < 0.001) than the other two clusters (Fig. 1C). We further analyzed other immune signatures to validate our microenvironment clustering. The distributions of Th1, Th2, and Th17 and MDSCs were also coincident with our established microenvironment clusters (Fig. 1D). Moreover, we compared the distributions of the TNBC intrinsic subtypes among clusters. Cluster 1 mainly consisted of the classical basal-like and immune-suppressed (BLIS) subtype; cluster 2 was mainly composed of tumors with mesenchymal stem–like (MSL) features, which were mostly classified as the nonbasal subtype in the PAM50 clustering; and clusters 3 was primarily made up of the immunomodulatory subtype, a nonclassical basal-like subtype. The luminal androgen receptor TNBC subtype was dispersed among all three clusters (Fig. 1E). Furthermore, we utilized the expression profiles of the METABRIC and TCGA cohorts to externally validate the repeatability of our clustering result (Supplementary Figs. S3 and S4). Overall, we demonstrated that TNBC had three heterogeneous microenvironment phenotypes that could not be fully explained by mRNA-based TNBC intrinsic subtyping. The features of the three clusters are summarized in Fig. 1F.

Considering the important role of the tumor microenvironment in prognosis, we investigated the clinical relevance of the microenvironment clusters. The three clusters had similar clinicopathologic characteristics (Supplementary Table S3). Cluster 3 had significantly better RFS (log-rank P = 0.04) and overall survival (OS; log-rank P = 0.04) than the other two clusters (Fig. 2A). The multivariate Cox proportional hazards model also revealed that cluster 3 independently predicted better RFS in TNBC (HR, 0.45; 95% confidence interval, 0.21–0.97; P = 0.04; Fig. 2B). The time-dependent AUC demonstrated that the addition of microenvironment clusters into the Cox proportional hazards model significantly increased the prognostic efficacy of 1- (AUC, 0.72 vs. 0.81; P = 0.0037) and 2-year (AUC, 0.74 vs. 0.81; P < 0.001) recurrence (Fig. 2C). Furthermore, we explored the prognostic significance of each cell subset (Fig. 2D; Supplementary Table S4). In the whole cohort, a higher infiltration of immune cells, even the infiltration of immune suppressive cells, predicted better prognosis. However, a within clusters analysis revealed that the prognostic significance of the cell subsets was diverse, even opposite, among the three clusters. Similar results were exhibited in METABRIC cohort (Supplementary Fig. S3).

The heterogeneity of the TNBC microenvironment phenotypes led to the question of whether different clusters of TNBC had distinct tumor immune escape mechanisms. On the basis of the immunoediting theory of previous publications (33), we first researched the extrinsic immune escape mechanism. This concept indicates that microenvironment components other than tumor cells contribute to the immune escape of tumor cells. The extrinsic immune escape mechanism consists of four major aspects: lack of immune cells, presence of immunoinhibitory cells [such as type 2 macrophages and regulatory T cells (Tregs)], high concentrations of immunoinhibitory cytokines (such as IL10 and TGF-β), and fibrosis (34).

We compared the estimated number of microenvironment cells between the tumor site and a paired normal site (Fig. 3A). Cluster 1 had almost no more microenvironment cells in the tumor site than in the paired normal site, suggesting an inability to attract innate immune cells (resulting in failure of adaptive immunity). The comparison of the expression of STING, the factor of spontaneous initiation of innate immunity, and other molecules potentially involved in initiation of innate immunity among the three clusters also suggested this theory (Fig. 3C; Supplementary Figs. S3F, S4D, and S5; ref. 35). Cluster 2 had more resting and activated innate immune cells in the tumor site, suggesting chemotaxis but inactivation of innate immune cells (also resulting in failure of adaptive immunity). Cluster 3 had not only abundant active innate and adaptive immune cells, but also immunosuppressive cells, such as Tregs, in the tumor site, suggesting a role of immunosuppressive cells in immune escape. The expression of chemokines was consistent with these results (Fig. 3B). Clusters 2 and 3 had higher expression of chemokines, including CCL4, CXCL9, and CXCL10, which have been proved to attract dendritic cells (DC) and CD8⁺ T cells. In addition, cluster 3 had both higher secreted immunostimulatory and immunoinhibitory cytokines; cluster 2 had increased secretion of immunoinhibitory cytokines only; and these cytokines were all relatively low in cluster 1. Notably, the difference in the expression level of these cytokines was not derived from SCNVs (all P > 0.05; Fig. 3B). Furthermore, the fibrosis existed in all three clusters, but it did not differ significantly among them (P = 0.22; Fig. 3D and E). Overall, the inability to attract innate immune cells, the chemotaxis but inactivation of innate immunity, the increase of immunoinhibitory factors after immune stimulation might contribute to the extrinsic immune escape of the three clusters, respectively.

We further investigated the potential intrinsic immune escape mechanisms of TNBC. Intrinsic immune escape indicates that tumor cells directly mediate their own immune escape. There are at least two aspects of intrinsic immune escape: tumor immunogenicity and immune checkpoint molecules expression (33).

We first compared some potential factors determining tumor immunogenicity among the three clusters: mutation load (Fig. 4A), neoantigen load (Fig. 4B; Supplementary Fig. S6), chromosomal instability level (Fig. 4C; Supplementary Fig. S7), CTA level (Fig. 4D; Supplementary Fig. S8), necrosis level (Fig. 4E), ITH (Fig. 4F) and tumor antigen-presenting capability (Fig. 4G, left). The first five factors were the main source of tumor antigens. In general, the difference in the tumor antigen burden among the three clusters of TNBC was not as large as that between microsatellite-instable and microsatellite-stable colorectal cancers. Cluster 1 (highest mutation load, HRD score, CTA burden, and necrosis level; all P < 0.05) and cluster 2 (lowest mutation load, HRD score, and CTA load; all P < 0.05) of TNBC had relatively high and low tumor antigen burdens, respectively. Moreover, both the clusters had lower MHC I–related antigen-presenting molecules expression than cluster 3 (all P < 0.001), contributing to their low immunogenicity. The difference in the expression of antigen-presenting molecules among clusters was also not explained by SCNVs (all P > 0.05; Fig. 4G, bottom left). In addition, the three clusters had no significant differences in ITH, although cluster 3 showed a nonsignificant tendency toward having lower ITH (mean fraction of subclonal mutations, 20.27% vs. 21.26% vs. 15.97%; P = 0.096; Fig. 4F). Overall, the difference in tumor immunogenicity among the TNBC clusters might be relatively small, with clusters 1 and 2 having relatively low immunogenicity.

The expression of immune checkpoint molecules after immune stimulation is another potential important intrinsic immune escape mechanism. Therefore, we referred to a database of costimulatory and coinhibitory molecules (https://www.rndsystems.com/cn/research-area/co–stimulatory-and-co–inhibitory-molecules) to compare these immunomodulators among clusters. We demonstrated that cluster 3 had a higher expression of costimulating (most P < 0.05) and immune checkpoint molecules (all P < 0.05) than the other clusters (Fig. 4G, right). This result suggested that cluster 3 tumors expressed immune checkpoint molecules to avoid immune killing after immune stimulation. The high expression of one costimulator in cluster 3, TNFSF8 (mean log₂ fold change relative to the paired normal tissue, 0.05 vs. 1.41 vs. 1.70; P < 0.001), could be explained by SCNVs (mean log₂ ratio of somatic copy numbers, −0.25 vs. −0.14 vs. 0.03; P = 0.002; Fig. 4I). In addition, we noticed that a coinhibitor, VTCN1 (B7-H7), exhibited higher expression in cluster 1 (mean log₂ fold change relative to the paired normal tissue, 0.57 vs. −1.10 vs. −0.67; P = 0.014) and was negatively correlated with immune infiltration (Spearman correlation = −0.27; P = 0.039; Fig. 4H). Furthermore, we investigated the relationship among immune infiltration [TILs and cytolytic activity (CYT); ref. 25], immunogenicity, and expression of immune checkpoint molecules. We demonstrated that immune infiltration and the expression of most checkpoint molecules were positively correlated, whereas the mutation load, neoantigen load, HRD score, necrosis, and ITH seemed not to be correlated with these factors (Fig. 4J).

We further investigated the genomic alterations that could be correlated with the low immune infiltration in clusters 1 and 2. We aimed to identify some potential targets to reverse the absence of adaptive immune infiltration in these clusters.

We first referred to published signatures (Supplementary Table S5) to calculate the enrichment scores of 10 common oncogenic pathways among the three clusters (31). The Hippo, MYC, PI3K and cell cycle–related pathways had higher scores in cluster 1 (all P < 0.01); the Notch, TGF-β, NRF2, and RTK/RAS pathways were enriched in cluster 2 (all P < 0.001); and the scores of the Wnt pathway were higher in cluster 3 (all P < 0.001; Fig. 5A). GSEA validated some of these results (Supplementary Fig. S9). When selecting cluster-specific mutated genes (Materials and Methods; Fig. 5B; Supplementary Table S6), we found that mutations among the PI3K-AKT pathway members [PI3KCA (30.8%), AKT1 (7.7%), PIK3R1 (6.1%), and PKD1 (6.1%)] were most frequently detected in cluster 2. After adjusting for the mutation load, PI3KCA and AKT1 mutations remained significantly increased in this cluster (all adjusted P < 0.05). Our cluster-specific SCNVs analysis (Materials and Methods; Fig. 5C; Supplementary Table S7) demonstrated that the amplification of 8q24.13–8q24.3 (MYC) was more frequent in cluster 1 (mean log₂ ratio of somatic copy numbers, 3.99 vs. 3.06 vs. 3.19; P_adj < 0.05) consistent with the results of the expression profile analysis. Other cluster 1–specific amplifications included 1p34.2 and 3q26.1-3q26.31. 1p36.22-1p36.21 (TNFRSF8, TNFRSF1B) was a cluster 1–specific deletion. Moreover, the amplification of 20p13-20p12.1 and 20q11.22 and the deletion of 9p24.3-9p22.3 (CD274, PDCD1LG2) were specific to cluster 2. We annotated genes located in copy-number peaks with the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathways. Notably, genes located in the cluster 2–specific amplification peaks were annotated to various pathways, including the innate immune response (Supplementary Table S8). Overall, our analysis revealed that some genomic alterations might drive the low immune infiltration in cluster 1 and cluster 2.

Using multi-omics data from the largest single-center TNBC cohort, our study revealed three heterogeneous TNBC microenvironment phenotypes and their clinical significance. Although cluster 1 and cluster 2 were both so-called cold tumors, they had distinct microenvironment phenotypes. Cluster 3 was the so-called hot tumor. We emphasized the following characteristics of the clusters that might contributed to immune escape: defects in the attraction of innate immune cells in cluster 1, chemotaxis but inactivation of innate immune cells and low tumor antigen burden in cluster 2, and high expression of immune checkpoint molecules in cluster 3. Furthermore, we found some genomic alterations that might drive the low immune infiltration of clusters 1 and 2 (Table 1). To the best of our knowledge, this study constitutes the first systemic analysis of the microenvironment heterogeneity of TNBC. The clustering results are in accordance with the immunologic principles described by previous articles (36).

Our study has important implications for clinical translations. First, our results might facilitate the selection of suitable patients for ICIs treatment. We revealed a cluster of “hot tumors” in TNBC (cluster 3, approximately 28% of all TNBCs) and demonstrated that high expression of immune checkpoint molecules might lead to the immune escape of this cluster. Although most clinical trials have shown that the efficacy of ICIs in TNBC was less than 10% (10, 11), these patients usually received several lines of treatment without detection of the real status of immune checkpoint proteins before immunotherapy. Notably, the efficacy of PD-L1 inhibitors in first-line monotherapy could be up to 25% (10), which was consistent with the percentage of “hot tumors” (28%) in our research. We hypothesized that some TNBC cells, which were sensitive to ICIs, were also more sensitive to chemotherapy. As a result, after several lines of chemotherapy, ICI-sensitive cells were eliminated. The level of immune checkpoint proteins after surgery could not represent this fact after several lines of chemotherapy. Considering the more targeted efficacy of ICIs than chemotherapy, we speculated that ICIs might need earlier lines application in “hot tumors” of TNBC. Second, we identified potential indicators to predict ICI efficacy in TNBC. Previous studies have suggested that the tumor mutation burden (TMB), checkpoint molecules expression, and TILs score might individually and collectively predict ICI efficacy in other tumor types (37). However, we illustrated that the TMB of TNBC was not associated with the TIL score or checkpoint molecules expression, whereas the latter two were highly correlated (Fig. 4J), which suggests that the latter two, but not the mutation load, might be ideal predictors. Recently, IMpassion130, the first phase III clinical trial of ICI in metastatic TNBC, also demonstrated that PD-L1 might be a suitable biomarker to predict ICI efficacy in TNBC (12). Moreover, the nonsignificant association between T-cell infiltration and TMB in TNBC has been reported in previous articles (25, 26). Several factors in addition to TMB, such as local tumor microenvironment, host immune state, and neoantigen quality, might influence the T-cell infiltration in the tumor microenvironment as well (38–40).

Our study could also help facilitate research on the relationship between oncogenic pathways and immune infiltration in TNBC. The transformation from a “cold tumor” to a “hot tumor” is currently a popular topic in cancer research. Recent studies have suggested that the activation of oncogenic pathways might decrease immune infiltration (41). However, few studies have focused on the mechanisms of impaired immune infiltration in TNBC. Our study demonstrated that at least two phenotypes of “cold tumors” exist in TNBC. Cluster 1 was characterized by almost no immune cell infiltration and a high percentage of MYC amplification. Previous studies have indicated that MYC amplification could induce the expression of CCL5, CCL23, IL1β, CD47, and PD-L1, inactivate DCs and macrophages; and limit natural killer, T, and B cells recruitment (42–44). Therefore, we speculated that MYC-induced low innate immune cells chemotaxis might be the reason for the poor immune infiltration in cluster 1. In addition, the features of cluster 2 were chemotaxis but inactivation of innate immunity and higher fibroblast and endothelial cells infiltration, with a higher mutation frequency in PI3K-AKT pathway members. The TGF-β–related and cancer-associated fibroblasts (CAF)-related pathways were enriched in cluster 2. Therefore, we speculated that hyper-activated PI3K-AKT pathway in tumor cells might play a major role in suppressing immune infiltration in cluster 2, such as through TGF-β secretion and CAF differentiation. In particular, previous studies have revealed that pan-PI3K inhibition enhances antitumor immunity and anti-PD1 response in breast cancer (45).

Our research had some limitations. First, there are more than 24 stromal cells types, and it is difficult to precisely include all phenotypes. We classified these cell subsets into three categories, namely, adaptive and activated innate immune cells, inactivated innate immune cells, and nonimmune cells, which solved this problem to some extent. In addition, immunogenomic analysis could not reflect the cause and consequence effect. Potential driver molecules in our research, such as MYC and the PI3K-AKT pathway, require further functional validation. Experimental studies are ongoing in our laboratory.

In conclusion, our study revealed that the microenvironment phenotypes of TNBC could be classified into three heterogeneous clusters with distinct potential immune escape mechanisms. Specific oncogenic pathways might drive the formation of these microenvironment phenotypes.

No potential conflicts of interest were disclosed.

The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the article.

Conception and design: W. Huang, K. Yu, Y.-Z. Jiang, Z.-M. Shao

Development of methodology: Y. Xiao, D. Ma, C. Suo, W.-T. Yang, Y.-Z. Jiang, Z.-M. Shao

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y. Xiao, D. Ma, J. Shi, M. Ruan, L. Shi, W. Huang, K. Yu, Y.-Z. Jiang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Xiao, D. Ma, S. Zhao, C. Suo, M.-Z. Xue, M. Ruan, H. Wang, J. Zhao, Q. Li, P. Wang, S. Huang, Y.-Z. Jiang

Writing, review, and/or revision of the manuscript: Y. Xiao, D. Ma, C. Suo, W. Huang, K. Yu, F. Bertucci, Y.-Z. Jiang, Z.-M. Shao

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): W. Huang, X. Hu, K. Yu, Z.-M. Shao

Study supervision: C. Suo, P. Wang, L. Shi, K. Yu, F. Bertucci, Y.-Z. Jiang, Z.-M. Shao

Other (binding prediction between peptides and MHC molecules based on bioinformatics): Q. Li

The authors wish to thank Prof. McGranahan for providing us with R code to calculate ITH score in TNBC. This work was supported by grants from the National Natural Science Foundation of China (81874112, 81874113, 81572583 and 81502278), the Training Plan of Excellent Talents in Shanghai Municipality Health System (2017YQ038), the “Chen Guang” project supported by Shanghai Municipal Education Commission and Shanghai Education Development Foundation (17CG01), Shanghai Pujiang Program (18PJD007), the Training Plan of Excellent Talents of Fudan University Shanghai Cancer Center (YJYQ201602), the Municipal Project for Developing Emerging and Frontier Technology in Shanghai Hospitals (SHDC12010116), the Cooperation Project of Conquering Major Diseases in Shanghai Municipality Health System (2013ZYJB0302), the Innovation Team of Ministry of Education (IRT1223), and the Shanghai Key Laboratory of Breast Cancer (12DZ2260100).

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.