Differential mRNA expression of SIX1 and its relationship with prognosis

We obtained the standardized pan-cancer dataset, TCGA TARGET GTEx (PANCAN, N = 19,131, G = 60,499), from the UCSC (https://xenabrowser.net/) database. Specifically, we extracted the expression data of the ENSG00000126778 (SIX1) gene from various samples, including solid tissue normal, primary solid tumor, primary tumor, normal tissue, primary blood derived cancer—bone marrow, primary blood derived cancer—peripheral blood samples. To ensure consistency, we performed a log2(x + 0.001) transformation for each expression value. We further filtered out cancer species with less than 3 samples, resulting in a final set of 34 cancer species with their corresponding expression data. To investigate the expression of SIX1 in the four major subtypes of breast cancer (luminal A, luminal B, HER2-positive, and triple-negative breast cancer), we downloaded RNA-sequencing expression (level 3) profiles and corresponding clinical information for breast cancer from the TCGA dataset (https://portal.gdc.com). Additionally, we obtained the current-release (V8) GTEx datasets from the GTEx data portal website (https://www.gtexportal.org/home/datasets). To explore the association between SIX1 mRNA expression and prognosis, we utilized two databases: Kaplan–Meier plotter (http://kmplot.com/analysis/) [21] and UALCAN (http://ualcan.path.uab.edu) [22]. Firstly, we used Kaplan–Meier plotter to analyze the overall survival (OS) curves of patients with high and low SIX1 expression in the aforementioned four major subtypes of breast cancer. Patients were classified as having high or low SIX1 expression based on the median SIX1 expression. Subsequently, we employed the UALCAN database to investigate the differential mRNA expression of SIX1 in various tumor histologies of breast cancer, as well as individual cancer stages and nodal metastasis statuses.

The relationship between the expression of SIX1 and ploidy and immune infiltration

The TCGA Pan-Cancer (PANCAN) dataset, consisting of 10,535 cases and 60,499 genes, was downloaded from the UCSC database (https://xenabrowser.net/). This dataset provides a standardized and comprehensive collection of pan-cancer data. To investigate the relationship between SIX1 expression, ploidy, and immune infiltration, we specifically extracted the expression data of the ENSG00000126778 (SIX1) gene from multiple samples. To ensure the reliability of the data, we focused on primary blood-derived cancers, including peripheral blood and primary tumor samples. Ploidy data for each tumor were obtained from a previous study [23]. We then integrated the ploidy data with the gene expression data. To facilitate analysis, a log2(x + 0.001) transformation was applied to each value. To further streamline the analysis, we excluded cancer species with fewer than 3 samples within each specific cancer type. As a result, we obtained expression data for 37 different cancer species. Subsequently, we explored the relationship between SIX1 expression and two key indicators: stromal score (which reflects the presence of stroma in tumor tissue) and immune score (which denotes the infiltration of immune cells in tumor tissue) [24]. For this analysis, we utilized the SangerBox website (http://vip.sangerbox.com/home.html), a valuable online platform specifically designed for TCGA data analysis. Additionally, using the SangerBox website, we further investigated the relationship between SIX1 expression and various immune cell types, including B cells, M1 macrophages, M2 macrophages, monocytes, neutrophils, natural killer cells, CB4 + T cells, CB8 + T cells, Tregs, and dendritic cells.

The gene alteration landscape of SIX1

We utilized cBioportal (https://www.cbioportal.org) [25] and TCGA (https://portal.gdc.cancer.gov/) databases to investigate gene alterations related to SIX1 in patients with breast cancer. To obtain the gene alterations from cBioportal, we first selected the Breast Invasive Carcinoma (TCGA, PanCancer Atlas) and selected genomic profiles, including mutations, putative copy-number alterations from GISTIC, and mRNA expression z-scores relative to all samples (log RNA Seq V2 RSEM). We selected all available samples as our patient/case set, then queried SIX1 to obtain information on genetic alterations, including their proportion and type, within breast cancer patients and subtypes. We then downloaded RNA-sequencing expression (level 3) profiles, genetic mutation data, and corresponding clinical information for breast cancer from the TCGA. The mutation data were downloaded and visualized using the maftools package in R software. Genes with higher mutational frequency detected in breast cancer patients were displayed as a histogram.

The mRNA and protein landscape of SIX1

RNA-sequencing expression profiles and associated clinical information for breast cancer were downloaded from the TCGA database. Differential mRNA expression was studied using the limma package in R software. Patients were grouped based on the expression levels of SIX1, with the experimental group comprising patients in the top 25% and the control group comprising patients in the bottom 25%. A threshold of "Adjusted P < 0.05 and Fold Change > 1.5 or Fold Change < − 1.5" was defined to identify differentially expressed mRNAs, and genes with differential expression were selected for further analysis. To further analyze the differentially expressed genes, the STRING database (https://string-db.org/) [26] was utilized to construct a protein–protein interaction network associated with SIX1. To investigate potential targets' function, gene functional enrichment analysis was conducted using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis. The ClusterProfiler package (version: 3.18.0) in R software was used to analyze GO functions and KEGG pathway enrichments to improve the understanding of mRNA carcinogenesis. Box plots were created using the ggplot2 package in R, while heatmaps were generated using the pheatmap package.

The relationship of SIX1 with cancer stem cells

We downloaded the TCGA Pan-Cancer (PANCAN, N = 10,535, G = 60,499) dataset, which is a standardized pan-cancer dataset, from the UCSC (https://xenabrowser.net/) database. Specifically, we extracted the expression data of the ENSG00000126778 (SIX1) gene from various samples. To ensure data quality, we selected only primary blood derived cancer—peripheral blood and primary tumor samples. Subsequently, we utilized the OCLR algorithm, as proposed by Malta et al. [27], to calculate the RNAss stemness score based on mRNA features. We then integrated the tumor stemness scores and gene expression data of the samples. Finally, to ensure statistical power, we removed cancer types with fewer than three samples, ultimately obtaining expression data for 37 different cancers. In each type of tumor, we computed their Pearson correlation coefficients. In the present study, we investigated the association between SIX1 and stem cell markers in breast cancer patients using data obtained from the TCGA database. Pearson correlation analysis was performed to explore the relationship between these factors.

Cell culture

MCF7 derivative cell lines were described previously [28]. The 66cl4 cell line, generously provided by Prof. Kongming Wu from Tongji Medical College, Huazhong University of Science and Technology, was employed in this study. We generated MCF7-SIX1 cell line by overexpressing SIX1 in MCF7 cells for convenience in description. MCF7-NC was used as the negative control for comparison. To knock down Six1 expression in 66cl4 cells, we targeted two sites of Six1 and generated two knockdown cell lines, named 66cl4-shSix1 KD1 and 66cl4-shSix1 KD2, respectively. The negative control for knockdown experiments was 66cl4-SCR. All cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin at 37 °C in 5% CO2.

Western blot analysis

We electrophoresed equal amounts of lysates, ranging from 30 to 50 μg, onto polyvinylidene difluoride membranes. Subsequently, the membranes were blocked using TBST with 5% milk and probed with primary antibodies, specifically β-actin (1:5000; Abmart, M20011F), SIX1 (1:1000, Cell Signaling, D5S2S), SOX2 (1:1000, abcam, ab92494), Oct4 (1:1000, abcam, ab181557), c-Myc (1:1000, Abmart, TA0358S), EPCAM (1:1000, huabio, EM1111), ALDH1A1 (1:1000, abcam, ab52492), ITGB1 (1:1000, abcam, ab179471), p-STAT3 (1:1000, Abmart, T56566F), and STAT3 (1:1000, Abmart, T55292F), overnight at 4 °C. After washing thrice with TBST, the membranes were incubated for 1.5 h at room temperature with secondary antibodies, including goat anti-rabbit IgG-HRP (1:10,000, Proteintech, SA00001-2) and goat anti-mouse IgG-HRP (1:10,000, Proteintech, SA00001-1).

Quantitative real-time PCR (RT-qPCR)

RNA was prepared using the RNAeasy™ Animal RNA Isolation Kit with Spin Column (Beyotime, R0026). cDNA was reverse transcribed from 1 μg total RNA using the HiScript® III All-in-one RT SuperMix Perfect for qPCR (Vazyme, R333-01). RT-qPCR was performed with the Taq pro Universal SYBR qPCR Master Mix (Vazyme, Q712-02). The primer sequences for RT-qPCR were as follows: Six1 (mouse): (forward) 5′-CAAGAACGAGAGCGTGCTCAAGG-3′, (reverse) 5′-GGTGATTGTGAGGCGAGAACTGG-3′; SIX1 (human): (forward) 5′-CAAGAACGAGAGCGTACTCAAGGC-3′, (reverse) 5′-GGTGGTTGTGAGGCGAGAACTG-3′; Oct4 (mouse): (forward) 5′-CATTGAGAACCGTGTGAGGTGGAG-3′, (reverse) 5′-GCGATGTGAGTGATCTGCTGTAGG-3′;OCT4 (human): (forward) 5′-GTGGTCCGAGTGTGGTTCTGTAAC-3′, (reverse) 5′-CCCAGCAGCCTCAAAATCCTCTC-3′; Sox2 (mouse): (forward) 5′-CAGCATGTCCTACTCGCAGCAG-3′, (reverse) 5′-CTGGAGTGGGAGGAAGAGGTAACC-3′; SOX2 (human): (forward) 5′-CAGCATGTCCTACTCGCAGCAG-3′, (reverse) 5′-CTGGAGTGGGAGGAAGAGGTAACC-3′; Aldh1a1 (mouse): (forward) 5′-ATGGTTTAGCAGCAGGACTCTTCAC-3′, (reverse) 5′-CCAGACATCTTGAATCCACCGAAGG-3′; ALDH1A1 (human): (forward) 5′-ACGCCAGACTTACCTGTCCTACTC-3′, (reverse) 5′-TCTTGCCACTCACTGAATCATGCC-3′; CD44 (mouse): (forward) 5′-CTCAAGTGCGAACCAGGACAGTG-3′, (reverse) 5′-ATCAGAGCCAGTGCCAGGAGAG-3′; CD44 (human): (forward) 5′-TCTACAAGCACAATCCAGGCAACTC-3′, (reverse) 5′-ATGGGAGTCTTCTTTGGGTGTTTGG′; The β-actin primers (Sangon Biotech, B661302) and β-ACTIN primers (Sangon Biotech, B661102) were purchased from Sangon.

Mammosphere formation and self‑renewal capability assay

Cells were dissociated into single cells by 0.05% trypsin‐EDTA solution and plated into Corning ultralow attachment culture dish (Corning, 3471) at a density of 2 × 103 viable cells per milliliter in primary culture. Cells were grown in a serum‐free DMEM medium supplemented with B27 (Gibco, 17,504,044), 20 ng/mL Epidermal Growth Factor (PeproTech, 315-09/AF-100-15), 20 ng/mL Basic fibroblast growth factor (PeproTech, AF-450–33/100-18B), 2 μg/mL heparin (MCE, HY-17567). For MCF‐7 cells, mammospheres were kept in culture 7 days. Whereas 66cl4 mammospheres were kept in culture 5 days. To assess self-renewal capacity, mammospheres (diameter > 50 µm) were manually enumerated and representative images captured using an OLYMPUS IX71 microscope (Tokyo, Japan). Mammosphere-forming efficiency was calculated as follows: (number of mammospheres per well/number of cells seeded per well) × 100.

Flow cytometry analysis

To detect the stem cell subpopulations, the following antibodies were used: APC anti-mouse-CD24 (Biolegend, 101,813, 1:167 dilution), PE anti-mouse-CD49f (Biolegend, 313,612, 1:167 dilution), APC anti-human-CD24 (Biolegend, 311,117, 1:167 dilution), PE anti-human-CD44 (Biolegend, 103,007, 1:167 dilution). A total of 1 × 106 cells were incubated with antibodies in the dark at 4 °C for 30 min. Cells were washed and re-suspended in 500 µl of PBS and analyzed using a flow cytometer (Beckman Coulter, CytoFLEX).

Aldehyde dehydrogenase (ALDH) activity

Cells were first placed on ice then ALDH was detected by an ALDH test kit (Solarbio, BC0755) as indicated by the manufacturer. All ALDH activities were evaluated using a microplate reader (PerkinElmer, Ensight) at 340 nm by measuring the production of NAD + . Higher optical density (OD) values indicate stronger activity.

Cell proliferation assay

The influence of SIX1 on cancer cells viability were determined with Cell Counting Kit 8 (CCK-8) assay (Biosharp, BS350B) according to the manufacturer’s instructions. For CCK-8 assay, 1000 cells of various 66cl4 and 2000 cells of various MCF7 were plated in 96-well plates. Day 0 time point was measured 6 h post plating. Following 24-, 48-, 72- or 96-h of incubation, day 1–4 time points were analyzed. 10 μl of CCK-8 was added to each well and incubated at room temperature for 60 min and luminescence was measured by using a microplate reader (PerkinElmer, Ensight) at 450 nm.

Tumor-bearing model and imaging

Prior to the commencement of experiments, female BALB/c mice (6–8 weeks old) were provided with ad libitum access to food and water. All animal studies were meticulously reviewed and granted ethical clearance by the Laboratory Animal Welfare & Ethics Committee of Renmin Hospital at Wuhan University (Issue No. 20200702). It should be noted that all animal experiments strictly adhered to the guidelines articulated in the Guide for the Care and Use of Laboratory Animals developed by the Institute of Laboratory Animal Research. Firstly, a total of 18 mice were randomly assigned to three distinct groups. Specifically, cells harvested from 66cl4-SCR-luc, 66cl4-shSix1 KD1-luc, and 66cl4-shSix1 KD2-luc were resuspended in serum-free medium at a density of 1 × 106 cells per 100 μl. Subsequently, tumor cells were administered orthotopically into the fourth mammary fat pad of the mice. The tumor volume and luminescence signals were monitored on a weekly basis via an in vivo imaging system (PerkinElmer, IVIS Spectrum). To assess the impact of Six1 knockdown on the tumorigenic capacity of 66cl4 cells in vivo, we performed gradient dilution experiments by injecting varying numbers of cells for comparison. A total of 54 mice were randomly divided into three groups. Cells obtained from 66cl4-SCR-luc, 66cl4-shSix1 KD1-luc, and 66cl4-shSix1 KD2-luc were suspended in serum-free medium at a density range of 1*104 to 1*106 cells per 100 μl. Subsequently, the tumor cells were orthotopically administered into the fourth mammary fat pad of the mice. Tumor sizes and weights were monitored starting from 10 days post-injection.

In vivo expression of cancer stem cell markers in tumors

Tumors from each aforementioned group were fixed with 4% paraformaldehyde at room temperature for 24 h, followed by paraffin embedding and sectioning into 5 µm thick slices. Deparaffinization of paraffin sections was carried out using standard techniques, and sodium citrate antigen retrieval was performed. The sections were permeabilized for 15 min with 0.1–0.25% Triton X-100 and blocked for 30 min with 10% goat serum. Subsequently, different sections were incubated overnight at 4 °C with specific primary antibodies for Oct4 (1:1000, abcam, ab181557), Sox2 (1:100, abcam, ab92494), Aldh1a1 (1:100, abcam, ab52492), Epcam (1:20,000, abcam, ab213500), and Itgb1 (1:1000, abcam, ab179471). Following this, goat anti-rabbit Immunoglobulin G H&L (HRP) (1:100, abcam, ab205718) was incubated at room temperature for 1 h. Finally, all slides were stained with 300 nM DAPI for 5 min at room temperature. Fluorescence microscopy at × 40 magnification was used to observe images in each sample.

Statistical analysis

Statistical analyses were performed using GraphPad Prism 9.0 software (GraphPad Software). The Student's t-test was employed for two-group comparisons, while one-way ANOVA non-parametric was used for comparisons involving three or more groups. Proliferation data and tumor growth curves were subjected to two-way ANOVA analysis. Statistical significance was indicated by P values < 0.05. Detailed P-values can be found in the figures.