Landscape of cancer-associated fibroblasts identifies the secreted biglycan as a protumor and immunosuppressive factor in triple-negative breast cancer.

Cancer-associated fibroblasts (CAFs) are essential for tumor microenvironment remodeling and correlate with tumor progression. However, interactions between CAFs and tumor cells and immune cells in triple-negative breast cancer (TNBC) are still poorly explored. Here, we investigate the role of CAFs in TNBC and potential novel mediators of their functions. The clustering of classic markers was applied to estimate the relative abundance of CAFs in TNBC cohorts. Primary fibroblasts were isolated from normal and tumor samples. The RNA and culture medium of fibroblasts were subjected to RNA sequencing and mass spectrometry to explore the upregulated signatures in CAFs. Microdissection and single-cell RNA sequencing datasets were used to examine the expression profiles. CAFs were associated with hallmark signalings and immune components in TNBC. Clustering based on CAF markers in the literature revealed different CAF infiltration groups in TNBC: low, medium and high. Most of the cancer hallmark signaling pathways were enriched in the high CAF infiltration group. Furthermore, RNA sequencing and mass spectrometry identified biglycan (BGN), a soluble secreted protein, as upregulated in CAFs compared to normal cancer-adjacent fibroblasts (NAFs). The expression of biglycan was negatively correlated with CD8 + T cells. Biglycan indicated poor prognostic outcomes and might be correlated with the immunosuppressive tumor microenvironment (TME). In conclusion, CAFs play an essential role in tumor progression and the TME. We identified an extracellular protein, biglycan, as a prognostic marker and potential therapeutic target in TNBC.

Introduction

Breast cancer is the most common malignant tumor in women worldwide.[1] Triple-negative breast cancer (TNBC) is a special molecular subtype that is characteristic of negative hormone receptors (estrogen and progesterone receptors) and human epidermal growth factor receptor 2 (HER2). Neither endocrine therapy nor routine targeted therapy is effective for TNBC. There are limited therapeutic selections for TNBC and some patients are still burdened by low efficacy of treatment response and high risk of recurrence or metastasis.[2] Hence, it is crucial to explore the mechanism underlying TNBC development and progression, which may lay a foundation for its diagnosis and treatment.

The crucial role of the tumor microenvironment (TME), which serves as the soil for seeds (cancer cells), has been proven in many studies.[3-5] Cells in the TME mainly include stromal cells and immune cells. Recently, increasing evidence has highlighted that appropriate stromal cells are crucial for the development of tumors.[6-8] Among them, cancer-associated fibroblasts (CAFs) represent the main fraction, and accumulating evidence has indicated their role in cancer proliferation, progression and invasion.[9,10] It is well acknowledged that CAFs interact with malignant cells and orchestrate the metastasis of breast cancer.[11-13] Although various clinical trials targeting CAFs have been performed in recent years, such as targeting surface markers, reducing CAF infiltration and normalizing CAF functions, most of them are still ongoing.[10] Previous studies have identified many CAF markers, but few of them have moved into clinical practice.[14] This may be due to the internal heterogeneity of CAFs. The CAFs seemed to originate from diverse cell types, such as fibrocytes, stellate cells, endothelial cells and mesenchymal stem cells.[14] It is well accepted that most activated fibroblasts are derived from fibroblasts of adjacent normal tissues and induced by oxidative stress or specific cytokines and chemokines from cancer cells.[15] Hence, distinct subclusters have been identified by previous studies. CAFs in human breast cancer can be divided into CAF-S1 to CAF-S4 according to specific signatures, such as CD29, FAP, FSP1, α-SMA, CAV1 and DPP4.[16] Furthermore, the CAF-S1 cluster contributes to a poor response to immunotherapy.[17] Immunotherapy has become a novel strategy for breast cancer treatment in recent years.[18] Therefore, focusing on the function and mechanism of fibroblasts in the tumor microenvironment may provide a strategy for breast cancer treatment, especially immunotherapy.

In this study, we explored the relative infiltration level of fibroblasts in triple-negative breast cancer and the correlation between CAFs and immune components in the TME. We further identified an upregulated secreted protein, biglycan, in CAFs compared with normal cancer-adjacent fibroblasts (NAFs) using RNA sequencing and mass spectrometry methods. Biglycan is encoded by the BGN gene and is mainly expressed in the stromal part of tumors. Biglycan is a protein that belongs to the small leucine-rich proteoglycan (SLRP) family. We found that the expression level of BGN is correlated with the extracellular matrix, lymphangiogenesis, epithelial-mesenchymal transition, angiogenesis and TGF-β signaling. Single-cell sequencing results show that BGN is mainly expressed in stromal fibroblasts. Moreover, BGN is highly expressed in CAFs of TNBC compared with other cell subpopulations. The role of BGN in the TME and its mechanism underlying how to affect the microenvironment remain unknown.

