Tumor draining lymph nodes connected to cold triple-negative breast cancers are characterized by Th2-associated microenvironment

Abstract

Tumor draining lymph nodes (TDLN) represent a key component of the tumor-immunity cycle. There are few studies describing how TDLNs impact lymphocyte infiltration into tumors. Here we directly compare tumor-free TDLNs draining "cold" and "hot" human triple negative breast cancers (TDLN_Cold and TDLN_Hot). Using machine-learning-based self-correlation analysis of immune gene expression, we find unbalanced intranodal regulations within TDLN_Cold. Two gene pairs (TBX21/GATA3-CXCR1) with opposite correlations suggest preferential priming of T helper 2 (Th2) cells by mature dendritic cells (DC) within TDLN_Cold. This is validated by multiplex immunofluorescent staining, identifying more mature-DC-Th2 spatial clusters within TDLN_Cold versus TDLN_Hot. Associated with this Th2 priming preference, more IL4 producing mast cells (MC) are found within sinus regions of TDLN_Cold. Downstream, Th2-associated fibrotic TME is found in paired cold tumors with increased Th2/T-helper-1-cell (Th1) ratio, upregulated fibrosis growth factors, and stromal enrichment of cancer associated fibroblasts. These findings are further confirmed in a validation cohort and public genomic data. Our results reveal a potential role of IL4⁺ MCs within TDLNs, associated with Th2 polarization and reduced immune infiltration into tumors.

Fig. 1. Overview of the study.

We initially collected 15 pairs of tumor-free (tf) tumor draining lymph nodes (TDLN) and primary tumors (PT) of triple-negative breast cancer (TNBC) patients as our discovery set and stratified them into “cold” (n = 7) and “hot” (n = 8) cohorts based on percentages of stromal tumor infiltrating lymphocytes (TIL). NanoString® nCounter® Pan-Cancer immune profiling panel was used to quantify the gene expression of 730 immune-related genes for both TDLN and PT samples, and NanoString® nCounter® Pan-Cancer cancer pathway panel was used to quantify the gene expression of 730 cancer-related genes for PT samples only. To further analyze TDLN gene expressions, intranodal self-correlation analysis was conducted on TDLN_Cold and TDLN_Hot. To validate the findings from self-correlation analysis, we investigated the spatial distributions of mature dendritic cells (DC) and T helper (Th) cells using multiplex immunofluorescence staining. CD117 was also stained on TDLN samples to quantify mast cells. In addition, we collected another 20 tumor-free TDLNs (with similar stromal TIL percentages from paired primary tumors as in our discovery cohort) as an independent validation set. We co-stained IL4, CD117, and CD3 on 22 TDLN slides from both discovery and validation cohorts and confirmed the existence and location of IL4⁺ mast cells within TDLNs. To further investigate the impact of Th2 polarization in TDLNs on paired tumors, we implemented differential expression analysis, correlation analysis, and gene deconvolution on the paired PT samples in the discovery set as well as the large-scale public breast cancer database METABRIC. Using another public dataset, we also compared the gene expression profiles and immune cell compositions within micro-dissected stroma regions between cold and hot tumors to support a more fibrotic stroma in cold tumors.

Fig. 2. Intranodal self-correlation analysis reveals dysfunctional transcriptome-based immune regulation in tumor draining lymph nodes (TDLN_Cold) with cold tumors.

A Hierarchical clustering of gene-to-gene correlation matrices for TDLN_Cold (upper panel) and TDLN_Hot (lower panel). The cluster order of gene pairs in each cohort was mapped to the other cohort to show the differences of the correlation coefficients. B Higher Silhouette scores (p = 6.0 × 10⁻⁴) based on the k-means clustering (k = 2–12) of correlation coefficients in TDLN_Cold (n = 10) compared to TDLN_Hot (n = 10). 10 different cluster numbers (from 2 to 12) for k-means clustering algorithm were tested as technical replicates to obtain the statistical significancy. Each dot represents a cluster number. Paired two-sided t test was used here. C Opposite TBX21-CXCR1 correlations (p < 10⁻⁴) between TDLN_Cold (n = 7) and TDLN_Hot (n = 8). D Opposite GATA3-CXCR1 correlations (p = 1.7 × 10⁻³) between TDLN_Cold (n = 7) and TDLN_Hot (n = 8). 95% confidence bands of the best fit line were shown as shaded regions in (C) and (D). Two-sided Fisher’s z test was used for (C) and (D). Each dot represents a TDLN sample in (C) and (D). Source data are provided as a Source Data file.

