MDA-MB-231 and SUM-149 were purchased from the American Type Culture Collection (ATCC, USA), and MDA-MB-231 was cultured in high glucose DMEM medium (PM150210, Procell) supplemented with 10% fetal bovine serum (164210, Procell) and 1% penicillin and streptomycin (PB180120, Procell). SUM-149 was cultured in Ham’s F12 medium (L410KJ, BasalMedia) supplemented with 5% fetal bovine serum, 5µg/ml of insulin (UNITED LABORATORIES), and 1µg/ml of hydrocortisone (R011855, RHAWN). MDA-MB-231-Re epirubicin-resistant cell lines were established as previously described (11). All cells were cultured in 5% CO₂, at 37°C.

MDA-MB-231-Re cells were seeded in six-well plates and incubated for 24 h. Lentivirus containing shRNA (shRUNX2, NC; MOI=10) was added when cell reached 20-35% confluency. Then, the supernatant was replaced with fresh medium after 10 h, followed by the observation of GFP fluorescence under a microscope at 72 h and stable line selection by puromycin addition. MDA-MB-231-Re cells with shRNA NC infection were used as a negative control.

Similarly, the RUNX2 over-expression (denoted as RUNX2) and the negative control (denoted as Vector) SUM-149 cell lines were generated by utilizing the aforementioned lentiviral strategy.