Establishing the relevance of the role of fibroblasts and biglycan in human triple-negative breast cancer using publicly available datasets and clinical samples raises the probability that targeting biglycan may yield clinical utility.

Methods

Datasets and tissue specimens

The Cancer Genome Atlas (TCGA) dataset was obtained using the TCGAbiolinks package in R software.[19,20] The transcripts per million (TPM) value was estimated at the transcript level. Patients who were diagnosed with breast cancer histologically and available for transcriptomic data were included. To further distinguish triple-negative breast cancer samples from the breast cancer cohort (BRCA), patients with immunohistochemically negative estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor 2 (HER2) were included. A total of 116 patients were enrolled in the TCGA-TNBC cohort. Overall survival (OS) was assessed using vital status and days from diagnosis to death or the last follow-up date. Only patients with active follow-up information were included in survival analysis.

The processed data of the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) dataset were downloaded from cBioPortal (http://www.cbioportal.org/) according to the website’s guidance.[21] No identification information of participants was involved during download and analysis. Patients diagnosed with breast cancer and available for active follow-up information were included. Triple-negative breast cancer samples were also enrolled in the METABRIC-TNBC cohort as described above.

For other TNBC datasets (GSE25066, GSE103091, GSE21653 and GSE88715), the expression matrices were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/).[22-25] The probes were mapped using the corresponding annotation platforms. The expression values were further normalized by the limma package if necessary.[26] The secreted genes were predicted as previously reported.[27] Single-cell sequencing (scRNA-seq) datasets were acquired and analyzed according to the guidance of previous studies. The cell types were annotated using the SingleR package if necessary.[28,29] GSE114727, GSE136206, GSE138536 were used to observe the expression of CAF markers in single-cell level.[30-32] The BGN expression was also analyzed in cohorts (Bassez A, et al.; GSE118389; Wu SZ, et al.).[33-35] These datasets are shown in Table S1. No identification information of participants was involved during download and analysis.

Tissues for immunohistochemistry and primary cell isolation were obtained at Sun Yat-sen University Cancer Center (SYSUCC) between January 2010 and June 2021. The clinical breast cancer specimens were collected with permission from the Institutional Review Board of the SYSUCC and informed consent was obtained from participants. The study was conducted according to the principles stated in the Declaration of Helsinki.

Tumor microenvironment estimation

The Estimate the Proportion of Immune and Cancer cells (EPIC), xCell, and Microenvironment Cell Populations-counter (MCP-counter) algorithms were applied to calculate the cancer-associated fibroblast scores in datasets.[36-38] To analyze the correlation among fibroblasts and immune cells, fractions of 22 immune cells were estimated using the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm.[39-41] The estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) was applied to calculate the overall stromal and immune scores in cancer.[42]

Molecular markers and CAF clustering

Representative immune-related genes (IRGs) and CAF signatures were obtained from previous studies.[3,43,44] The IFN-γ signature (IFNG1, STAT1, CXCL10, IDO1, CXCL9), myeloid lineage (CD14, CD163, ARG1, CD68, OLR1, NOS2, CD33), inhibitory immune ligands/receptors (HAVCR2, CTLA4, LAG3, PDCD1, CD274), immune modulators (ENTPD1, NT5E) and activating immune receptors (TNFRSF4, TNFRSF9, CD40, CD80, CD27, ICOS) were analyzed as previously reported.[45] The edgeR and limma packages were used to calculate the fold change of genes between groups.[26,46]

Several classic CAF markers were adopted from previous studies to estimate the relative abundance of fibroblasts in cancer samples, including PDGFRA, PDGFRB, ACTA2, THY1, PDPN, COL1A1 and FAP.[47-50] Clustering was performed in individual datasets, and the samples were further classified into high-, medium-, and low-infiltration groups using the ComplexHeatmap package.[51] To further confirm the clustering results, the principal components analysis (PCA) method was applied as previously desribed.[52] Inflammatory and myofibroblastic CAF features were also included to assess the internal characteristics. Comparisons of biological markers among different CAF infiltration groups are shown by the ggheatmap and ComplexHeatmap packages.

Functional analysis