Fig. 3. Differential spatial distributions between tumor draining lymph nodes (TDLN) with cold and hot tumors.

A, B Representative regions from mIF stained images showing the spatial distributions of mature dendritic cell (DC), T helper 1 (Th1), and T helper 2 (Th2) cells in TDLN_Cold (A) and TDLN_Hot (B). Whole slide images for these two samples are shown in Fig S4A, B. C Higher ratios of Th2 to Th1 cell numbers (p = 6.0 × 10⁻⁴) in TDLN_Cold (n = 7) compared to TDLN_Hot (n = 7). D Higher Th2 percentages to total CD4⁺ T cells (p = 6.0 × 10⁻⁴) in TDLN_Cold (n = 7) compared to TDLN_Hot (n = 7). E Similar cluster percentages of mature-DC-Th1 cluster (relative to total cluster number, p = 0.53) between TDLN_Cold (n = 7) and TDLN_Hot (n = 7). F Higher cluster percentages of mature-DC-Th2 cluster (relative to total cluster number, p = 0.05) in TDLN_Cold (n = 7) compared to TDLN_Hot (n = 6) for mature-DC-Th1 clusters. G The averages of cluster percentages across Th1/2 cell types and TDLN_Cold/Hot. Data in (C–F) were presented as mean values ± the standard errors of the mean. All the outliers were identified by ROUT method with Q = 1%. The statistical significances based on the two-sided Mann–Whitney test was shown in (C–F). Calculated p values are displayed as “ns”, p > 0.05; *, p ≤ 0.05; **, p ≤ 0.01; ***, p ≤ 0.001. Source data including the outliers are provided as a Source Data file.

Fig. 4. IL4 from mast cells (MC) skews naïve CD4⁺ T cell polarization to T helper 2 (Th2) cells in tumor draining lymph nodes (TDLN) with cold tumors.

A, B Upregulation of IL4 (A, p = 0.02) and TPSAB1 (B, p = 2.3 × 10⁻³) in TDLN_Cold (n = 7) compared to TDLN_Hot (n = 7). C–F Representative regions (100 μm scale bar) from H-DAB stained CD117 (dark color) in TDLN_Cold (C, discovery cohort; D, validation cohort) and TDLN_Hot (E, discovery cohort; F, validation cohort). G More MCs (p = 2.3 × 10⁻³) in TDLN_Cold (n = 7) compared to TDLN_Hot (n = 7) represented by CD117⁺ cell densities in discovery cohort. H Receiver operating characteristic (ROC) curve of CD117⁺ cell densities to predict the cold and hot primary tumor in the validation cohort (n_Cold = 10, n_Hot = 10). I, J Representative regions of interests (100 μm scale bars for larger views) with IL4 and CD117 co-staining and TDLN subregion annotations in TDLN_Cold (I) and TDLN_Hot (J). K Higher overall IL4⁺ MC densities (p = 0.01) in TDLN_Cold (n = 11) compared to TDLN_Hot (n = 11). L Higher overall IL4⁺ MC percentages relative to total MC numbers (p = 6.4 × 10⁻³) in TDLN_Cold (n = 11) compared to TDLN_Hot (n = 10). M Comparison of the IL4⁺ MC densities crossing B cell zones (n = 26), T cell zones (n = 29), and sinuses regions (n = 29) regardless of cold/hot cohort. The Kruskal–Wallis test with Dunn’s multiple comparison was used here. The multiplicity adjusted P-value was reported here for each comparison (p_{B-vs-T-cell-zone} = 8.0 × 10⁻⁴, p_{B-cell-zone-vs-Sinus} < 1.0 × 10⁻⁴, p_T-cell^-_{zone-vs-Sinus} = 0.97). N Comparison of the IL4⁺ MC densities between TDLN_Cold and TDLN_Hot within B cell zones (n_Cold = 12, n_Hot = 14, p_adj = 1.0), T cell zones (n_Cold = 14, n_Hot = 15, p_adj = 4.2 × 10⁻³), and sinuses (n_Cold = 14, n_Hot = 15, p_adj = 0.03) regions. Multiple Mann–Whitney tests with Holm–Sidak p-value adjustment were used here. O Comparison of the IL4⁺ MC percentages between TDLN_Cold and TDLN_Hot within T cell zones (n_Cold = 15, n_Hot = 11, p_adj = 0.41) and sinuses (n_Cold = 14, n_Hot = 14, p_adj = 9.3 × 10⁻³). Multiple Mann–Whitney tests with Holm–Sidak p-value adjustment were used here. To improve the reliability of the IL4⁺ MC percentage, all the data points whose total MC numbers were smaller than 3, were excluded in the comparisons of Fig. 4O (excluded sample number: n_Cold = 1, n_Hot = 3). Data in (A, B, G, K–O) were presented as mean values ± the standard errors of the mean. The statistical significances based on the two-sided Mann–Whitney test(s) were shown in this figure except the subfigures with specific explanations. All the outliers were identified with Grubbs’ method (alpha = 0.05). Calculated p values are displayed as “ns”, p > 0.05; *, p ≤ 0.05; **, p ≤ 0.01; ***, p ≤ 0.001; ****, p ≤ 0.0001. Source data including outliers are provided as a Source Data file.

