The expression profiles and clinical information of BCa patients were downloaded from TCGA via GDC Legacy Archive (113 normal cases and 1098 BCa cases with FPKM data). The expression profiles and clinical information of the METABRIC BCa cohort were downloaded from cBioPortal (https://www.cbioportal.org/). The gene set of cytokine activity was obtained from GSEA (https://www.gsea-msigdb.org/gsea/index.jsp). The GEPIA (http://gepia.cancer-pku.cn/) website was used to detect gene expression levels in TCGA.

A gene co-expressing analysis was performed by screening mRNA profiles in the METABRIC dataset using the WGCNA R package (10). BCa cases in the METABRIC cohort were classified into six groups based on clinical information (claudin-low: 199 patients; basal: 199 patients; HER2+: 220 patients; Luminal A: 679 patients; Luminal B: 461 patients). The top 200 genes in the most significant module were extracted to perform protein–protein interaction network (PPI) construction via the STRING online tool (Search Tool for the Retrieval of Interacting Genes, https://string-db.org/, version: 11.0). The PPI network was input to Cytoscape software (11). The most significant module of the PPI network was detected by the MCODE which was constructed for clustering a network based on topology to determine intensively connected Regions (12), and then hub genes in this module were found by using CytoHubba (13).

Differentially expressed genes (DEGs) were analyzed in four-pair-wise comparisons by the R package of linear models for microarray data (“limma”) (14). Pathway enrichment analysis was conducted on the DEGs with |log2-fold change (FC)| ≥ 1 using Kyoto Encyclopedia of Genes and Genome (KEGG) via the “clusterProfiler” R package. DEGs with log2-fold change (FC) ≥ 1.5 intersected with the gene sets associated with cytokines from Molecular Signatures Database (MSigDB). The TCGA-BRCA dataset was divided into two groups by the median of the expression levels of CLDN3/4/7, OCLN and CDH1 (172 cases in claudin-low TP and non-claudin-low TP, respectively). Claudin-low BCa were identified in the Oslo2 and TCGA-BRCA databases by the nine-cell line claudin-low predictor using the R package of “Genefu” (15), and then expression levels of interested cytokines and claudin proteins were also confirmed in CL-BCa of the Oslo2 dataset and TCGA-BRCA database (Oslo2 dataset: 13 CL-BCa cases and 352 non-CL-BCa cases; TCGA-BRCA: 33 CL-BCa cases and 1058 non-CL-BCa cases).

Infiltration of immune cells in 1,758 BCa patients of the METABRIC dataset were investigated using the R package of “immunedeconv” (16) and web server of TIMER 2.0 (http://timer.cistrome.org/) which is an online tool for assessing the specific gene(s) associated with tumor immune infiltrating cells (17). Relative proportions of tumor-infiltrating immune cells were extracted and classified according to the different algorithms. Combined with expression levels in the TCGA dataset, all possible target genes were selected for survival analysis. Survival analysis was conducted by the “survminer” and “survival” R package. Kaplan–Meier survival curves were conducted for each group, respectively, and compared statistically using the log-rank test for testing the null hypothesis of no difference in survival between two or more independent groups. Multivariable Cox regression analyses were conducted by the “forestplot” R package. Pearson correlations were calculated with the “cor.test” package in R, and correlation plots were prepared with the “ggpubr” and “corrplot” packages in R. The Wilcoxon test and Kruskal–Wallis H test were used to determine significance for non-parametric data; p value <0.05 was considered statistically significant. All analyses were conducted using R version 4.0.2.

The relative expression of DEGs was examined by quantitative polymerase chain qPCR (Bio-Rad CFX Connect™ Optics Module) and was normalized to the expression of β-actin using the 2^-ΔΔCt method according to a previous study (18). Total RNA was extracted with TRIzol reagent (TaKaRa, Japan, Cat#9109), and cDNA was reverse transcribed subsequently with the PrimeScript RT Reagent Kit (TaKaRa, Japan, Cat#RR037A). qPCR was then performed using a SYBR Premix Ex Taq II (TaKaRa, Japan, Cat#RR820A) according to the manufacturer’s instructions. The sequences of qPCR primers used were presented as follows: CCL5, 5′-CAGCAGTCGTCCACAGGTCAAG-3′ (forward) and 5′-TTTCTTCTCTGGGTTGGCACACAC-3′ (reverse); CCL19, 5′-AGCCTGCTGGTTCTCTGGACTTC-3′ (forward) and 5′-AGGGATGGGTTTCTGGGTCACAG-3′ (reverse); CXCL9, 5′-TCTTGCTGGTTCTGATTGGAGTGC-3′ (forward) and 5′-GTCCCTTGGTTGGTGCTGATGC-3′ (reverse); CXCL10, 5′-AACTGTACGCTGTACCTGCATCAG-3′ (forward) and 5′-ACGTGGACAAAATTGGCTTGCAG-3′ (reverse).

The expression of DEGs (CCL5 and CCL19) in BCa patients was examined by immunohistochemistry (IHC). Antibodies used for IHC were CCL5 (Cat#abs136939, Absin, Shanghai, China), CCL19 (Cat#abs149041, Absin), iNOS (Cat#GB11119, Servicebio, Wuhan, China), and CD8 (Cat#GB13068, Servicebio). The IHC method was followed as in our previous publication (19). IHC staining was examined with microscopy. IHC optical density scores were calculated to quantify protein expression in IHC images (20). CL-BCa samples were identified by the low expression of CLDN3/CLDN4/CLDN7/CDH1 in patients’ tissue biopsies. All patients were diagnosed and treated at The First Affiliated Hospital of Chongqing Medical University, and the informed consent was obtained from all patients.

To identify the modules of gene expression in claudin-low BCa using the WGCNA package in R, patients in METABRIC were divided into two groups, claudin-low (CL) group (199 cases) and non-claudin-low (NCL) group (1,559 cases). The top 5,000 genes with similar patterns of expression were categorized into modules via hierarchical average linkage clustering. Totally, 15 modules were identified between the CL group and NCL group and the relevance of each module was calculated ( Figures 1A, B ). We found that the brown module was the most relevant module to the trait of the claudin-low subtype (Cor = 0.83, p = 2.2e-130) ( Figure 1B ). Hub genes of the brown module also tend to be highly correlated with weight, which suggested that these hub genes might be the featured genes of claudin-low subtype ( Figure 1C ). Genes identified in the brown module belonged to gene sets associated with cytokine interaction ( Figure 1D ). The top 200 genes ranked by co-expressing correlation in the brown module were chosen to perform PPI network (protein–protein interaction network) construction using the STRING online database. The top-scoring module in the PPI network (score: 22.074) was analyzed using MCODE in Cytoscape software. Specially, we found that the hub genes in this module were CXCL9, CXCL10, CCL5, TLR8, CD8A and EOMES, which were relevant to the immune pathway and were detected by CytoHubba ( Figure 1E ).

To uncover the characteristics of claudin-low BCa, differentially expressed genes (DEGs) were used to portray the difference at the transcription level between CL-BCa and other subtypes of BCa. We found that the overall gene expression patterns at the mRNA level were significantly distinct between claudin-low and other subtypes ( Figure 1F ). Compared to the basal subtype, 67 DEGs were downregulated and 180 DEGs were upregulated in the claudin-low subtype. Relative to the HER2+ subtype, 333 DEGs were identified in claudin-low subtype and 204 DEGs were upregulated. Additionally, 515 and 584 differently expressed genes were identified in the claudin-low subtype compared with luminal A or luminal B, respectively. Since the genes identified in the brown module were enriched in the pathway associated with cytokine interaction ( Figure 1D ), we further investigated the cytokines which were upregulated in the claudin-low group by intersecting upregulated DEGs with the gene sets of cytokine-associated pathways. Compared with the basal or HER2+ subtype and using the 1.5 fold change as the threshold, three cytokines upregulated in claudin-low were identified (CXCL12, CCL15, CCL19) ( Figure 1G ). Compared with the luminal A subtype, nine cytokines (TNFSF13B, CCL8, LTB, CCL2, CCL5, CXCL9, CCL19, CXCL10, CXCL13) were upregulated in claudin-low; similarly, we found 10 cytokines in claudin-low compared to the luminal B subtype (CCL19, CX3CL1, LTB, CCL2, CCL5, CCL21, CXCL12, CCL15, CXCL13, CXCL9) ( Figure 1G ). These findings indicated that the brown module in WGCNA analysis was the most relevant to CL-BCa subtype and cytokines genes were remarkably DEGs in CL-BCa.

Since the result of the WGCNA analysis hinted that cytokines were distinctive feature genes of claudin-low BCa subtype, we intended to explore the signal pathways that these differentially expressed cytokines were involved. In order to portray the signal pathways of these DEGs in claudin-low BCa, enrichment analysis on the DEGs was performed to clarify DEG-related pathways in the claudin-low subtype of the METABIRC dataset. Generally, pathways of enriched DEGs in CL-BCa were mainly associated with cytokine signals and cell adhesion ( Figures 2 , 3 and Tables 1 , 2 ). Compared to the basal subtype, upregulated genes were mainly involved in pathways of hematopoietic lineage and chemokine signaling ( Figure 2A ). Especially, the CCL15, CCL19, CXCL12 and CXCL14 in chemokine signaling pathways were the most upregulated chemokine genes in the claudin-low subtype ( Figure 2B ). Relative to the HER2+ subtype, upregulated DEGs mainly focused on signaling pathways of cytokine–cytokine receptor interaction and chemokine ( Figure 2C ). Among these DEGs, CXCL12, CCL2, CCL15 and CCL19 were top upregulated chemokine genes, while CX3CR1, CSF1R, IL33, IL7R and LTB were the most enhanced genes on cytokine–cytokine receptor interaction ( Figure 2D ). Next, enrichment analysis on DEGs between the claudin-low and luminal A or luminal B subtype was also performed by METABRIC. Specifically, compared to the luminal A subtype, the DEGs in the claudin-low subtype were mostly enriched in the chemokine signaling pathway, phagosome, and cell adhesion ( Figure 3A ). As upregulated genes, CXCL10, CXCL9, CCL19, CCL5 TAP1, SELL and VCAM1 were enhanced in the claudin-low subtype ( Figure 3B ). Similarly, compared to the luminal B subtype, genes related to virus infection and viral protein with the cytokine and cytokine receptor were mostly enhanced ( Figure 3C ). For instance, CCL19, CCL5, CDK6, CXCL10, IL6, IL33 and IL10RA were upregulated in the claudin-low subtype ( Figure 3D ). Meanwhile, the expression of downregulated genes in claudin-low tissues versus the other four subtypes was also identified ( Supplementary Figures 1A–F ). Taken together, based on the result of the WGCNA analysis and featured DEGs and their associated pathways in CL-BCa, cytokine/chemokine signaling might play a pivotal role in claudin-low BCa.

The preliminary result in our study has clarified that mostly upregulated DEGs in the claudin-low subtype were related to the cytokine–cytokine receptor interaction and chemokine signaling pathway compared to other BCa subtypes. Considering various interactions of cytokine/chemokine with immune cells, we sought to clarify the composition of tumor-infiltrating immune cells in claudin-low and non-claudin-low BCa tissue using the R package of “immunedeconv”. Generally, the fraction of tumor-infiltrating immune cells in claudin-low BCa tissues was significantly enhanced relative to non-claudin-low BCa. We found that T cells in claudin-low BCa strikingly increased compared with other BCa subtypes, especially T cell CD4+ memory, T cell CD8+ central memory and T cell CD4+ Th2 ( Figure 4A ). Specifically, the fraction of these T cells displayed distinctive differences between claudin-low BCa tissue and basal/HER2+ tissue, which indicated that increased T cell infiltration was a trait for claudin-low BCa ( Figure 4A ). In addition, claudin-low BCa tissue retained a higher fraction of immune cells than luminal A/luminal B tissues, including T cell CD4+ memory, T cell CD4+ naïve, T cell central/effector memory and T cell CD4+ Th2 ( Figure 4A ). Further, we analyzed the composition of macrophage in claudin-low and non-claudin-low tissues assessed by XCELL ( Figure 4B ). Compared with the basal or HER2+ subtype, fractions of both M1 and M2 macrophage were higher in tissues of the claudin-low subtype, which suggested that macrophages were probably more engaged in claudin-low BCa development ( Figure 4B ). Similarly, compared with the luminal A or luminal B subtype, M1 and M2 macrophage also showed a higher fraction in the claudin-low subtype ( Figure 4B ). Based on these observations, CD8+ T cell, CD4+ T cell, and macrophages were the distinctively enriched immune cells in claudin-low BCa.

In order to verify the enriched expression of featured DEGs identified by KEGG, we examined the expression levels of featured DEGs (cytokines and claudins) in METABRIC and other BCa databases. Firstly, the expression levels of featured cytokines were examined in the METABRIC dataset. Seven differentially expressed cytokines (CCL15, CCL2, CCL19, CCL5, CXCL12, CXCL9 and CXCL10) had a higher expression level in the claudin-low subtype than non-claudin low subtypes ( Figures 4C, D ). Then, to confirm that these featured cytokines and claudins were also distinctively expressed in claudin-low BCa rather than in normal breast, we sought to screen the expression levels of these DEGs in the TCGA-BRCA dataset to detect the differences in their expression between claudin-low BCa tissue and normal breast tissue. Due to the lack of the claudin-low cohort in the TCGA-BRCA dataset, the claudin-low group was self-defined by dividing the TCGA-BRCA cohort according to the expression levels of CLDN3 and CLDN4 and CLDN7 and Occludin (OCLN) and E-cadherin (CDH1) (172 cases in the Claudin-low TP group and 113 cases in the normal group). We found that compared to the normal breast group, 6,270 downregulated genes and 4,770 upregulated genes were identified in the claudin-low BCa group ( Figure 4E ). Among these differentially expressed genes, CCL5, CCL19, CXCL9 and CXCL10 were upregulated in the claudin-low BCa group which were also identified as overexpressed genes in the METABRIC dataset; six claudin proteins (CLDN1, CLDN2, CLDN12, CLDN19, CLDN11 and CLDN8) were downregulated claudin proteins in the claudin-low BCa group compared to the normal group in addition to CLDN3&4&7 ( Figure 4F ). CLDN8, CLDN11 and CLDN19 were top three of the significantly downregulated claudins in our self-defined claudin-low group compared to the normal group ( Figure 4F ).

To further confirm the expression pattern of featured cytokines and claudins, we identified claudin-low patients in the Oslo2 database and TCGA dataset by the nine-cell line claudin-low predictor with the use of the R package “Genefu” (Oslo2 dataset: 13 CL-BCa cases and 352 non-CL-BCa cases; TCGA-BRCA: 33 CL-BCa cases and 1058 non-CL-BCa cases). Although the CL-BCa cases in the Oslo2 and TCGA databases were limited, CLDN11, CLDN19 and CLDN6 were downregulated claudins in CL-BCa compared to non-CL-BCa and CCL5, CXCL9 and CXCL10 were upregulated cytokines in CL-BCa compared to non-CL-BCa ( Supplementary Figures 2A, B ). Additionally, we investigated correlations between the expression of CLDN3/4/7 and cytokines (CCL5, CCL19 and CXCL9) expression in the GSE25066 dataset ( Supplementary Figure 4 ). The negative correlations were shown between expressions of CLDN3/CLDN4/CLDN7 and CCL5/CCL19/CXCL9, which suggested that these upregulated cytokines (CCL5, CCL19 and CXCL9) could be used as potential biomarkers of CL-BCa patients with the feature of low-expressed CLDN3/CLDN4/CLDN7 ( Supplementary Figure 4 ).

Finally, the potential functions of these cytokines and claudins in TCGA-BRCA were investigated by CANERSEA. These four cytokines were mostly correlated with inflammation, differentiation, and angiogenesis, which had an effect of promoting tumorigenesis ( Figure 4G ). It was worthwhile to mention that CCL5, CXCL9 and CXCL10 were the hub genes in the significant PPI module calculated by the DMNC algorithm, which indicated that these three cytokines probably played a crucial role in immune network ( Figure 1E ). The three featured claudin proteins were also negatively correlated with DNA damage, apoptosis, invasion, and metastasis ( Figure 4H ). Taken together, in claudin-low BCa, four cytokines (CCL5, CCL19, CXCL9 and CXCL10) were upregulated and three claudin proteins (CLDN8, CLDN11 and CLDN19) were downregulated, which were consequently selected as the genes of interest for the following analysis.

To further investigate the relation of featured cytokine expression with immune cell recruitment, correlations of featured cytokines and immune cells were assessed by XCELL, CIBERSORT, and QUANTISEQ via TIMER2.0 in the BCa patient cohort of TCGA-BRCA. Firstly, we examined the correlation between each of the featured cytokines and tumor purity in TCGA-BRCA cohort which reflected the proportion of tumor cells in this sample cohort. The results of tumor purity showed that three cytokines (CCL5, CCL19 and CXCL9) were negatively correlated with tumor purity, which indicated that tumor cells were reduced but immune cells increased in tissues with high expression of these cytokines ( Figure 5A ).

Next, the results showed that there were positive correlations between CCL5 and M1 macrophage infiltration assessed by XCELL and CIBERSORT (Rho = 0.487 and 0.62, respectively) but no positive correlation between CCL5 and M2 macrophage by CIBERSORT (Rho = -0.366), which indicated that a high expression of CCL5 might facilitate infiltration of M1 macrophage into BCa tumor tissue ( Figures 5B, C ). Similar correlations were shown for CXCL9, which indicated that increased CXCL9 might enhance M1 macrophage infiltration into BCa tumor tissue ( Figures 5D, E ). Besides, CCL5, CCL19 and CXCL9 were positively correlated with CD8+ T cell infiltration assessed by QUANTISEQ (Rho = 0.818, 0.578, and 0.716, respectively), which suggested that a high expression of CCL5, CCL19, and CXCL9 might promote the tumor infiltration of CD8+ T cells in CL-BCa ( Figures 5F–H ). Finally, CXCL10 was positively correlated with CD4+ T cells, especially Th2 CD4+ T cell (Rho = 0.58) and memory CD4+ T cell (Rho = 0.544), which indicated that CXCL10 was associated with CD4+ T cell infiltration ( Figures 5I, J ).

To further confirm the connection between immune cells (CD8+ T cell and M1 macrophage) and cytokines in claudin-low BCa, we analyzed immunohistochemistry staining images from BCa. Among these featured cytokines, the expression of CCL5 was more relevant to the expression level of the M1 macrophage marker (iNOS) and CD8+ T cell marker (CD8). IHC results showed that a higher expression of iNOS and CD8 was observed in CCL5 highly expressed CL-BCa tissues, compared to non-CL-BCa tissues with a low CCL5 expression ( Figure 5K ). The results indicated that more M1 macrophages and CD8+ T cells existed in CL-BCa tissues with CCL5 high expression. In conclusion, we found that the protein expressions of CCL5, iNOS, and CD8 in CL-BCa tissues were higher than those in non-CL-BCa tissues and that more infiltration of M1 macrophages and CD8+ T cells was found in CCL5 highly expressed CL-BCa tissue. Taken together, the expression of these featured cytokines was associated with immune cell recruitment in CL-BCa tissues, which possibly altered the composition of immune cells in the tumor environment and thereby affected the prognosis of CL-BCa patients.

To clarify the relation between featured cytokines and claudin proteins with patients’ prognosis, Kaplan–Meier survival analysis was performed on the four upregulated cytokines (CCL5, CXCL9, CCL19 and CXCL10) and three downregulated claudin proteins (CLDN19, CLDN11 and CLDN8). Generally, high expression levels of CCL5, CCL19, CXCL9 and CXCL10 implied better prognosis, especially for the BCa patients with a negative ER status ( Figures 6A–H ). Firstly, expression levels of CCL5, CCL19, CXCL9 and CXCL10 were examined in TCGA-BRCA dataset for their relation with overall survival of BCa patients. Except CXCL10, high expression levels of CCL5, CCL19 and CXCL9 were significantly associated with better prognosis in breast cancer patients ( Figures 6A–D ). Additionally, multivariable Cox regression analyses were also performed in TCGA-BRCA dataset. A low expression of CCL5, CCL19 and CXCL9 was associated with poor outcome of BCa ( Supplementary Figure 5 ).

Meanwhile, the relation between expressions of these cytokines and overall survival of BCa was also checked in the METABRIC dataset. Considering the lower expression of ER and ER-related genes in claudin-low BCa tissues than luminal tissues, we further investigated the association of these cytokines with overall survival of ER-positive BCa patients or ER-negative BCa patients. The result of Figures 6E, F indicated that patients with high expressions of CCL5 or CCL19 had a significantly longer overall survival in the cohort of ER-negative BCa patients rather than ER-positive BCa patients ( Figures 6E, F ). Moreover, we refined the claudin-low cohort in the METABRIC dataset by the nine-cell line claudin-low predictor and then evaluated the relation between CCL5/CCL19 with overall survival of CL-BCa patients (8, 15). A higher expression level of CCL5 was related with enhanced survival probability of ER-negative BCa patients, despite of the limited number of CL-BCa patients (23 ER-negative cases and 66 ER-positive cases) ( Figure 6G ). Likewise, higher CCL19 was also associated with better survival probability of ER-negative BCa patients ( Figure 6H ).

In addition, to further confirm the overexpression of these featured cytokines in claudin-low BCa, we detected their expression levels by qPCR in CL-BCa cell lines MDA-MB-231 and HS578T (ER-negative) and non-CL-BCa cell line MCF7 (ER-positive). Generally, the qPCR result indicated that these cytokines (CCL5, CXCL9, CCL19 and CXCL10) had a higher expression level in CL-BCa cell lines compared with the non-CL BCa cell line ( Figure 6I ). Specifically, compared to MCF7, the expression levels of CCL5, CXCL9, CCL19 and CXCL10 were significantly enhanced in MDA-MB-231, while the expression levels of CCL19 and CXCL10 were significantly raised in HS578T ( Figure 6I ). Finally, to further confirm the overexpression of CCL5 and CCL19 in BCa patients, we also examined the expression of CCL5 and CCL19 by IHC in both ER-positive BCa tissues and ER-negative CL-BCa tissues. Compared to in ER-positive BCa tissue, a higher expression of CCL5 and CCL19 was identified in ER-negative CL-BCa tissue, which was consistent with their expression in cell lines ( Figures 6J, K ). Therefore, regarding their significantly overexpression in CL-BCa, these featured cytokines had the potential as the prognosis markers for CL-BCa.

A Kaplan–Meier survival analysis on featured claudin proteins was also performed in TCGA-BRCA and METABRIC datasets. In TCGA-BRCA dataset, unfortunately, we found that none of CLDN3, CLDN4 and CLDN7 showed a statistical relationship with survival of BCa patients, which suggested that CLDN 3/4/7 has no potential to be used as a predictive marker for overall survival ( Figures 7A–C ). However, a low expression of CLDN8, CLDN11 and CLDN19 was significantly related to shorter survival of BCa patients ( Figures 7D–F ), which was consistent with their lower expression levels in claudin-low BCa tissues ( Figure 4F ). Further, multivariable Cox regression analyses were also performed in TCGA-BRCA dataset. A low expression of CLDN8, CLDN11 and CLDN19 was related to poor outcome of BCa ( Supplementary Figure 5 ). In the METABRIC dataset, for the cohort of overall BCa patients, we surprisingly found that low-expressed CLDN11 was related to shorter overall survival in ER-positive BCa patients, but there was no significant difference in ER-negative patients ( Figures 7G, H ). By evaluating the association of CLDN8 and CLDN11 with overall survival in the refined claudin low cohort in METABRIC, we found that patients with low expression of CLDN8 or CLDN11 showed shorter overall survival time; however, low expression of CLDN11 was more significantly associated with worse overall survival in CL-BCa patients ( Figures 7I, J ). Combining the results of the survival analysis and expression levels of these claudins in databases, CLDN8 and CLDN11 may be the featured members of the claudin family which were associated with the overall survival of CL-BCa, rather than CLDN 3&4&7. However, more studies are needed to explore the specific function of these featured claudins in BCa.