The GSVA package was used for gene set variation analysis (GSVA).[53] The GSVA results were compared between the high- and low-CAF infiltration groups and are displayed. Gene set enrichment analysis (GSEA) was used to explore the biological functions of BGN and performed using GSEA 4.1.0.[54] Hallmark and gene ontology gene sets were obtained from the MSigDB Collections (http://www.gsea-msigdb.org/gsea). The TGF-β response signatures (TBRS) of T cells (T-TBRS), fibroblasts (F-TBRS), endothelial cells (End-TBRS) and macrophages (Ma-TBRS) were calculated and compared between different CAF infiltration groups as previously reported.[55]

Primary cell culture, RNA sequencing and mass spectrometry

Primary fibroblasts were isolated from immunohistochemically confirmed breast cancer samples obtained from surgery. For normal cancer-adjacent fibroblasts (NAFs), tissues >2 cm away from cancer were used. Primary fibroblasts were isolated as described in previous studies.[56] Briefly, tissues were digested with hyaluronidase (Sigma), collagenase type I (Sigma) and collagenase type III (Worthinton) at 37°C with agitation for 3–4 hours. The digested tissues were then incubated without shaking for 10 min, and the supernatant enriched with stromal cells was moved to new tubes and centrifuged at 250 g for 10 min. The supernatant was discarded, and the pellet was harvested and cultured in DMEM with 10% FBS. Fibroblasts were passaged within 5 population doublings after isolation. The RNA of the primary cells was extracted using TRIzol reagent (Invitrogen). More than 2 μg of high-quality RNA with a 260/280 absorbance ratio of 1.8–2.2 per sample was used for subsequent experiments. The rRNA-depleted RNA was prepared by the Ribo-off™ rRNA depletion kit (Vazyme, N406-01) and further applied to synthesize the first cDNA using Total RNA-seq Library Prep Kit (Vazyme, NR603). After library construction, fragments were further enriched by PCR amplification (NEST Biotechnology) and screened before sequencing. The Illumina NovaSeq 6000 sequencing platform was used for sequencing and raw reads were generated following the recommended protocol from the vendor. The transcripts per million (TPM) were calculated at the transcript level based on counts data.

The culture medium of primary cells was collected and centrifuged at 500 g for 10 min at 4°C and 10,000 g for 30 min at 4°C to remove cellular debris. The supernatant was further filtered using a 0.22 μm filter. Briefly, the protein was extracted, digested by trypsin and further applied for label-free mass spectrometry. The results were further searched in Maxquant software and quantified using LFQ method.

Immunohistochemistry

The samples were fixed with paraformaldehyde and embedded in paraffin. The samples were then sectioned into 3 μm slides. Antigen retrieval was performed using a pressure cooker with citrate buffer (CWBIO) for 10–15 min. Goat serum (ZSGB-BIO) was used to block nonspecific binding to the samples. The samples were further incubated with specific anti-biglycan (1:150, Proteintech) or anti-CD8 (ready to use, ZSGB-BIO) antibodies at 4°C overnight and a secondary antibody (ZSGB-BIO) at room temperature for 1 hour. The expression of protein was detected using DAB (Dako) following the manufacturers’ instructions. The density of infiltrating CD8 + T cells was estimated using HALO software (Indica Labs).

Western blot

For cell lines, cells were lysed using RIPA buffer (Beyotime) at 4°C directly. Phenylmethanesulfonyl fluoride (PMSF, Beyotime) was added to reduce protein degradation during extraction. Proteins were separated in SDS-polyacrylamide gels and transferred to PVDF membranes, and nonfat milk was used to block the nonspecific binding sites on the membrane. The membranes were incubated with primary antibodies against biglycan (1:1000, Proteintech) and GAPDH (1:3000, SAB) at 4°C overnight. The secondary antibody (1:5000, Bioss) was applied on the following day, and the reaction was detected using enhanced chemiluminescence solution (ECL, Affinity).

Statistical methods

The best cutoff values for specific markers in each cohort were determined using the survminer package. The survival package was used for Kaplan-Meier overall survival analysis, and the log-rank test was applied for comparison. The hazard ratio (HR) was calculated via univariate Cox regression. Immune signatures were divided into two groups according to the median value and calculated by Cox regression in. Student’s t-test or Wilcoxon rank-sum test were used for comparison of normally and non-normally distributed variables in unpaired groups, respectively. The paired Student’s t-test was performed for paired samples. Chi-square test and Fisher’s exact test were applied for comparison of clinical features. The Spearman method was applied for correlation analysis. All P values were two-tailed. Statistical analysis was performed using R software (Version 3.5.3, https://www.r-project.org).

Results

CAF marker-based classifier identified three clusters with different CAF infiltration in TNBC cohorts

The relative abundance of fibroblasts in TNBC was estimated by expression profile clustering using classic CAF markers reported in previous studies, including ACTA2, FAP, PDGFRA, PDGFRB, PDPN, THY1 and COL1A1. The expression level of the above markers in CAFs of breast cancer was examined using single-cell sequencing datasets (Figure S1). Samples from the TCGA-TNBC, METABRIC-TNBC, GSE25066-TNBC, GSE21653 and GSE103091 cohorts were analyzed. The clustering results of TCGA-TNBC, METABRIC-TNBC and GSE25066-TNBC are shown in Figure 1a, S2a and S2c, and patients were divided into different CAF infiltration groups, namely, low, medium and high infiltration. The distributions of clinical or biological features, such as the ESTIMATE stromal score, immune score, mutation type, age at diagnosis, TNM stage, subtypes, amongst others, were displayed under the heatmap. PCA plots were employed to verify the clustering efficacy in the cohorts (Figure 1b, S2b, S2d). To further examine whether the clustering groups reflect the relative abundance of CAFs in the TME, we compared the CAF scores reported in previous studies and demonstrated an ascending tendency of MCP-counter and xCell scores in the low, medium and high CAF infiltration groups of representative cohorts (Figure 2a, 2d, 2g). Moreover, the expression of iCAF and myCAF markers was also compared. The results showed that lots of CAF markers were upregulated in the high and medium infiltration groups (Figure 2b, 2e, 2h). The results indicated that samples in these groups showed abundant infiltrating CAFs. The CAF infiltration clusters of all TNBC cohorts in this study are shown in Figure S3a. The overall survival did not show a significant difference among the CAF infiltration groups (Figure S3b-f).

Figure 1.

Clustering of classic CAF markers in the TCGA-TNBC cohort. a, Heatmap showing the expression level of CAF markers and clustering of CAF infiltration in the TCGA-TNBC (N = 116) cohort. Bars under the heatmap display the distributions of clinical features and ESTIMATE scores. b, Scatter plots showing the PCA result of the TCGA-TNBC cohort based on the classic CAF markers (x and y axes). External circle, 98% confidence interval; Internal circle, 95% confidence interval; Eclipse, 80%. Arrows, the tendency of marker profiling; low, medium and high infiltration levels are shown as green, blue and red, respectively.

Fibroblastic and immune features in CAF infiltration groups. Box plots comparing the fibroblastic scores from MCP-counter (upper) and xCell (lower) in different infiltration groups of the METABRIC-TNBC (a), TCGA-TNBC (d) and GSE25066-TNBC (g) cohorts. P values from Student’s t-test. Heatmap showing the expression profile of myCAF and iCAF features in different infiltration groups of the METABRIC-TNBC (b), TCGA-TNBC (e) and GSE25066-TNBC (h) cohorts. Heatmap showing the fold change of the immune features in high and medium groups compared with the low CAF infiltration group (high vs. low, medium vs. low) in the METABRIC-TNBC (c), TCGA-TNBC (f) and GSE25066-TNBC (i) cohorts. P values are reported as: ns, nonsignificant; *, P < .05; **, P < .01; ***, P < .001; ****, P < .0001.

Continued.

Continued.

Correlation between CAFs infiltration and immune features

Previous studies have indicated the crucial role of CAFs in immune modulation of the microenvironment. Here, we sought to examine the correlation between CAFs and the immune landscape. First, we calculated the fold change of genes between different infiltration groups (high vs. low, medium vs. low) in METABRIC-TNBC, TCGA-TNBC and GSE25066-TNBC cohorts and the immune-related genes were shown (Figure 2c, 2f, 2i). The results indicated that some of them were upregulated in the high and medium infiltration groups of different cohorts, such as CD14, TGFB3, ENTPD1, NT5E and BMP1. Other upregulated markers were not consistent in different cohorts, such as LAG3, CD33, CD68, CXCL9, and IL6. Second, the GSVA scores of hallmark signaling pathways in cancer were compared between the high- and low-infiltration groups (Figure 3a-c). Some of the immune-related signaling pathways were significantly enriched in the high infiltration group, such as TGF-β, IL2-STAT5, complement, TNF-α-NF-κB and IL6-JAK-STAT3 signaling. Finally, as CAFs are one of the producers of TGF-β, we compared the expression levels of the TGF-β response signatures of fibroblasts (F-TBRS), T cells (T-TBRS), macrophages (Ma-TBRS) and endothelial cells (End-TBRS). Most signatures showed an ascending tendency from the low- to high-infiltration groups (Figure S4).

Figure 3.

Correlation between CAFs and hallmark signaling and immune cells. a-c, Bar plots showing the GSVA scores of cancer hallmark signaling in the high infiltration group compared with the low CAF infiltration group. Significantly upregulated and downregulated signaling pathways are shown in red and green, respectively. Nonsignificant signalings are shown as gray. d-f, Network showing the correlation among CAFs and CIBERSORTx-derived immune cells. Significantly positive and negative correlations are shown as red and blue lines, respectively. The color and size of the nodes indicate the type of cells and P values from Cox regression. Prognostic signatures for overall survival are marked with dark dots in the nodes.

Based on the results above and previous studies, CAFs might exert their function by affecting tumors and immune cells. We evaluated the correlation between previously reported CAF scores and CIBERSORTx immune cells (Figure 3d-f). The CAF scores estimated by xCell, EPIC and MCP-counter were analyzed. CAFs were negatively correlated with cytotoxic cells, including CD8 + T cells and activated NK cells. In the METABRIC-TNBC cohort, MCP-counter-CAF and xCell-CAF were negatively correlated with activated NK cell, respectively. But only MCP-counter-CAF was negatively correlated with CD8 + T cell. In the TCGA-TNBC cohort, EPIC-CAF, MCP-counter-CAF and xCell-CAF were all negatively correlated with CD8 + T cell. But only EPIC-CAF and MCP-counter-CAF were negatively correlated with activated NK cell. In the GSE25066-TNBC cohort, only a negative correlation between xCell-CAF and CD8 + T cell was observed. Moreover, CAFs might exhibit significantly positive correlations with M2 macrophages and this result was only observed in METABRIC-TNBC and GSE25066-TNBC cohorts. Furthermore, xCell-CAF and EPIC-CAF were risk factors for overall survival in METABRIC-TNBC and TCGA-TNBC cohorts, respectively.

BGN is upregulated and mainly expressed in stromal CAFs

Here, we sought to investigate effective biomarkers derived from CAFs and their corresponding prognostic value. Normal cancer-adjacent fibroblasts and cancer-associated fibroblasts were isolated from tissue samples after surgical resection. RNA sequencing and label-free mass spectrometry were further applied to explore the transcriptomic expression profile and differentially expressed genes (DEGs) in the NAF and CAF groups (Figure 4a). The top 40 expressed DEGs in CAFs are shown in Figure 4b. The label-free mass spectrometry identified the top 40 secreted proteins (Figure 4c). BGN, an extracellular protein, was upregulated in CAFs compared to NAFs. Previous studies indicated the roles of BGN in therapy resistance and immune activity.[57-60] Overexpression of biglycan induced resistance to 5-FU and rapamycin in cancer cells via NF-κB signaling pathway. Besides, most mechanistic studies mainly focused on the interaction between biglycan and toll-like receptors (TLR2/4).[61] However, the mechanism underlying how biglycan affects tumor microenvironment in other ways is still poorly understood despite its abundance in TME. The TPM values of BGN from RNA sequencing were compared (Figure 4d). We further compared the BGN level in the TCGA-TNBC cohort and observed a significantly higher level of BGN in tumor tissues (Figure 4e). We further verified our results at the protein level by Western blotting and observed a significantly higher level of biglycan protein in CAFs of patient 1, 2 and 3 (Figure 4f).

Figure 4.

BGN is upregulated in cancer-associated fibroblasts. a, Illustration showing the procedure of RNA sequencing and mass spectrometry analysis for tissue-derived fibroblasts and culture medium, respectively. Representative image of primary fibroblasts is shown. b, Heatmap showing the top 40 expressed genes upregulated in CAFs. c, Heatmap showing the top 40 secreted proteins identified by label-free mass spectrometry in culture medium derived from CAFs and the corresponding NAFs. d. Box plot comparing the RNA level of biglycan in NAFs and CAFs. P value from paired Student’s t-test. e, Box plot comparing BGN expression in normal and tumor tissues in the TCGA-TNBC cohort. P value from Student’s t-test. f, Western blot showing the expression of biglycan protein in fibroblasts isolated from normal and cancer tissues of TNBC. g, Violin plot comparing the BGN expression level in the epithelial and stromal components derived from microdissection (GSE88715). P value from paired Wilcoxon rank-sum test.

To further confirm the expression profile of BGN in different tumor components, we analyzed its expression pattern using microdissection and single-cell sequencing datasets. Microdissection analysis (GSE88715) indicated a higher level of BGN expression in the stromal part of the tumor (Figure 4g). Moreover, BGN expression showed an ascending tendency from the low to high CAF infiltration group described above (Figure S5a). The expression level of BGN was also positively correlated with classic CAF markers in each cohort (Figure S5b).

To further confirm its expression feature in TNBC, we sought to analyze the BGN expression level in TNBC scRNA-seq datasets. BGN was mainly expressed in the overall fibroblasts of TNBC cohorts (Bassez A, et al.; GSE118389; Figure 5a, 5b). In another dataset with more detailed cell clusters, BGN was mainly expressed in CAFs and perivascular-like (PVL) subpopulations (Wu SZ, et al.; Figure 5c-e). The BGN expression level was significantly correlated with the scores of fibroblasts and endothelial cells estimated by MCP-counter and xCell in these cohorts, respectively (Figure 6a-c-). This result provided additional evidence for the derivation of biglycan. These results further indicated that BGN might be a CAF-specific biomarker in TNBC and exert its function in the extracellular matrix.

Figure 5.

BGN expression in single-cell RNA sequencing datasets of TNBC. a-c, Violin plots showing the expression level of BGN in cells of scRNA-seq datasets (Bassez A ; GSE118389; Wu SZ). d, Scatter plot showing the cell clusters in c (Wu SZ). e, Scatter plot showing the distribution of cells expressing a high level of BGN in c.

Figure 6.

The correlation of BGN and CAF scores and pathway enrichment. a-c, Scatter plots showing the correlation between BGN and CAFs and endothelial cells. P and r values from Spearman correlation analyses. d-f, Enrichment scores of the representative enriched signaling pathways in the high-BGN group compared with the low-BGN group. Enrichment Score from the Gene Ontology enrichment.

Biglycan correlates biological features and clinical outcomes

The BGN gene encodes biglycan, an extracellular soluble protein that might exert its function via intercellular contact. Here, we sought to explore its role in tumor development and its prognostic value in clinical practice. GO signaling enrichment was performed using BGN-correlated genes, and the representative results are shown in Figure 6d-f. Representative signaling pathways involved in cancer progression and TME modulation were significantly enriched in METABRIC-TNBC, TCGA-TNBC, GSE25066-TNBC cohorts (FDR<0.05). The high BGN level indicates the extension of new blood vessels from the existing vessels (SPROUTING_ANGIOGENESIS), lymph vessel formation (LYMPHANGIOGENESIS), extracellular matrix remodeling (COLLAGEN_CONTAINING_ECM) and promotion of tumor cell metastasis (POSITIVE_REGULATION_OF_EMT). Furthermore, BGN was also correlated with the formation of an immunosuppressive microenvironment (TGF_BETA_RECEPTOR_SIGNALING). These results indicated that biglycan might be associated with tumor angiogenesis, TME remodeling and tumor metastasis in TNBC. Hence it might serve as a protumor and immunosuppressive factor in these cohorts. Therefore, we further analyzed the association between BGN and the immune signatures. First, the correlation between BGN and immune-related genes was analyzed (Figure 7a). BGN was positively correlated with some immune-inhibitory molecules, including TGFB1, VEGFA, VEGFB and CD276, especially in METABRIC-TNBC and GSE21653 cohorts. However, we did not observe significant association with classic checkpoint molecules except for CTLA4 and CD274 (PD-L1) in GSE103091 and BTLA in GSE21653. These results imply that BGN might affect immune cell infiltration or modulate immune activity. Hence, we calculated the correlation between BGN expression and ESTIMATE scores (Figure 7b). The BGN was significantly positively correlated with the stromal scores in each cohort. We observed significant correlations between BGN and the immune scores in METABRIC-TNBC (negative, P < .05), GSE21653 (negative, P < .05) and GSE103091 (negative, P < .1). Finally, we examined the correlation between BGN and immune components. BGN was significantly correlated with both CD8 + T cells (negative, P < .05) and activated NK cells (negative, P < .05) in the METABRIC-TNBC and TCGA-TNBC cohorts, respectively. Moreover, BGN was correlated with M0 macrophages in all cohorts (positive, P < .05) except for GSE21653 and M2 macrophages in METABRIC-TNBC, GSE25066-TNBC and GSE103091 (positive, P < .05) (Figure 7c).

Figure 7.

The correlation between BGN and immune features. a, Heatmap showing the correlation between BGN and immune-related genes. b, Bubble plot showing the correlation between BGN and ESTIMATE stromal and immune scores. c, Bubble plot showing the correlation between BGN and CIBERSORT-derived immune cell scores. Positive and negative correlations are shown as red and blue, respectively. Correlation coefficients and P values from the Spearman method.

Since BGN plays an essential role in the biological process of cancer, we further sought to examine its clinical relevance. First, we observed a significantly lower level of BGN in TNBC than in non-TNBC samples of the TCGA and METABRIC cohorts (Figure S6a). BGN expression was significantly upregulated in tumors compared with paired normal tissues in the TCGA-TNBC cohort (Figure S6b). Second, the Cox regression model was conducted according to the best cutoff value for overall survival and indicated that BGN was a poor prognostic predictor in the METABRIC-TNBC, GSE103091 and GSE21653 cohorts (Figure S6c). Besides, FAP predicted poor prognosis in TCGA-TNBC and GSE21653 cohorts. PDGFRB and COL1A1 predicted poor outcomes in TCGA-TNBC and GSE103091 cohorts. However, PDGFRA was predictive of favorable outcomes in cohorts except for TCGA-TNBC and PDPN predicted favorable prognosis in METABRIC-TNBC and GSE25066-TNBC cohorts. Kaplan-Meier analysis demonstrated that a higher BGN expression level predicted worse overall survival outcomes in each cohort (Figure 8a). Finally, we also validated the prognostic value and immunosuppressive function of biglycan in the SYSUCC cohort. Patients in the high-biglycan group showed a lower level of total infiltrating CD8 + T cells (Figure 8b, 8c). Survival analyses demonstrated that biglycan might serve as a poor prognostic marker in TNBC (Figure 8d, Table S2). The univariate and multivariate Cox regression model for the SYSUCC cohort indicated that biglycan was an independent risk factor for overall survival (Table S3).

Figure 8.

Survival analysis of BGN in datasets and the SYSUCC cohort. a. Kaplan-Meier plots comparing the overall survival between the high-BGN and low-BGN groups according to the best cutoff value. P values from log-rank tests. High-BGN and low-BGN groups are shown as red and blue. b, Representative IHC staining images of stromal biglycan and CD8 in TNBC of SYSUCC (N = 100). Scale bar, 100 μm. c, Boxplot comparing the density of CD8 + T cells in high-biglycan and low-biglycan groups. P value from Student’s t-test. d, Kaplan-Meier plot comparing the overall survival between high-biglycan and low-biglycan groups. P value from log-rank tests. High-biglycan and low-biglycan groups are shown as red and blue.

Discussion

Recent studies have indicated the crucial roles of cancer-associated fibroblasts in breast cancer.[15] Many of them focus on the heterogeneity and corresponding biological features of different CAF clusters in cancer.[16,17,62] However, the detailed mechanism underlying how CAFs affect the TME has been a prominent topic in recent years.

Here, we have included the expression profile of the tumor bulk and single-cell RNA sequencing to explore the characteristics of fibroblasts in triple-negative breast cancer. We evaluated the relative abundance of fibroblasts in TNBC and further identified an extracellular secreted protein, biglycan, as a biomarker for CAFs and a predictor for poor prognosis in TNBC.

CAFs have been reported to be a tumor-promoting component of the stroma in most cancers.[10] In previous studies, fibroblasts exerted their function by producing excreted factors, remodeling the extracellular matrix, influencing cancer cell metabolism and direct cell-cell interactions. For example, CAF-derived PLGF and BDNF promote NR4A1 expression in triple-negative breast cancer cells and enhance invasiveness.[11] In addition, CCL2/7 are produced by lung-resident fibroblasts to modify cancer cell metabolism by activating cholesterol synthesis in TNBC cells in the lung and promoting angiogenesis.[63] Furthermore, CAFs serve as leading cells for cancer cells during cell migration.[64] Fibroblasts pave the way for subsequent malignant cells and lead to tumor invasion by direct cell-to-cell contact. Here, we sought to explore the effect of the relative abundance of CAFs in TNBC. The CAF infiltration level affects tumor cells and the TME with different CAF, immune and tumor cell expression profiles. Although some immune markers show an ascending tendency from low to high CAF infiltration levels (CD14, NT5E and ENTPD1, etc.), some other markers do not. This is consistent with a previous study which indicates that fibroblast-enriched samples might be clustered into fibrotic (F) and immune-enriched/fibrotic (IE/F) groups with distinct features.[65] Furthermore, the simple relative abundance of CAF infiltration did not necessarily seem to be a prognostic predictor. This is consistent with previous studies showing that CAFs are highly heterogeneous and different subpopulations display various features. Tumor-infiltrating fibroblasts are commonly divided into inflammatory CAFs (iCAFs) with lower α-SMA and higher cytokine production, and myofibroblastic CAFs (myCAFs) with higher levels of α-SMA.[34] Recently, a novel cluster of antigen-presenting CAFs (apCAFs) was identified, and CAFs were confirmed to directly participate in immune reactions.[44] The reason why high infiltration group exhibited a higher level of BGN might be partially due to a larger amount of overall CAFs. But it is not capable of predicting the dominant role of BGN-expression CAFs. These results indicate that the heterogeneity of CAFs is more important than simply cell abundance. A single-cell sequencing cohort of TNBC with a large population and active follow-up data will help to clarify the role of the amounts of overall CAFs and specific dominant subclusters in predicting clinical outcomes in future studies.

Furthermore, internal heterogeneity in the cellular population poses a new challenge in modern medicine. In recent years, single-cell sequencing technology has helped to explore single-cell expression profiles and identify internal heterogeneous clusters of cells. Infiltrating CAF-S1 is enriched in TNBC, and further clustering by single-cell sequencing reveals 8 subclusters with distinct features, in which Clusters 0/3 are linked to immunotherapy resistance in cancer.[17] Together, these results imply that specific CAF subpopulations might truly matter instead of simply the number of infiltrating fibroblasts in TNBC. Hence, it is more important to explore the expression profile of CAFs.

As CAF-derived proteins represent the most common way for intercellular crosstalk and might serve as biomarkers for cancer, we used RNA sequencing and mass spectrometry technology, and identified the coding gene for biglycan, BGN, to be mainly expressed in CAFs. The biglycan protein belongs to the SLRP family and plays a role in the extracellular matrix (ECM).[66] Interestingly, although BGN was found to be upregulated in CAFs compared to NAFs, it also seemed to be expressed in endothelial cells. A previous study reported that BGN was approximately 100-fold more abundant in tumor endothelial cells than in corresponding normal endothelial cells.[67] In another study, the results from scRNA-seq indicated that BGN is also expressed in perivascular-like subclusters except for CAFs.[34] However, this is not contradictory to our findings, as endothelial cells and pericytes are commonly recognized as crucial origins for CAFs and transform into CAFs under certain circumstances.[14] Furthermore, the number of fibroblasts was far greater than that of endothelial cells in scRNA-seq datasets, which indicates that CAFs are the major origin of biglycan proteins. BGN expression correlates with poor overall survival, cancer biological processes and microenvironmental components. Biglycan is a nonfibrillar component in the ECM, suggesting that extracellular biglycan might not only constitute the ECM but also exert other functions. BGN encodes a secreted protein that could bind to receptors on immune cells and it might also potentially interact with tumor cells, pericytes or endothelial cells.[60] The Gene Ontology enrichment in this study implies the potential role of BGN in tumor vasculature and our findings are consistent with previous studies that CAFs promoted tumor angiogenesis mainly by releasing secreted factors, such as SDF1, MMPs and VEGFs.[68,69] Biglycan is a vital component in bone and might promote the proliferation of osteosarcoma cells through an LPR6/β-catenin/IGFR-IR axis.[70] Other studies indicate that biglycan might enhance drug resistance in cancer cells.[58,59] These studies imply that biglycan promotes cancer progression by cell crosstalk. In previous studies, soluble biglycan was found to interact with CD14 and Toll-like receptor 2/4 (TLR2/4), leading to the activation of downstream NF-κB and ERK1/2 signaling.[57,60] TLR2/4 are membrane receptors for various immune cells and modulate the immune response in the TME.[71] Previous studies have confirmed the role of TLR4 in the polarization of tumor-associated macrophages and the conversion of CAFs, tumor-associated dendritic cells and myeloid-derived suppressor cells.[60,61,72] On the basis of previous studies, we might hypothesize that biglycan exerts its function as a secreted ligand and modulates the downstream signalings by interacting with corresponding receptors on other cells. These results suggest that biglycan might play a crucial role in immune modulation and tumor angiogenesis.

We also realize the limitations of our study. First, the relative abundance of fibroblasts was evaluated by a clustering method using classic CAF markers with higher specificity instead of all well-known fibroblast markers. Comprehensive clustering method and single-cell sequencing data were applied to enhance the specificity, but the similarity of expression profiles in CAFs and vascular cells might still confound our analysis. Besides, the results are only comparable within each cohort. Second, sampling bias might occur because the fraction of the stromal part varies in different samples due to the internal heterogeneity of the tumor bulk. Bulk sequencing might not truly reflect the overall landscape of tumors. Third, the mechanical and chemical effects during isolation of primary cells from tissues and in vitro cells culture might change the phenotype of fibroblasts, leading to detection bias in the experiments. Finally, the mechanism underlying how biglycan affects the TME was not explored in this study. Further investigations about how biglycan affects tumor cells and immune components, such as CD8 + T cells and macrophages, will advance our understanding of the roles of CAFs in TNBC.

Notably, we observed that biglycan was upregulated and secreted in cancer-associated fibroblasts compared with normal cancer-adjacent fibroblasts. The expression of BGN was significantly correlated with overall clinical survival and biological processes in TNBC, suggesting that biglycan might play a crucial role in cell-to-cell crosstalk between CAFs and cancer cells. Moreover, we confirmed the prognostic value of biglycan in TNBC using the cancer tissue microarray. Interestingly, BGN was not significantly correlated with classic checkpoint molecules in most cohorts. But we observed an evident positive correlation between BGN and TGF-β, a well-known immunosuppressive molecule, indicating that BGN might exert its effect on immune activity in a certain way. Besides, BGN might participate in tumor angiogenesis which is also a prominent target to synergize with immunotherapy.

In conclusion, we identified that CAF-derived biglycan is crucial for TNBC progression. The upregulation of biglycan and the mechanism underlying how biglycan exerts its function in the TME may serve as a promising diagnostic biomarker and may provide more promising strategies for cancer treatment.

Click here for additional data file.