Fig. 5. Impact of preferred T helper 2 (Th2) cell polarization to microenvironment of primary tumors (PT).

A, B The ratios of Th2/Th1 cell type scores were negatively correlated with CD45 (A) and CD3D (B) expression from our paired primary tumor cohort (n = 15). Higher ratios of Th2/Th1 cell type scores in PT_Cold (n = 7) compared to PT_Hot (n = 8). C Th2 and tissue-repairing macrophages (TsRpM) cell type scores were significantly correlated from our paired PT cohort (n = 15). D TsRpM and cancer associated fibroblast (CAF) cell type scores were significantly correlated from our paired tumor cohort (n = 15). E–G Higher ratios of Th2/Th1 cell type scores (E, p = 6.0 × 10⁻⁴), higher relative IL4 expression from T cells (CD3D) (F, p = 1.2 × 10⁻³), and higher relative PDGFC expression from macrophages (CD68) (G, p = 3.0 × 10⁻⁴) in cold paired PTs (n = 7) compared to hot ones (n = 8). H, I The ratios of Th2/Th1 cell type scores were negatively correlated with CD45 (H) and CD3D (I) expression in the PTs from METABRIC-TNBC cohort (n = 299). Higher ratios of Th2/Th1 cell type scores. J Th2 and TsRpM cell type scores were significantly correlated in the PTs from METABRIC-TNBC cohort (n = 299). K TsRpM and CAF cell type scores were significantly correlated in the PTs from METABRIC-TNBC cohort (n = 299). L–N Higher ratios of Th2/Th1 cell type scores (L, p < 1.0 × 10⁻⁴), higher relative IL4 expression from T cells (CD3D) (M, p < 1.0 × 10⁻⁴), and higher relative PDGFC expression from macrophages (CD68) (N, p < 1.0 × 10⁻⁴) in cold PTs (n = 105) compared to hot ones (n = 105) in METABRIC-TNBC cohort. O T cell scores in cold tumors (n_Epithelium = 16, n_Stroma = 16) were significantly higher than the ones in hot tumors (n_Epithelium = 22, n_Stroma = 22) in both epithelium (p_adj = 2.8 × 10⁻³) and stroma (p_adj = 1.7 × 10⁻³) regions. Multiple Mann–Whitney tests with Holm–Sidak p-value adjustment were used here. P Average T cell scores crossing cold/hot tumors and epithelium/stroma regions. Q CAF scores in cold tumors (n_Epithelium = 16, n_Stroma = 16) were significantly higher than the ones in hot tumors (n_Epithelium = 22, n_Stroma = 22) in stroma regions (p_adj = 2.2 × 10⁻³) but not epithelium regions (p_adj = 0.08). Multiple Mann–Whitney tests with Holm–Sidak p-value adjustment were used here. R Average CAF scores crossing cold/hot tumors and epithelium/stroma regions. Data in (E–G, L–O, Q) were presented as mean values ± the standard errors of the mean. All the outliers were identified by ROUT method with Q = 1%. 95% confidence bands of the best fit line were shown as shaded regions in (A–D) and (H–K). The null hypothesis testing of no correlation (p values shown in the subfigures) was implemented in (A–D) and (H–K). Calculated p values are displayed as “ns”, p > 0.05; *, p ≤ 0.05; **, p ≤ 0.01; ***, p ≤ 0.001; ****, p ≤ 0.0001. Source data including outliers are provided as a Source Data file.