For RUNX2 interference with the specific siRNAs (siRUNX2#1, siRUNX2#2) (GenePharma), MDA-MB-231 cells were cultured in 60-mm dishes for 24 h and the medium was replaced with fresh medium containing 200 pmol siRNA and 10µl Lipo6000™ Transfection Reagent (C0526, Beyotime) for 5 h when cells reached 30-50% confluence. Then, the siRNA-Transfected MDA-MB-231 cells was cultured in fresh medium, and the efficiency of interference for RUNX2 (siRUNX2#1, siRUNX2#2) and negative mimic (NC) was verified 48 h later by utilizing the qRT-PCR assay.

To verify the cytotoxicity of the indicated anti-cancer drug at various concentrations, the cells were seeded in 96-well plates at a density of 5000 cells/well in 100 µl culture medium and incubated at 37°C, 5% CO₂ for 24 h. Then, the cells were incubated with anti-cancer drug epirubicin (Pfizer) for 24 h, and 10 µl/well CCK-8 solution (GK10001, GlpBio) was added for another 2 h. Finally, the absorbance value was read at 450 nm by using a microplate reader (Thermo).

Total RNAs were isolated by using TRIzol Reagent (CW0580S, CWBIO) according to the manufacturer’s instructions. Then, cDNAs were synthesized with the Reverse Transcriptase kit (E047-01B, Novoprotein) as we previously reported (12). After that, qRT-PCR analysis was conducted by using the SYBR Green Pro Taq HS premixed qPCR kit (AG11701-S, Accurate Biotechnology) and the ABI PRISM 7900 (Applied Biosystems). In details, an initial denaturation step of 95°C for 30s, followed by 40 cycles of 95°C for 5s and 60°C for 30s. The results were calculated by using 2^-ΔΔCt method. The primer sequences were shown in Supplementary Table S1 . GAPDH was used as internal reference.

Total protein was extracted with RIPA lysis buffer from the indicated cells, and the concentration was measured with the BCA Protein Assay Kit (P0012, Beyotime). The samples were subjected to 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and proteins were transferred to a PVDF membrane (IPVH00010, Millipore). After that, the membranes were incubated with primary antibody overnight at 4°C and secondary antibody for 1 h at room temperature, respectively. Finally, the protein bands were developed by using an enhanced chemiluminescence detection reagent (D045, Bridgen) and imaged by using an Image Lab (BIO-RAD). Primary antibodies against RUNX2 (12556, CST, dilution 1:1000), MMP1 (DF6325, Affinity, dilution 1:1000), E-CADHERIN (20874-1-AP, Proteintech, dilution 1:5000), VIMENTIN (10366-1-AP, Proteintech, dilution 1:5000), and β-ACTIN (AF7018, Affinity, dilution 1:5000) were used. Three independent experiments were conducted.

Cells preconditioned in DMEM basal medium (without 10% FBS) for 12 h were seeded in the upper chamber of the 24-well plate transwell insert (353097, Corning) at a density of 10⁵ cells/well, and 600 µl of DMEM complete medium (with 10% FBS) were added to the lower chamber. Then, the cells were fixed with 10% formalin (YULU) for 20 min after incubation for 24- 48 h, following by staining with 0.1% crystal violet (C8470, Solarbio) for 20 min. After wiping the non-migratory cells in the inner side of the upper chamber, images of migrated cells were captured from three randomly selected fields under the microscope, and the number of cells in each field was calculated.

The invasion assay was similar to the migration assay, where the upper chamber was coated with matrix gel (356234, Corning).

Breast cancer cells were seeded in 6-well plates at 1000 cells/well, and cultured for 14 days at 37°C, 5% CO₂. The cells were fixed with 10% formalin for 20 min, followed by staining with 0.1% crystal violet for 20 min. Typical colonies (foci > 100 µm) were counted. Three independent experiments were conducted.

To evaluate the malignant proliferative capacity of the indicated breast cancer cells, we conducted the soft agar colony formation assay. Briefly, the bottom layer was composed of 1.2% low melting point agarose (BS144, Biosharp) for preventing cell adhesion, and the upper layer was made up of 0.7% low melting point agarose containing 2000 cells/well. After 14 days, the morphology and count of colonies were recorded under the microscope.

5-week-old female BALB/c nude mice were applied for xenograft assay according to the Guide for the Care and Use of Laboratory Animals and were approved by the Experimental Animal Ethics Committee of Anhui Medical University (Approval No.: LLSC20210808). 5×10⁶ NC or shRUNX2-transduced MDA-MB-231-Re cells were subcutaneously injected into 5-week-old female BALB/c nude mice (N=4) in each group. Tumor long diameter (L) and short diameter (W) were measured every 3 days. The volume of the xenograft was calculated according to the formula: volume (mm³(=L×W²×0.5. Finally, the morphology and weight of the xenograft in nude mice in the indicated groups were recorded at day 33 of the xenograft assay. Besides, mice were sacrificed and tumors were resected for immunohistochemical staining.

For preparation of RNA-SEQ samples, total RNAs were isolated by using the TRIzol Reagent according to the manufacturer’s instructions. Then, the RNAs were sequenced by Personal Bio Company with the Illumina mRNA-Seq sample preparation kit and the Illumina Hiseq sequencing platform (Illumina, Inc.). The multifaceted bioinformatic analyses, including HeatMap, volcano plots, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were accomplished as we recently reported with several modifications (13).

Cells were crosslinked with 37% formaldehyde and lysed in SDS buffer. Then, the DNA was fragmented by using sonication. For RUNX2 binding sites on MMP1 promoter, ChIP assay was performed as we previously reported (14).

The eluted DNA fragments were turned to qRT-PCR assay, and the values generated from immunoprecipitated samples were normalized to those of the corresponding input samples. Besides, qRT-PCR products were stored at -20°C in a refrigerator until subsequent agarose gel electrophoresis (AGE). Primers sequences were shown in Supplementary Table S2 .

For the purpose of exploring the correlation of MMP1 in breast cancer, we analyzed the expression level of MMP1 according to the TCGA datasets. In details, we verified the MMP1 expression level in subtypes of breast cancer by the UALCAN database (http://ualcan.path.uab.edu). To assess the potential clinical significance of MMP1 in breast cancer, we took advantage of the Kaplan-Meier Plotter database (http://kmplot.com/analysis/index.php). The potential binding sites of RUNX2 in the promoter region of MMP1 were predicted via utilizing the JASPAR database (http://jaspar.genereg.net/).

Student’s t-test was used to analyze the differences between the two groups. All statistical analyses were performed with the GraphPad Prism (GraphPad Software) and SPSS 20.0 software as we described before (14). Only when P<0.05 was considered statistically significant differences. All data were shown as mean ± SD (n=3 independent experiments). *,P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001.

To generate the triple negative breast cancer resistant cell line MDA-MB-231-Re cells, we took advantage of the chronic-induced MDA-MB-231 cells for 8 months by step-up therapy (11). According to the CCK-8 values of the half of inhibition rate (IC₅₀), we found the drug resistance of MDA-MB-231-Re cells was higher than that of MDA-MB-231 parental cells (15.10 nmol/l vs. 2.035 nmol/l). In details, the Epirubicin drug resistance index of MDA-MB-231-Re cells was 7.42-fold over that of the MDA-MB-231 cells ( Figure 1A ). Meanwhile, we noticed the visible differences in cytomorphology between the MDA-MB-231-Re cells and MDA-MB-231 parental cells ( Figure 1B ). Distinguish from the epithelial-like shape in the MDA-MB-231 group, MDA-MB-231-Re cells showed typical mesenchymal-like morphology instead ( Figure 1B ). To investigate whether MDA-MB-231-Re cells could facilitate EMT, qRT-PCR and western-blotting analyses were performed to detect the expression pattern of EMT-associated biomarkers in the MDA-MB-231-Re cells and MDA-MB-231 parental cells. As show in Figures 1C, D , the expression of epithelial-related biomarker E-CADHERIN was reduced, while the mesenchymal-related biomarker VIMENTIN was up-regulated. Taken together, the triple negative breast cancer resistant MDA-MB-231-Re cells revealed enhanced drug-resistance against Epirubicin and facilitated EMT compared to the parental MDA-MB-231 cells.

To assess the influence of RUNX2 in the migration and invasion of triple negative breast cancer, we transduced lentivirus-mediated shRUNX2 into MDA-MB-231-Re cells (denoted as shRUNX2) to knockdown the endogenous expression of RUNX2. By conducting qRT-PCR and western-blotting analyses, we found the expression of RUNX2 both at the transcription and translation level was consistently impaired in the shRUNX2 group compared with the negative control group (denoted as NC) ( Figures 2A, B ). Strikingly, MDA-MB-231-Re cells with stable RUNX2 knockdown restored the epithelial-like cytomorphology, which was distinguish from that in the NC group ( Figure 2C ).

Considering the switching action in orchestrating the cytomorphology, we supposed that RUNX2 might play a pivotal role in the progression of epithelial-mesenchymal transition (EMT). Therefore, qRT-PCR and western-blotting analyses were performed to verify the expression pattern of EMT-associated biomarkers in the shRUNX2 group and NC group, respectively. Intuitively, we found that RUNX2 knockdown increased the expression of E-CADHERIN, while reduced the expression of VIMENTIN ( Figures 2D, E ). Consistently, we took advantage of the transwell assay to evaluate the influence of RUNX2 knockdown upon MDA-MB-231-Re cells. Notably, the migration and invasion ability of the MDA-MB-231-Re cells with RUNX2 knockdown was consistently impaired ( Figures 2F, G ). Collectively, RUNX2 knockdown was adequate to inhibit EMT and reduce the in vitro migration and invasion in the triple negative breast cancer resistant MDA-MB-231-Re cells.

Having clarified the inhibitory effect of RUNX2 knockdown upon cytomorphology as well as migration and invasion, we were next curious about the potential impact upon proliferation, drug resistance and anchorage-independent growth in MDA-MB-231-Re cells. With the aid of plate cloning assay, we found that fewer clones were formed in the shRUNX2 group compared to that in the NC group ( Figure 3A ). CCK-8 assay showed that the IC₅₀ of MDA-MB-231-Re cells with shRUNX2 transfection was much lower than that in the NC group (4.875 nmol/l vs. 17.74nmol/l) as well, which indicated the weakened drug resistance according to the epirubicin resistance index (RI=27.48%) ( Figure 3B ). In consistence, soft agar colony formation assay further revealed the reduced anchorage-independent growth capacity of shRUNX2 cells compared to NC cells ( Figure 3C ). Taken together, the proliferation, drug resistance, and anchorage-independent growth in triple negative breast cancer MDA-MB-231-Re cells were consistently impaired by RUNX2 knockdown.

To further verify the oncogenic potential of RUNX2 in vivo, the aforementioned NC cells and shRUNX2 cells were subcutaneously injected into BALB/c nude mice for xenograft tumor formation. Distinguish from those formed in the NC group, tumor sizes in the shRUNX2 group were significantly decreased at day 33 after transplantation ( Figures 4A, B ). Furthermore, as shown by the spatio-temporal growth curves, the tumor volume and tumor growth rate of shRUNX2 cells were less than that of the NC cells ( Figure 4C ). Similarly, the tumor weight in the shRUNX2 group was also sharply declined when compared with the NC group ( Figure 4D ). Overall, our data indicated the crucial role of RUNX2 in tumorigenesis during triple negative breast cancer.

Having verified the cellular characteristics of shRUNX2 cells and NC cells for triple negative breast cancer, we next turned to dissect the underlying molecular mechanisms. Therewith, we took advantage of the RNA-SEQ analysis to explore the potent similarities and differences at the transcriptomic level. As shown by the principal component analysis (PCA), the three independent duplicate sample (rep1, 2, 3) in the shRUNX2 group (in pink) or the NC group (in green) were clustered, respectively ( Figure 5A ). Hierarchical cluster analysis further showed the discrepancy in the profiling of differentially expressed genes (DEGs) between the indicated two groups ( Figure 5B ). Furthermore, the volcano and MA plots intuitively revealed the distribution of the DEGs including 896 up-regulated genes and 217 down-regulated genes (|log2Fold Change|>1, P <0.05) ( Figures 5C, D ).

Aiming to clarify the biological significance of the DEGs, we performed top gene ontology (topGO) analysis and observed the enrichment of certain datasets such as extracellular matrix (cell components, short for CC), protein binding (molecular function, short for MF), and regulation of cell population proliferation (biological processes, short for BP) ( Figure 6A ). Meanwhile, EMT-associated subsets were enriched by conducting GO enrichment analysis, including extracellular region, cell periphery, cell communication, regulation of cell population proliferation ( Figure 6B ).

Subsequently, we took advantage of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database to further dissect the underlying variations accompanied with the DEGs. Notably, a certain number of signaling pathways associated with proliferation and tumorigenesis were enrich based on the -log10 (P value) such as PI3K-Akt, MAPK and Jak-STAT signaling pathways ( Figure 6C ). In consistence, KEGG pathway enrichment assay based on FDR values (FDR<0.001) further confirmed the specific enrichment of pathways in cancer, Jak-STAT, MAPK and PI3K-Akt signaling pathways, together with cell adhesion molecules(CAMs) ( Figure 6D ).

Therewith, we detected the mRNA expression pattern of extracellular matrix-, cell adhesion-, cell invasion and metastasis-associated DEGs in triple negative breast cancer MDA-MB-231-Re cells with stable NC or shRUNX2 transfection. Of the indicated candidates (e.g., TGF-β, MDR1, NTN4, NOTCH1), we noticed that the expression of MMP1 with the most downregulation in shRUNX2 cells when compared with the NC group ( Figure 6E ). Collectively, our data indicated the biological significance of DEGs were mainly related with EMT process (e.g., migration, invasion) and resistance. MMP1 might play a critical role in mediating the tumor-promoting effect of RUNX2 in triple negative breast cancer cells.

For the purpose of further exploring the correlation of MMP1 in breast cancer, we detected the MMP1 expression level in 1098 breast cancer samples and 113 normal breast tissues according to the TCGA datasets. Notably, we found breast cancer tissues with higher level of MMP1 expression over that in normal breast tissue, which was confirmed by the UALCAN database (http://ualcan.path.uab.edu) ( Figures 7A, B ). Furthermore, compared with other types of breast cancer, we noticed that triple negative breast cancer with the highest level of MMP1 expression ( Figure 7C ).

Meanwhile, we turn to assess the potential clinical significance in breast cancer. Based on the Kaplan-Meier Plotter database (http://kmplot.com/analysis/index.php), we found that MMP1 expression showed positive correlation with the overall survival, relapse free survival and distant metastasis free survival of patients with breast cancer ( Figures 7D–F ). Taken together, our data revealed the positive correlation of MMP1expression with poor prognosis of breast cancer.

To further illustrate the molecular mechanism of MMP1 and RUNX2 in TNBC, we initially examined the expression pattern of RUNX2 and MMP1 in the SUM-149 and MDA-MB-231 triple negative breast cancer cell line and the MCF-7 non-triple negative breast cancer cell line. Consistently, we found that the mRNA expression levels of RUNX2 and MMP1 revealed similar expression pattern, and in particular, those in triple negative breast cancer cell line were significantly higher than in non-triple negative breast cancer cell line ( Figure 8A ). Furthermore, by conducting qRT-PCR and western-blotting assays, we found the expression levels of MMP1 and RUNX2 both at transcription and translation level were higher in MDA-MB-231-Re cells compared with those in MDA-MB-231 cells ( Figures 8B, C ). Then, with the aid of shRNA- and siRNA- mediated knockdown of RUNX2, we observed the sharp decrease of MMP1 expression in MDA-MB-231-Re and MDA-MB-231 cells compared to the NC groups ( Figure 8D ). Simultaneously, the expression of MMP1 was significantly upregulated in SUM-149 breast cancer cells with RUNX2 overexpression ( Figure 8E ). Besides, the expression levels of RUNX2 and MMP1 in tumor sections from both groups of xenograft tumor were detected by immunohistochemical staining, which showed that they were obviously lower in tumors of the shRUNX2 group compared with those of the NC group ( Figure 8F ).

Having preliminarily verified the potential regulatory relationship between RUNX2 and MMP1, we assumed whether RUNX2 functions via direct activation of MMP1 expression. Therefore, we predicted the potential binding sites of RUNX2 in the promoter region of MMP1 according to the JASPAR database (http://jaspar.genereg.Net/) ( Figures 9A, B ; Table 1 ). Of the seven candidates, the top 3 binding sites (site 1, 2, 3) with the highest score were selected, for further identification ( Figure 9C ). Then, with the aid of chromatin immunoprecipitation (ChIP) and qRT-PCR analyses together with agarose gel electrophoresis, we found only site 2 was specifically enriched by RUNX2 in MDA-MB-231-Re cells with high level of RUNX2 expression ( Figures 9D, E ). This was finally verified by agarose gel electrophoresis ( Figure 9E ). Overall, these data collectively indicated the direct modulation of MMP1 by RUNX2 via specific binding motif within the promoter region.