Cellular heterogeneity map of diverse immune and stromal phenotypes within breast tumor microenvironment

Data curation and cohort characteristics

The RNA-seq data and clinical parameters from 1,092 patients with breast cancer were obtained from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). Of the 1,092 breast cancer patients included in this study, 112 paired normal tissues were identified and recurrent tumor tissues were excluded. The clinical characteristics of the cohort are listed in Table 1. 276 known Cancer/Testis (CT) genes were downloaded from CTDatabase (http://www.cta.lncc.br/).

Bioinformatics analysis

xCell (http://xcell.ucsf.edu/) is a high-resolution gene-signature-based method for cell type enrichment for up to 64 cell types, including immune and stroma cells. We used xCell R package (Aran, Hu & Butte, 2017a) (Beta version) from GitHub in R (version 3.3.1) to deconvolute the cellular heterogeneity within the breast tumor microenvironment from RNA sequencing data. We determined the cellular heterogeneity of 1,092 breast tumor tissues and 112 normal tissues using the xCell method. The 64 cell types were divided into four groups, including 34 immune cells, 13 stromal cells, 9 stem cells, and 8 other cells (Table 2). Over half of the 64 cell types were immune cells, providing a full view of the innate and adaptive immune status with detailed cell subtypes, including CD4+ naive cells, CD4+ T-cells, CD4+ Tcm, and CD4+ Tem. Stromal cells, including fibroblasts, osteoblasts, and pericytes were also included. ImmuneScore and StromaScore were generated by the xCell package using the sums of fractions of certain cell types (Aran, Hu & Butte, 2017b). Tests for differences and correlations were performed. We used the t-Distributed Stochastic Neighbor Embedding (t-SNE) method with tsne package (version 0.1-3) to perform cluster analysis based on cell fraction types.

Survival analysis

Univariate and multivariate COX regressions were performed using the survival package (version 3.1-7) to search for survival-associated genes. The best cutoff value for each factor was determined using the Survminer package (version 0.4.6). Significant prognostic factors were displayed in a forest plot and the most significant factors were further evaluated using multivariate analysis. The final prognostic model was built with five independent factors including CD8+ T cells, mesangial cells, NKT, keratinocytes, and class-switched memory B cells. We used the survivalROC package (Version 1.0.3) in R, which uses a time-dependent ROC curve estimation with censored data (Heagerty, Lumley & Pepe, 2000) to compare the aptitude of the individual prognostic factors. The final prognostic model was used to generate the area under the curve (AUC) of the receiver-operator characteristic (ROC) curve for each parameter.

Statistical analysis

Differentially enriched cell types between groups were compared using the Student’s t-test (two groups) or one-way ANOVA analysis (three groups). Correlation analyses were performed using the Spearman method. The survival curves were compared using the Kaplan–Meier method and log-rank test. All tests were two-sided and p < 0.05 was considered to be statistically significant unless otherwise noted. Data were analyzed using R (version 3.4.4).

Breast tumor tissues had higher fractions of immune cells than normal tissues

The median fractions for each cell type were calculated for normal and breast tumor tissues and the proportions of the 64 cells were found to differ between breast tumor and normal tissues (Fig. 1A, Table 2). Breast tumor tissue had higher fraction of immune cells with red to light blue markers whereas normal tissue had larger proportions of stem and stromal cells with blue to red markers (Fig. 1A). Unsupervised cluster analysis revealed that breast tumor tissues and the adjacent normal tissues were clustered into different groups. Immune cells were also clustered into several subgroups (Fig. 1B), indicating that the cellular heterogeneity in tumor vs. normal tissues was much greater than that in a single sample. Dimensionality reduction and visualization by t-Distributed Stochastic Neighbor Embedding (t-SNE) also suggested clear difference between the tumor and normal tissues (Fig. 1C).

We compared the fractions of each cell type between breast tumor and normal tissues, revealing dramatic differences in the number of cell types between the tumor and normal tissues (Fig. 1D–1W). Immune cells tended to be more diverse compared to stem or stromal cells. For innate immune cells, neutrophils were more prevalent in normal tissues whereas eosinophils were higher in tumor tissues (Fig. 1D–1E). There was not a significant difference in DC cells between normal and tumor tissues. However, iDC was significantly lower in tumor tissues but pDC and aDC were significantly higher (Fig. 1F–1I). This phenomenon was also seen in macrophages, in which macrophage M1 was higher in tumor tissues while macrophage M2 was lower (Fig. 1J–1K). CD4+ Tcm was found to be significantly lower in tumor tissues, while CD4+ Tem, CD8+ naïve T cells, and CD8+ Tcm were significantly higher (Fig. 1L–1O). Plasma cells, pro B cells, Tgd, Th1, Th2 cells, and Tregs were also found to be significantly higher in tumor tissues (Fig. 1P–1U). Representative stromal cells, such as endothelial cells and fibroblasts, were found to be significantly lower in tumor tissues (Fig. 1V–1W). Differential analysis with paired tumor and normal tissues showed similar patterns (Table 2).

Inhibitory receptors were co-expressed on certain subtypes of T cells

Inhibitory receptors, including PD1, CTLA4, LAG3, and TIM3, expressed on T cells often led to T-cell exhaustion allowing tumors to evade the immune response (Huang et al., 2017; Nirschl & Drake, 2013). The use of specific antibodies to inhibit CTLA4 or PD1 and overcome immune suppression and tumor regression is promising (Brahmer et al., 2012; Callahan, Wolchok & Allison, 2010).

We investigated the correlations between these inhibitory receptors and CD4+/CD8+ T cells. Heatmaps suggested that the expression patterns of the inhibitory receptors correlated with specific subsets of T cells with distinctions between tumor tissues and normal tissues (Fig. 2A and 2B). Correlation analyses also demonstrated that CD8+ T-cells, CD8+ Tcm, CD8+ naive T-cells, CD4+ memory T cells, and CD4+ naïve T cells were all positively correlated with expressions of these inhibitory receptors in tumor tissues, especially with PD1 and CTLA4 (P < 0.05) (Fig. 2C and 2D). CD8+ Tem, CD4+ Tcm, and CD4+ T-cells were not strongly correlated with the expressions of the inhibitory receptors (Fig. 2C). Only a few T cells were significantly correlated with these inhibitory receptors in normal tissues; TIM3 expression was negatively correlated with CD4+ Tcm (Fig. 2C and 2F). We observed a significant correlation among the expression of inhibitory receptors (Fig. 2E and 2G).

Cancer/testis genes TDRD6 and TTK show promise as breast cancer targets

Cancer/Testis (CT) genes are a cluster of tumor-associated proteins normally expressed in germ cells and different cancers. However, they are not typically seen in normal somatic cells (Scanlan et al., 2002). The limited expression of CT genes makes them ideal cancer and immunotherapy biomarkers.

We studied the antitumor immunity response to antigens generated by CT genes by examining 276 known CT genes (obtained from the CTDatabase) for their association with immune components. The significant associations between immune cells and CT genes (P < 0.001) are shown in Fig. 3A. Most of the adaptive immune cells were significantly correlated with CT genes. T cells, such as CD8+ T cells and aDC, which belonged to adaptive and innate immune responses, were positively correlated with most of the CT genes (Fig. 3A–3C). Moreover, two CT genes, TDRD6 and TTK, were positively correlated with a number of immune cells, especially the CD4+/CD8+ T cells (Fig. 3D and 3E), implying strong host immune reactions to these two cancer antigens.

Cellular heterogeneity correlated with clinical pathology of breast cancer

Cellular heterogeneity is an important part of the tumor microenvironment and is necessary for the growth and development of a tumor. We studied whether certain cell types were significantly correlated with certain clinical parameters, including age, sex, T stage, N stage, M stage, and TNM stage. A number of cell types were significantly correlated with clinical parameters, especially T stage and M stage (Fig. 4A). Astrocytes, mesangial cells, and mast cells were negatively correlated with the T stage and plasma cells were positively correlated with the T stage (Fig. 4B–4E). CD4+ Tcm, CD4+ Tem, microvascular (mv) endothelial cells, NKT, and MSC were all significantly higher at the M stage in patients with distal metastasis (Fig. 4F–4J). CLP was significantly higher at the N stage in patients with lymph node metastasis (Fig. 4K). Th1 cells and MSC were both positively correlated with TNM stage (Fig. 4L–4M). The 12 male breast cancer patients studied tended to have higher proportion of CLP and NKT compared with the female breast cancer patients (Fig. 4N–4O).

Prognostic model with survival associated cell types

Emerging evidence suggests that the number of tumor-infiltrating lymphocytes (TILs) of primary tumors consistently predicts favorable outcomes for a number of tumor types, including breast cancer. Therefore, survival analyses were performed to find survival-associated cell types within the tumor microenvironment (Fig. 5A). Immune cells were more strongly associated with overall survival, especially CD4+ and CD8+ T cells (Fig. 5A). Most T cells, including CD8+ T cells, CD8+ Tcm, CD8+ naïve T cells, and CD4 Tem, were favorable prognostic factors. However, high CD4+ naïve T cells were associated with worse overall survival (Fig. 5A–5M). NKT, class switched memory B cells, NK cells, cDC, and pDC were also significantly associated with overall survival (Fig. 5A–5M). A number of stromal cells, including mesangial cells and keratinocytes, were favorable prognostic factors but myocytes were adverse prognostic factors (Fig. 5A–5M). Multivariate COX regression revealed that CD8+ T cells, mesangial cells, keratinocytes, NKT, and class switched memory B cells were independent prognostic factors. We built a prognosis predictor model with five independent prognostic factors. Our model more reliably determined the survival of breast cancer patients with the highest AUC of ROC of 0.708, versus when the factors were analyzed separately (Fig. 5N, 5O).

Subtypes of breast cancer had diverse phenotypes of cellular heterogeneity

Emerging evidence suggests that the breast cancer transcriptome has a wide range of intratumoral heterogeneity, as well as genomic heterogeneity based on ER, PR, and Her2, which are determined by the tumor cells and immune cells in the surrounding microenvironment (Chung et al., 2017). We explored cellular heterogeneity among different subtypes of breast cancer by Her2, ER, and PR status. 1,092 breast cancer patients were classified into five groups according the clinicopathological parameters provided by TCGA, including 30 Her2+_HR- patients, 59 Her2+_HR+ patients, 426 Luminal type (Her2-_HR+) patients, 97 triple negative (Tri-negative) patients, and 480 unknown patients (Table 1). The relative proportion of different cells varied widely among these five subtypes (Fig. 6). Tri-negative patients had the highest fraction of immune cells while luminal type patients had the smallest fraction of immune cells, especially CD4+ and CD8+ T cells (Fig. 6A, Fig. S1). Cluster analyses from ImmuneScore, StromaScore, and heatmap suggested the absence of certain cell types used to distinguish these five subtypes (Fig. 6B). Furthermore, t-SNE cluster analysis suggested a large difference in tumor-infiltrating cells among these five subtypes (Fig. 6C and 6D). B cells, T cells, macrophages, Th cells and stromal cells, including keratinocytes, were significantly differentially enriched in these subtypes (Fig. 6E–6P and Fig. S2). Tri-negative breast cancer tissues had the highest fractions of plasma cells, pro B cells, macrophages M1, Th1, and Th2 cells, but M2 cells had the lowest fraction of macrophages (Fig. 6E–6L). Keratinocytes, sebocytes, and pericytes were found frequently in Tri-negative breast cancer whereas MSC cells were found in low amounts (Fig. 6M–6P). Survival analysis revealed interesting differences between the five subtypes (Fig. 7). Each subtype of breast cancer had a unique pattern of survival-associated tumor-infiltrating cells. Different cell types may have different functions in the prognosis of breast cancer patients. Keratinocytes had a favorable effect on the prognostic factors while neurons were associated with adverse prognosis factors in luminal type patients (Fig. 7A). However, in Tri-negative patients, keratinocytes predicted a worse overall survival and neurons predicted a better overall survival (Fig. 7D). Taken together, the diversity of cellular heterogeneity among the different subtypes of breast cancers suggested that tumor-infiltrating cells within the tumor microenvironment were essential in shaping the intratumor heterogeneity of breast cancer.