The fragments per kilobase of per million (FPKM) of breast cancer transcriptome were downloaded from TCGA database (https://portal.gdc.cancer.gov). A reference gene set of cancer extracellular matrix associated (C-ECM) genes including 30 significant upregulated genes and 28 significant downregulated genes in pan-cancer tissues was derived from Ankur Chakravarthy et al (14). ECM-up and –down scores were carried out by using ssGSEA method of R software Gene Set Variation Analysis (GSVA) package. According to the clinical data of breast cancer samples in the TCGA, the associations between ECM scores (ECM-up and –down scores) and breast cancer molecular subtypes were evaluated. Then the ECM up and -down scores are classified according to the median. We take the intersection of the samples that are greater than the median in the up-regulation score and those are less than the median in the down-regulation score, and define it as the ECM-high group. At the same time, we define the ECM-low group by taking the intersection of the samples that are less than the median in the up-regulated score and those that are greater than the median in the down-regulated score, while those samples, in which ECM-up scores higher than the median coexist with ECM-down scores lower than the median or ECM-up scores lower than the median coexist with ECM-down scores higher than the median, were defined as ECM-median group.

We compared the overall survival (OS) and disease free survival (DFS) of breast cancer patients in ECM-high, -median and –low subtype by R software “survival” package. The log-rank test was performed to test the survival differences. Kaplan-Meier curves were used to visualize the OS and DFS of different ECM subtypes. Univariate and multivariate Cox regression analysis regarding on OS was performed to determine the hazard ratio (HR) of each ECM subtype by using the “survival” package of the R software. The sensitivity and specificity of the ROC curve were conducted to evaluate the prognostic value of the ECM subtypes.

The “edgeR” package calculation in R software was used to perform a differential analysis of the mRNAs of ECM-high and –low groups according to the cutoff screening of |log2FC| > 1 and p <0.05. Gene ontology (GO) terms including biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were analyzed by R “clusterProfiler” package to perform biological functions of differentially expressed genes (DEGs) between ECM-high and –low groups. To explore the possible signaling pathways DEGs enriched, KEGG enrichment analysis was carried out utilizing the “clusterProfiler” package of R software with a statistical threshold of p < 0.05.

The relationship between ECM dysfunction and microenvironment status was analyzed on TIMER2.0 (http://timer.cistrome.org/). This freely available web provides more robust estimation of immune infiltration levels for TCGA or user-provided tumor profiles using six state-of-the-art algorithms. Microenvironment Cell Populations-counter (MCP-counter) algorithm was selected for quantification of the absolute abundance of eight immune and two stromal cell populations in each group (16).

In order to clarify the characteristics of the ECM corresponding to each ECM phenotype, we performed Masson’s trichrome staining (Solarbio, China) and Sirius red/Fast Green staining (Solarbio, China) to assess the ECM characteristics of breast cancer patients of different nodal status and pathological stages according to the manufacturer’s instructions. The density of collagen was estimated using red-stained collagen area of Sirius red/Fast green staining and blue-stained collagen area of Masson’s trichrome using Image J in ten randomly selected microscopic fields per specimen. The degree of collagen fiber alignment was quantified using an image J plugin Orientation J.

All the animal studies conformed to the guidelines set by the Institutional Animal Care Committee of Xi’an Jiaotong University. Collagen fibers were harvested from rat tail tendons and labeled with fluorescein isothiocyanate (FITC) (Yeasen, China) as described previously (17, 18). The human breast cancer cell line MDA-MB-231, the mouse embryo fibroblast cell line NIH/3T-3 and human embryonic lung fibroblast cell line MRC-5 were obtained from American Type Culture Collection (ATCC, Manassas, VA) and cultured in DMEM supplemented with fetal bovine serum (FBS; Gibco, Grand Island, NY, USA), penicillin, and streptomycin (Sigma, St. Louis, MO, USA). Primary human T cells were obtained as a gift from Dr Jinteng Feng and Jurkat T cells were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). The culture medium of primary T and Jurkat T cells were RPMI 1640 (Hyclone, Logan, UT, USA) supplemented with FBS, penicillin, and streptomycin. The supernatant of MB-231 cells was added to fibroblasts and T cell culture medium at a 1:1 ratio for 12 hours during 3D culturing. Primary T and Jurkat T cells were stimulated by Dynabeads™ Human T-activator CD3/CD28 (ThermoFisher, 11161D) for 24h before co-cultured with the supernatant of breast cancer cells. T cells and fibroblasts were capsulated in low and high density gels contained 1mg/ml and 4 mg/ml type I collagen at a concentration of 6×10⁵ cells/gels respectively. Both cell types were used between passages 2–20 and all cells were maintained at 37°C in a humidified 5% CO2 atmosphere.

Cells capsulated in collagen gels were injected into a dog bone-like mold. 10 iron oxide microspheres with diameter of ~200 μm were placed into both edges of collagen suspension. To obtain aligned collagen fibers, two NdFeB magnets (6 mm × 1.5 mm) were used to apply 10% uniaxial strain to the gels for 24h.

Cells were encapsulated in collagen gels with different collagen concentration or fiber organization. Then rinsed, fixed in 4% formaldehyde and permeabilized for 10 min with 0.5% Tritonx-100, and incubated in 1% BSA for 1 hour followed by incubation with 1:500 dilution of anti-alpha smooth muscle actin antibody (ab7817, Abcam, 4°C, overnight) and fibronectin monoclonal antibody (14-9869-82, invitrogen, 4°C, overnight), then incubated with secondary antibody (Alexa Fluor^® 647, Invitrogen, room temperature, 2hs), and 4′,6-Diamidino-2-phenylindole (DAPI, Invitrogen) was used to stain cell nuclei. Confocal fluorescence microscopy was performed using Olympus FV300. The fluorophores were excited by 405, 488 and 647 nm laser lines.

Supernatant samples were collected 48h after encapsulation of cells in collagen gels and assayed for IL-2 levels using a commercial IL-2 ELISA Reagent Kit (Fisher scientific) following the product protocol. Measure absorbance on an ELISA plate reader set at 450nm with microplate reader (Thermo Fisher) and use the provided Standard Diluent to prepare standard curve serial dilutions.

All data were analyzed by R version 3.6.3 or Graphpad Prism 8, and all cell culture experiments were performed at least three times independently. These results were presented as mean ± standard deviation (SD). Student’s two-sided t-test was used to compare the differences between two groups. Differences in survival between different risk groups were compared by Kaplan-Meier curves followed by log-rank test. P < 0.05 was considered as statistically significant.

Initially, to study whether ECM characteristics showed significant effects in invasion and prognosis of breast cancer, a total of 1109 breast cancer samples and 113 paracancerous samples from TCGA were obtained. The entire flow chart of the research is presented in Figure 1 , and the clinical information of patients in TCGA database was list in Table 1 . Upon summarization using single sample Gene Set Enrichment Analysis (ssGSEA) scores to evaluation ECM characteristics, C-ECM-up and -down scores show broad variation across breast cancer molecular subtypes ( Figures 2A, B ). Also, we found that triple negative breast cancer (TNBC) which is considered as the most aggressive molecular subtype showed highest ECM-up scores and lowest ECM-down scores ( Figures 2A, B ). Specifically, compared with the paracancerous samples, cancer samples had higher ECM-up scores and lower ECM-down scores in paired and unpaired breast cancer tissues ( Figures 2A – D ).

Subsequently, to estimate whether ECM characteristics contribute to the prognosis of breast cancer, breast cancer patients were divided into 154 ECM-high, 662 ECM-median and 155 ECM-low patients according to ECM-up and –down scores. Kaplan-Meier curves showed that the samples in the ECM-high group exhibited worse overall survival (OS) and disease free survival (DFS) than those in the ECM-median and -low groups, while ECM-low group exhibited the best OS and DFS, indicating the association between ECM dysregulation with breast cancer prognosis ( Figures 3A, B ). Beside, to further explore the prognostic value of the ECM score in each molecular subtypes of breast cancer, we compared the OS and DFS between different ECM subgroups for different molecular subtypes separately. The results of DFS and OS in luminal breast cancer patients were consistent with the general samples ( Figures S1A, B ), however, in TNBC subtype, the ECM-medium group showed the best OS and DFS trends among the three groups with no statistically significant p-values ( Figures S1C, D ). Considering the small number of patients in the TNBC subtype, we reclassified these patients into ECM-high and -low groups, then the OS and DFS results began to show a similar trend with ECM-high subgroup showed the worse OS and DFS, while the p-values still not statistically significant ( Figures S1E, F ).

We also performed univariate ( Figure S2 ) and multivariate Cox regression analyses regarding OS based on ECM group with TNM stage and age of onset to examine the prognostic impact of this dysregulation. ECM-high group were significantly associated with poor prognosis ( Figure 3C , hazard ratio (HR)=2.404, p=0.022). In order to compare the sensitivity and specificity of risk score on the prognosis of patients with breast cancer, time-dependent receiver operating characteristics (ROC) analysis was performed. The area under the ROC curve (AUC) of the ECM group was 0.755 ( Figure 3D ), suggesting the ECM dysregulation prognostic signature for breast cancer was highly reliable. These results together identify ECM dysregulation as a key player in the prognosis of breast cancer.

Since our data suggest that the ECM dysregulation was associated with the prognostic of breast cancer, we try to select possible pathways under this dysregulation. To do this, the mRNA expression profiles were compared across the best prognosis (ECM-high) and the worse prognosis (ECM-low) groups. A total of 793 differentially expressed genes including 403 upregulation and 390 down regulation genes were identified ( Figure 4A ). Subsequently enrichment factor of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses referred to the significance of specific functions. In the GO analysis, DEGs were significantly enriched in ECM associated process such as extracellular structure organization (GO:0043062), collagen-containing extracellular matrix (GO:0062023), extracellular matrix structural constituent (GO:0005201) and so on ( Figure 4B ). Notably, enrichment for inflammatory processes and adaptive immune responses such as humoral immune response (GO:0006959), leukocyte migration (GO:0050900), T cell activation (GO:0042110) immunoglobulin complex (GO:0019814), were also identified ( Figure 4B ). For the KEGG pathway analysis, DEGs were significantly enriched in ECM-receptor interaction, PI3K-Akt signaling, IL-17 signaling, human T-cell leukemia virus 1 infection pathway ( Figure 4C ).

As we identified ECM associated and immune response pathways, we then asked whether those two process were related. Given the previously identified role of CAFs in determining the composition and organization of ECM (19, 20), as well as T cell infiltration and activation in clinical efficacy of immunotherapies (21), we then try to evaluate the association between CAF induction, T cell infiltration and ECM dysregulation. To further support this hypothesis, Micro-environment Cell Populations-counter (MCP-counter) algorithm were conducted by TIMER to examine the abundance of T cell infiltration and CAF populations (16). Importantly, samples in the ECM-high group exhibited lower infiltration of T cells (1.153 vs 2.136, p=0.0001) and higher quantification of CAFs (232.6 vs 207.2, p=0.0343) ( Figures 4D, E ). To further confirm the results, two sets of T cell function and CAF activation biomarker genes were generated. Consistently, T cell function biomarkers were significantly lower in the ECM-high group, while CAF biomarkers were significantly higher ( Figures 4F – I ). Collectively, these results indicate that ECM characteristics affect activities of T cells and CAFs.

To understand the corresponding pathological characteristics of ECM dysregulation in human breast cancer samples, we performed Masson staining and Sirius red/Fast green staining on 30 breast cancer samples to evaluate the characteristics of collagen fibrils that constitute the main component of ECM. the clinical information of the real breast cancer patients was list in Table 1 . We compared the collagen density and organization of breast cancer patients with nodal status negative and positive as well as grade I-II and grade III ( Figures 5A, B ). Interestingly, the results showed increased collagen fibril alignment in the nodal status positive ( Figure 5C ) and grade III groups ( Figure 5E ) contrast to nodal negative or grade I-II breast cancer tissues. In addition, the density of collagen showed higher levels in the nodal status positive ( Figure 5D ) and grade III groups ( Figure 5F ). These results suggest an essential role of ECM characteristics in breast cancer invasion and metastatic dissemination. Give that ECM-high groups showed worse prognosis, taken together, these results strongly support increasing ECM alignment and density promoting invasion and predict worse prognosis in breast cancer.

Considering that the combined results suggests breast cancers with different invasiveness and prognosis showed different ECM characteristics with specific immune status and CAF activity, we tested the hypothesis that ECM characteristics would affect the activity of T cell activation and CAF induction. To evaluate the potential effects of CAF induction and T cell activation by collagen organization, we visualized the α-SMA and fibronectin expression of MRC-5 and NIH/3T-3 cells as well as the secretion of IL-2 of primary T cells and Jurkat T cells in 3D random and aligned collagen matrix supplement with supernatant samples of breast cancer cell MB-231. Visually, through magnetic stretching like Figure S3A , collagen fibrils with magnetic stretching remained aligned, whereas the control group remained random organized ( Figures 6A, B ). The α-SMA and fibronectin expression was significantly higher in the aligned collagen matrix ( Figures 6C – E and Figures S3C, D ). Conversely, the IL-2 levels in activated primary T cells and jurkat T cells exhibited significant decrease in the aligned group ( Figure 6F ). Remarkably, as activated primary T cells and jurkat T cells secreted IL-2 at levels significantly higher than the assay range, we detected IL-2 levels in a 10-fold dilution ( Figure S3B ). These observations indicate that increased ECM alignment significantly promote CAF induction and inhibit T cell activation, thus increase breast cancer invasiveness.

To further investigate the CAF induction and T cell activation regulated by another ECM characteristic, ECM density, we compared MRC-5 cells, NIH/3T-3 cells, primary T cells and Jurkat T cells cultured in 3D high- and low-density collagen matrix (4mg/ml vs 1mg/ml). We observed that α-SMA and fibronectin expression were significantly upregulated in high-density groups ( Figures 7A – C and Figures S4A, B ). Consistent with the published study (15), the increased collagen density was accompanied by decreased secretion of IL-2 ( Figure 7D ). Altogether, these results suggest that the T cells acquire a decreased activation in the high density than low density of collagen gels. This is consistent with the TCGA database analysis ( Figure 4G ).

Altogether, our study reveals that breast cancer patients with high ECM scores have increased ECM density and aligned organization, which promotes breast cancer invasion and metastasis via immunosuppression of T cells as well as enhancing CAF induction ( Figure 7E ).