SIX2 Mediates Late-Stage Metastasis via Direct Regulation of SOX2 and Induction of a Cancer Stem Cell Program

The molecular mechanisms behind the metastatic progression of breast cancer remain poorly understood. Although advances in therapeutics are extending the lives of patients, metastatic spread remains the major cause of death. Understanding the mechanisms driving metastasis, particularly establishment and maintenance of metastasis, could lead to development of more effective therapeutics and increased survival.

The metastatic cascade is a multistep process involving escape from the primary tumor, dissemination, extravasation, survival/establishment, and colonization. The rate-limiting step of metastasis is outgrowth at the secondary site, as newly arriving tumor cells must survive a harsh new microenvironment to establish and maintain proliferation at secondary sites (1). Outgrowth of metastatic cells is complicated by the fact that these sites initially lack sufficient vasculature to provide nutrients needed by the cells (2). Because the microenvironment at the secondary site presents challenges for newly arriving tumor cells, only a small percentage of these cells will go on to develop overt metastases (3). Despite increased evidence around the importance of metastatic cell outgrowth at secondary sites, the vast majority of metastasis studies continue to focus on identifying molecular mechanisms that drive the earlier stages of metastasis, such as migration and invasion. However, recent studies have established that tumor cells exit the primary site at very early stages of the diseases (4). Thus, targeting early stages of metastasis may be less effective, because by the time the primary tumor is diagnosed, many cancer cells would have likely already left the primary site. Because the later stages of metastasis are less efficient, and often not completed at the time of diagnosis, targeting this aspect of the metastatic cascade may be more effective (1).

Recently, advancements in imaging, mouse models, and genomic and single cell genomics analyses have allowed for more accurate modeling of the later stages of metastasis. This technological progress has led to the identification of several proteins that specifically alter metastatic colonization. For example, tumor cells can indirectly activate osteoprotegerin, allowing for bone degradation and release of growth factors into the bone microenvironment, which provides signals to metastatic cells to proliferate (5). Additionally, tumor cell–derived tenascin C promotes the survival and outgrowth of lung metastases through regulation of stem cell signaling components musashi homolog 1 and leucine-rich repeat–containing G protein–coupled receptor 5 (6). However, even with these promising discoveries, metastatic colonization remains the major cause of mortality in patients with breast cancer. In part, this continued mortality can be attributed to a dearth of targeted therapies inhibiting metastatic colonization and/or the outgrowth of metastatic clones, which are often resistant to the chemotherapy regimens to which the primary tumor responded (1).

SIX2 is a homeodomain-containing transcription factor that plays a key developmental role in the kidney. KO of six2 in mice demonstrates a critical role for the gene in maintenance of the mesenchymal phenotype, self-renewal, and survival of progenitor cells that will eventually give rise to all epithelial cell types of the developing nephron (7). Not surprisingly then, Six2 loss in mice leads to postnatal lethality due to defects caused by premature differentiation of cells in the kidney and kidney hypoplasia (7). Importantly, mutations in SIX2 have been found in humans with congenital abnormalities of the kidney, such as renal hypodysplasia (7), underscoring an important role for this gene in human kidney development.

Tumor cells often hijack developmental processes to promote metastasis (8). Recent studies have thus begun to examine the role of SIX2 in tumor progression. Overexpression (OE) of SIX2, as compared with normal tissue, has been observed in esophageal and lung carcinomas (9). SIX2 also plays a critical role in tumor cell invasion and drug resistance in colorectal cancer (10), and regulates proliferation and epithelial-to-mesenchymal transition (EMT) in hepatocellular cancer (11). Most recently, whole exome sequencing of Wilms' tumor patient samples identified that almost a fifth (18.1%) of patients with the most aggressive blastemal subtype of Wilms tumor possess a gain-of-function mutation (Q177R) in the homeodomain of SIX1 or SIX2 that correlates with high proliferation (12). This mutation has been shown to enhance the SIX1-mediated Warburg effect, suggesting that a similar function could be attributed to SIX2 (13). Importantly, SIX2 marks the cancer stem cell (CSC) population in Wilms tumor (14), suggesting conservation of its developmental function in cancer. Collectively, these data support a critical role for SIX2 across multiple different cancer types and underscore the need to better understand the molecular mechanisms by which SIX2 promotes tumor progression and metastasis.

Previous studies from our lab have demonstrated that Six2 is critical for late-stage metastasis in triple-negative mouse mammary carcinoma models (9). Knockdown (KD) of Six2 in the triple-negative, metastatic 66cl4 mouse mammary carcinoma line results in decreased metastatic burden when injected orthotopically (9). Interestingly, Six2 KD does not affect primary tumor growth or tumor-associated lymphangiogenesis, contributors to early-stage metastasis that are known to be regulated by a related Six family member, Six1 (8). Conversely, Six2 OE in the 4T07 mammary carcinoma cell line results in increased lung metastases after tail vein injection in BALB/c mice when compared with control cells (9). This result is of particular interest, as the 4T07 cells, which are syngeneic to the 66cl4 cells, can reach secondary sites but are unable to efficiently colonize them (15). These data demonstrate that Six2 promotes the ability of cancer cells to colonize secondary sites. Further analysis demonstrated that Six2 OE leads to a decrease in E-cadherin (Cdh1) expression, suggesting that Six2 may mediate at least a partial EMT (9). Importantly, restoration of E-cadherin in 4T07 cells downstream of Six2 OE inhibited metastasis induced by Six2 (9). However, KD of Six2 inhibited metastasis of 66cl4 cells without restoration of E-cadherin expression, likely due to epigenetic silencing of Cdh1 (9). Taken together, these data suggest that Cdh1 repression is necessary downstream of Six2 to mediate metastasis, but that additional genes regulated by Six2 are also required.

Herein, we demonstrate that SIX2 plays a critical role in metastatic colonization by promoting stemness-associated properties in triple-negative breast cancer (TNBC) and TN mouse mammary carcinoma cells. We show that SIX2 is more highly expressed in patients with TNBC as compared with other subtypes, and that it regulates metastasis not only of mouse mammary carcinoma lines, but also in a human TNBC model. Using three different models of TNBC, we demonstrate that SIX2 enhances cancer stem-cell associated phenotypes and regulates a genetic stem cell program. We identify Six2 as a direct transcriptional regulator of a master pluripotency factor, Sox2, further resulting in upregulation of a second master pluripotency factor, Nanog. Importantly, we show that Sox2 is the critical mediator downstream of Six2 in promoting stemness-associated phenotypes in vitro and late-stage metastasis in vivo. Finally, we demonstrate for the first time that SIX2 and SOX2 expression positively correlate in human breast cancer, including TNBC, and that a SIX2-mediated gene signature is associated with significantly shortened distant metastasis-free survival and relapse and recurrence-free survival. Collectively, our data suggest that a novel SIX2/SOX2 axis may promote stem cell characteristics in newly arriving breast cancer cells at secondary sites, providing insight into the regulation of metastatic colonization and outgrowth.

The 4T07 and 66cl4 mouse mammary carcinoma cells were generously provided by Dr. Fred Miller (15). Stable Six2 OE and KD cell lines were generated and grown in media as described previously (9). MDA-MB-231 cells have been cultured long term in the Ford laboratory, and were grown in MEM/EBSS, FBS, HEPES, nonessential amino acids, insulin, pen/strep, and l-glutamine. Authentication of the human cell line MDA-MB-231 was performed using STR analysis by the University of Colorado Cancer Center Tissue Culture Shared Resource in March of 2017. All cells were regularly monitored for Mycoplasma contamination every 3 months and only Mycoplasma-free lines were used for studies. All cell lines were cultured at 37°C in HEPA-filtered humidified air in 5% CO₂. Stable SIX2 KD in the MDA-MB-231 cells was performed by transduction of two different lentiviral shRNAs, or with a nontargeting (NT) shRNA. After transduction, cells were selected for 1 week using 3 μg/mL of puromycin to obtain stable KD of SIX2. SIX2 levels in KD and OE cells were examined using RT-qPCR and Western blot analysis. Transient KD of Sox2 and SOX2 was performed using ON-TARGETplus SMARTpool siRNA constructs (Supplementary Table S1). A NT siRNA pool was used as a negative control. KD levels were assessed 72 hours posttransfection via RT-qPCR.

4T07 control (CTL) and Six2 cells were plated at 3 × 10⁵ cells/mL in 10 cm plates and allowed to attach and grow for 48 hours. Cells were collected and lysed by removing media, rinsing with PBS and adding 1 mL of TriReagent directly to the plate. After vortexing to ensure homogeneity, total RNA was isolated using the Zymo Research RNA Isolation Kit (Light Labs, R2053). RNA samples were evaluated for purity and concentration using a spectrophotometer and submitted to the Genomics and Microarray Core Facility at the University of Colorado Denver Anschutz Medical Campus. Sequencing sample preparation was performed using standard Illumina HiSeq protocols and reagents. Directional mRNA Seq was performed using Illumina HiSeq 4000 HT Mode 1 × 150 cycles.

The read quality and mapping were done as described previously (16). Significantly differentially expressed genes were determined after using DESeq2 in R with batch correction, and a significance cutoff of adjusted P value <0.1. MA plots, PCA plots and heatmaps were generated using Matplolib from Python plotting library. To find overlap between significantly differentially expressed genes and stemness-associated gene signatures, DESeq genes were input into the StemChecker database (http://stemchecker.sysbiolab.eu). Significantly enriched gene signatures were determined using a Hypergeometric test and adjusted by Bonferroni correction. Heatmaps were generated by converting reads per kilobase per million (RPKM) values to z scores and using ComplexHeatmap and Circularize library R packages to visualize z-score normalized gene expression. Genes were grouped by hierarchical clustering of z-score normalized gene expression. The RNA-sequencing (RNA-seq) data were deposited into the GEO repository with the accession number GSE123489.

Gene-set enrichment analysis (GSEA) was performed by utilizing the GSEA Preranked Module of the GenePattern Webserver (17). All significantly expressed genes were ranked by adjusted log₂-fold change values (Six2 OE vs. CTL) and tested against the C2: Curated gene sets and C5: Gene Ontology gene sets to test for enrichment.

Whole cell and nuclear lysates were isolated as done previously (9). Protein concentration was determined using the Lowry Protein assay (Bio-Rad). After gel electrophoresis of protein in 10% polyacrylamide gels, proteins were transferred to a polyvinylidene difluoride membrane at 100 v for 1 hour 30 minutes, blocked in 5% Milk-TBST and incubated at 4°C in primary antibodies overnight (see Supplementary Table S1, for list of primary antibodies). Blots were imaged using either film radiography or the OdysseyFc imaging system after horseradish peroxidase substrate chemiluminescence incubation and detection.

RNA was isolated using the Zymo Research RNA Isolation Kit (Light Labs, R2053). cDNA was made using the QuantaraBio cDNA Synthesis Kit (VWR Scientific, 95047-100). PCR was performed using the SYBR Green Master Mix per manufacturer's instructions and protocols (Qiagen; 330500). Primers used for this study are contained in Supplementary Table S1 (IDT). Gene expression was normalized to peptidylpropyl isomerase B (cyclophilin B) (Ppib/ppib) mRNA expression and fold-change is relative to the CTL, nonsilencing, or NT cells. Significance was determined by t test or one-way ANOVA on means ± SEM of triplicate samples for a representative experiment (n = 3).

Tumorsphere assays were performed as previously described (18). Briefly, cells were trypsinized and plated at either 300, 1,000, or 20,000 cells/well in six-well, ultra-low attachment plates (Corning) in 2 mL of serum-free DMEM/F12 media (Hyclone), supplemented with 20 ng/mL bFGF (BD Biosciences), 20 ng/mL EGF (BD Biosciences), 4 μg/mL heparin (Sigma), penicillin–streptomycin (Hyclone), and B27 (Gibco) for all cell lines tested. F500 μL of media was added every 3 to 4 days. On day 7, all cells/spheres were collected, washed with PBS, digested using 0.05% trypsin in Hank's Balanced Salt Solution with 0.2 g/L EDTA (Fisher), and single cells were plated as above to perform secondary tumorsphere assays. After an additional 7 days, tumorspheres were imaged at ×4 magnification and counted by the Incucyte Zoom Live-Cell Analysis System.

Flow cytometry analysis was done as previously described (18). Briefly, cells were incubated for 30 minutes in either CD24-biotin/CD49f-APC or Epcam-FITC/CD44-APC, and in the case of CD24 were additionally incubated in Strep-FITC or Strep-V450 secondary for 30 minutes (Supplementary Table S1). Positive control ultracomp beads were also stained and used as a positive control for consistent gating. Flow cytometry was performed by submitting samples to the Flow Cytometry Core Facility at the University of Colorado Anschutz Medical campus or using a C6 Flow Cytometer. Single-stained compensation controls (UltraComp Beads; eBioscience) were stained with individual antibodies (EpCAM-FITC, CD44-APC, CD24-biotin-strep-FITC, CD24-biotin-strep-V450, or CD49f-APC) to define gating for each cell-surface marker. Gates were set at the lowest edge of fluorescence for the positive stain populations relative to the negative stain populations. Additionally, negative control (unstained) cell populations were checked to make sure they overlapped <1% with the single or double positive quadrants after gating with single-stained control beads. Flow cytometry analysis was performed using the FlowJo Analysis Software.

Tail vein injections were performed in 6- to 8-week-old female immunodeficient NOD scid gamma (NSG; for MDA-MB-231) and BALB/c mice (for 4T07; The Jackson Laboratory). For the MDA-MB-231 cells, 4 × 10⁵ cells in 100 μL of Hank's Balanced Salt Solution were injected. For the 4T07 cells, 1 × 10⁵ cells in 100 μL of DMEM medium were injected. Both the MDA-MB-231 and 4T07 cells were luciferase tagged to enable imaging of metastases in vivo. IVIS imaging was used to image the mice as described previously (9). Quantification of luciferase signal was determined by Living Image Software. Limiting dilution assays were performed by injecting 5 × 10⁵, 2.5 × 10⁵, or 5 × 10⁴ luciferase-tagged MDA-MB-231 cells in 100 μL of a 1:1 Matrigel/PBS (Corning, 354234) underneath the nipple of fourth mammary fat pad of nude mice. Tumor formation was monitored weekly via IVIS imaging. TIC capacity was determined using the ELDA software website application (http://bioinf.wehi.edu.au/software/elda/) with a 95% confidence interval setting. All animal work was approved and performed in compliance with the Institutional Animal Care and Use Protocol at the University of Colorado Anschutz Medical Campus.

Cells were plated in 12-well plates to achieve equal confluence (4T07: 4 × 10⁴ cells/mL, 66cl4: 7.5 × 10⁴ cells/mL, and MDA-MB-231: 1 × 10⁵ cells/mL) and then cotransfected with the indicated firefly-luciferase constructs (pGF-mSox2-ONL and pGF-mSox2-T27G-ONL) and Renilla luciferase constructs (pCNL-N1_3xNLS) using lipofectamine 2000 according to manufacturer instructions (Thermo Fisher Scientific). After 72 hours, cells were collected and analyzed for luminescence using the Dual-luciferase Reporter Assay System (Promega). To assess relative SRR2 activity, luciferase luminescence was normalized to Renilla luciferase activity.

Cells were plated in 15 cm plates to similar confluency and then fixed with 1% formaldehyde to crosslink DNA and proteins. After terminating the crosslinking with glycine, cells were collected with ice-cold PBS containing protease inhibitors and pelleted by centrifugation. Cell pellets were then processed to obtain chromatin immunoprecipitation (ChIP) DNA by using the Zymo Research ChIP Kit (Zymo Research, D5210) according to manufacturer instructions. After obtaining ChIP DNA, qPCR was performed by diluting the samples 1:5 with nuclease-free water and using primers to amplify regions of SRR2 (Supplementary Table S1). Ct (quantification) values from qPCR were normalized to input values for each condition. Results are presented as relative enrichment by subtracting the %input for IgG negative control from %Input for Six2/SIX2 IP to adjust for background signal.

Publicly available datasets were analyzed using cBioportal, Oncomine, SurvExpress, or KMplotter. For survival data analyzed through KMPlotter, we used the 206510_at probe for SIX2 and the 228038_at probe for SOX2 (mean expression) to analyze the relationship between SIX2/SOX2 gene expression and distant metastasis-free survival and relapse-free survival. After removal of biased array data (patient and clinical heterogeneity, different outcome measures, and size effects) and redundant samples, results were split by the “auto select best cut off” and graphed to show the approximate low versus high expression, FDR of <20% and stated P values.

The experiments performed were analyzed using GraphPad Prism 6. Error bars were generated by SEM calculations. For experiments with two conditions, an unpaired one-tailed Student T test was performed. For experiments with three or more conditions, one-way ANOVA followed by a Bonferroni comparison was used. For the animal experiments measuring luciferase signal over time, fitting of a mixed effects model for testing the main effect of KD/OE and KD/OE–time interactions were performed. Given the statistical significance of the overall mode (P < 0.05), we proceeded to carry out the multiple comparisons (with Bonferroni correction) to determine significant differences in the increase of luciferase signal over time between the KD/OE groups and their corresponding CTL groups. Asterisks represent the significance of difference from control group: *, P < 0.05; **, P < 0.01; ***, P < 0.001, and ****, P < 0.0001.

Our previous studies demonstrated that SIX2 is more highly expressed in human breast carcinomas when compared with normal breast tissue, and that its expression is also increased in invasive ductal carcinomas when compared with ductal carcinoma in situ (9). Building on these results, analysis of SIX2 mRNA expression from multiple patient datasets, including The Cancer Genome Atlas (TCGA) and Tabchy datasets (19), reveals that it is more highly expressed in TNBC compared with all other subtypes of breast cancer (Fig. 1A). Patients with TNBC present with higher tumor grade at diagnosis, have the highest rate of recurrence, and possess the poorest prognosis overall compared with the other subtypes, making it one of the more aggressive breast cancer subtypes (20). Given that we previously established a role for Six2 in late-stage metastasis of murine triple-negative mouse models (9) and its high expression in human TNBC, we next sought to understand the molecular mechanism through which SIX2 regulates TNBC metastasis. However, as Six2 has previously only been shown to regulate murine triple-negative mammary carcinoma metastasis, we first examined whether SIX2 can also increase metastasis of human TNBC. To this end, we knocked down SIX2 in luciferase tagged triple-negative MDA-MB-231 cells (Fig. 1B). Tail-vein injections of NT control and SIX2 KD MDA-MB-231 cells into immunodeficient NSG mice demonstrate that mice injected with SIX2 KD cells have significantly reduced metastatic burden when compared with mice injected with NT control cells (Fig. 1C and D). These data, in conjunction with previous data in mouse mammary triple-negative models (9), demonstrate a critical role for SIX2 in the later stages of metastasis in human TNBC.

Because Six2 plays a critical role in the maintenance and self-renewal of nephron progenitor cells during kidney development, and because it promotes metastasis of a cell line (4T07) that is capable of reaching secondary sites, but cannot efficiently colonize those sites (15), we hypothesized that SIX2 mediates late-stage metastasis, particularly the survival and establishment of tumors at secondary sites, partly due to its role in stem/progenitor-like phenotypes. To test this hypothesis, we first performed flow cytometry to determine if there are differences in the stemness-associated population in response to the presence of SIX2. Because MB-MDA-231 cells express little to no CD24, we measured cell surface expression of EPCAM and CD44 (21), which is known to more accurately enrich for cells with tumor-initiation and metastatic capacity (22). Interestingly, SIX2 KD results in a decrease in the EPCAM⁺/CD44⁺population in MDA-MB-231 cells (Fig. 2A). We additionally examined whether Six2 alters the CSC population within triple-negative mouse mammary carcinoma cells where we had previously found it to mediate metastasis (9). In this case, flow cytometry was performed to examine the CD24⁺/CD49f⁺ population, which enriches for cells with the ability to repopulate the mouse mammary gland (23), in 4T07 cells overexpressing Six2 (4T07-Six2), and in 66cl4 cells with KD of Six2, as compared with their control counterparts. Our data show that Six2 loss decreases, whereas Six2 OE increases, the percentage of CD24⁺/CD49f⁺ cells within the population (Fig. 2B and C).

To assess the role of SIX2 in stem cell–associated phenotypes using a more functional assay, we next performed tumorsphere formation experiments in all three models. The ability of cells to form secondary tumorspheres after dissociation and dilution to single cells is a measure of self-renewal capability, and thus serves as a surrogate CSC assay (24). Similar to what was observed in our flow cytometry analysis, in both the MDA-MB-231 and 66cl4 cells, SIX2/Six2 KD significantly decreases the ability of cells to form secondary tumorspheres (Fig. 2D and E). In line with these findings, Six2 OE in the 4T07 cells significantly increases secondary tumorsphere formation (Fig. 2F).

To further test whether SIX2 regulates cancer stem cell frequency and functionality, we performed limiting dilution assays in vivo using the MDA-MB-231 cells with SIX2 KD. We found that compared with NT control cells, SIX2 KD cells have decreased tumor-initiating cell (TIC) capacity, as demonstrated by a significant decrease in TIC frequency from 1:146,096 for NT control cells to 1:813,444 and 1:357,622 for SIX2 KD1 and KD2 cells, respectively (Fig. 2G). Taken together, these data suggest that SIX2 may play a similar role in mediating stemness characteristics during tumor progression as it does during development (7), and that this role may allow it to mediate late-stage metastasis by enhancing the ability of cells to self-renew and survive at secondary sites.

Because SIX2 is a transcription factor, we hypothesized that it directly transcriptionally regulates downstream effectors to promote CSC-like phenotypes and metastatic burden. We have previously shown that 4T07 cells express very little Six2 endogenously and that the OE levels are similar to other murine mammary triple-negative cell lines with endogenously high Six2 expression, such as the 66cl4 cells (9). Additionally, although we observe a range of SIX2 expression in human TNBC cell lines, we found that Six2 OE in 4T07 cells is comparable to, and does not artificially exceed, the levels of SIX2 observed in multiple TNBC cell lines (Supplementary Fig. S1). Therefore, we chose this system to perform RNA-seq analysis comparing 4T07-CTL and 4T07-Six2 cells. Six2 OE leads to a robust transcriptional response as evidenced by the significant increased expression of 1,997 genes and decreased expression of 2,002 genes (FDR q value <0.1 and log₂-fold change ±1.0; Supplementary Fig. S2A). Additionally, principal component analysis (PCA) based on the top 500 most-variable genes from each sample revealed a stark difference between 4T07-CTL and 4T07-Six2 cells, indicating that Six2 OE in 4T07 cells drives distinct changes in the transcriptional landscape (Supplementary Fig. S2B). Given the well-established role of Six2 as a critical transcriptional regulator of progenitor and stem cell populations (7), we used the web-server StemChecker (25) to examine whether differentially expressed genes identified from the RNA-seq analysis overlap with previously established and curated stemness-associated gene signatures. We first queried the 1,997 significantly upregulated and 2,002 significantly downregulated genes in response to Six2 OE (adjusted P-value <0.1 and log₂-fold change ±1.0) from the RNA-seq analysis to determine the overlap with gene signatures from StemChecker that represent specific normal stem cell types. We found that 77% of the Six2-upregulated genes from the RNA-seq analysis overlapped with genes from StemChecker signatures representing various stem cell types, whereas the Six2-downregulated genes contained a 72% overlap (Fig. 3A). Interestingly, gene signatures representing embryonic stem cells (ESC) were among the most highly enriched and represented signatures in the analysis when examining all significantly differentially expressed genes together (Supplementary Fig. S2C).

We then determined whether Six2/SIX2 KD or Six2 OE (Fig. 3B shows expression levels of Six2 in all three cell lines as determined in the same RNA samples used in panels C and D) could affect the expression of genes from the ESC-associated signatures. For instance, when examining a dataset that used a shRNA functional screen to identify crucial factors for ESC proliferation and pluripotency (26), we found that the genes contained in this signature were dramatically differentially expressed between the 4T07-CTL and 4T07-Six2 cells (Fig. 3C, ESC_Wang dataset; ref. 26). The same was true for a signature containing genes that are differentially expressed between mouse ESCs and adult tissues (Fig. 3D, SC Wong dataset; ref. 27), as well as signatures associated with maintaining pluripotency (Supplementary Fig. S2D; refs. 28, 29). To validate that SIX2 mediates gene expression alterations that are linked to stem-associated phenotypes, we performed RT-qPCR for stem cell genes identified in the different signatures on independent sets of RNA from all three cell line models. These genes were chosen based on their established roles in normal embryonic/cancer stem cell function, tumor progression, and/or metastatic colonization (26, 30–34), and were among the top quartile of upregulated genes after DESeq analysis. This analysis confirmed that SIX2/Six2 KD in both MDA-MB-231 and 66cl4 cells leads to a corresponding decrease in the mRNA expression of numerous genes associated with ESCs (Fig. 3C and D), whereas Six2 OE leads to an increase in many of the same genes (Fig. 3C and D).

To further confirm that SIX2 plays a role in a genetic stem cell program, we also performed GSEA after converting mouse gene symbols to their human equivalents (35). Notably, GSEA analysis using the C2-Curated gene sets revealed that Six2 OE results in negative enrichment of pathways associated with E-cadherin function (ONDER CDH1 TARGETS 2 DN) and cell focal adhesion (KEGG FOCAL ADHESION; FDR q-value <0.25 and normalized P-value <0.05; Fig. 3E and F; Supplementary Table S2). These data support previous studies from our lab demonstrating that Six2 OE in 4T07 cells results in increased metastasis in part through downregulation of E-cadherin and increased anchorage-independent growth (9). Interestingly, we also observed enrichment of pathways associated with stemness. For instance, pathways that are upregulated to promote stemness (BENPORATH PROLIFERATION, WONG EMBRYONIC STEM CELL CORE, MUELLER PLURINET) were among the most highly positively significantly enriched pathways across all pathways analyzed in the C2-Curated gene sets (Fig. 3E and F; Supplementary Table S2). Additionally, we identified significantly negatively enriched pathways that contain genes that are downregulated in ESCs (BENPORATH ES WITH H3K27ME3; Fig. 3E and F; Supplementary Table S2). When performing GSEA analysis using the C5-Gene Ontology gene sets, Six2 OE led to negative enrichment of pathways involved in adhesion and differentiation (Supplementary Fig. S2E and S2F; Supplementary Table S2). Together, these data demonstrate that SIX2 expression broadly enhances a stem-cell gene program in both triple-negative breast and mammary cancer cells.

Based on the large number of stem-related genes altered by SIX2, we asked whether SIX2 promotes the expression of master regulators of stemness. To address this question, we first examined the mRNA and protein expression of well-established master regulators of stemness and pluripotency. Intriguingly, KD of SIX2/Six2 in the MDA-MB-231 and 66cl4 cells resulted in corresponding decreases in both SOX2 and NANOG mRNA as measured by RT-qPCR (Fig. 4A and B). In contrast, Six2 OE in the 4T07 cells resulted in upregulation of mRNA expression of both pluripotency factors (Fig. 4C). Additionally, KD of SIX2 leads to a decrease in both SOX2 and NANOG protein (Fig. 4D and E), whereas OE of Six2 leads to corresponding increases in Sox2 and Nanog levels (Fig. 4F). These data suggest that SIX2 may promote stemness-associated phenotypes through regulation of SOX2 and NANOG.

Previous studies have shown that in the context of breast cancer, SOX2 regulates cancer stem cell properties upstream of NANOG (36). To determine whether Sox2 is upstream of Nanog in Six2 OE cells, we established stable KD of Sox2 in the 4T07-Six2 cells (Fig. 5A). Indeed, we observed that Sox2 KD in 4T07-Six2 cells resulted in a corresponding decrease in Nanog mRNA and protein levels (Fig. 5B), suggesting that Six2 may act more directly on Sox2 to mediate alterations in Nanog levels.

Sox2 transcriptional regulation is, in part, under the control of two enhancer regions, Srr1 and Srr2. The Srr2 domain, located ∼3 kb downstream of the Sox2 gene, functions as a transcriptional enhancer of Sox2 expression in ESCs (37). Intriguingly, further analysis of the sequence of the Srr2 enhancer revealed the presence of a Six2 binding region (38). Therefore, we performed luciferase reporter assays containing the Srr2 regulatory region (39) to determine whether Six2 may regulate Sox2 through the Srr2 domain (Fig. 5C shows the reporter construct used in assays). After transfection of a Srr2 luciferase construct (pGF-Sox2-ONL, 5C; ref. 39), along with a Renilla luciferase plasmid (pCNL-N1_3xNLS) to normalize for transfection efficiency, we found that Six2 KD decreases, whereas Six2 OE increases, activation of the Srr2 region (Fig. 5D). To determine if Six2 directly regulates the Srr2 enhancer via this binding region, we introduced a single point mutation (T27G) within the putative Six2 binding region of the Srr2 enhancer based on a Six2 binding site that was previously identified and characterized in the Gdnf gene (Fig. 5E; ref. 38). Interestingly, although introduction of the wild-type Srr2 region resulted in a similar increase in Srr2 enhancer activity in 4T07-Six2 cells compared with 4T07-CTL cells, introduction of the T27G mutation into the Srr2 enhancer-luciferase construct inhibited the activation of luciferase activity (Fig. 5F). We next performed ChIP-qPCR and found that Six2 binding at the Srr2 enhancer region is enriched in 4T07-Six2 cells compared with control cells (Fig. 5G). In contrast, the enrichment of Six2 observed in the Sox2 regulatory regions is lost when performing a pulldown at a region approximately 6 kb downstream of the Srr2 region that contains no known Six2 binding sites (Fig. 5G). Together, our results suggest that Six2 plays a critical role in altering CSC characteristics through direct transcriptional activation of Sox2, resulting in upregulation of Nanog and induction of a stem cell-like program.

We then asked whether Sox2 is required downstream of Six2 to mediate CSC characteristics and late-stage metastasis. To this end, we performed flow cytometry to measure the mouse mammary stem cell enriched population in response to Sox2 KD downstream of Six2 in the 4T07 cells. We found that Sox2 KD in 4T07-Six2 cells leads to a significant decrease in the CD24⁺/CD49f⁺ enriched CSC population when compared with control cells (4T07-Six2 shNT; Fig. 6A). Similarly, secondary tumorsphere formation was also decreased in response to Sox2 KD in 4T07-Six2 cells (Fig. 6B). Additionally, to better understand why we did not observe a complete reversal in the in vitro stemness-associated phenotypes after Sox2 KD, we performed ChIP-qPCR in 4T07-CTL and -Six2 OE cells to examine whether Six2 was directly bound to, and therefore could transcriptionally regulate, the promoter regions of the putative targets identified in the StemChecker analysis. We found that Six2 was bound to a subset of these targets, suggesting that Six2 regulates multiple genes that a play a critical role in stemness (Supplementary Fig. S3). These data suggest that in order to completely reverse the in vitro stemeness-associated phenotypes regulated by Six2, multiple downstream targets, in addition to Sox2, may need to be inhibited.

To examine whether Sox2 is a critical mediator of Six2 in promoting late-stage metastasis, we performed tail-vein injections using 4T07-Six2/Sox2 KD cells (as compared with corresponding controls). As had been previously observed (9), Six2 OE in 4T07 cells enables metastatic outgrowth in vivo (Fig. 6C, top). However, KD of Sox2 in 4T07 cells dramatically decreases Six2-induced metastatic outgrowth (Fig. 6C, bottom and D for quantitation), suggesting that in vivo Sox2 is the key target downstream of Six2 in promoting late-stage metastasis. Thus, our data demonstrate that Six2 regulates metastatic outgrowth in a Sox2-dependent manner.

On the basis of our experimental findings that Sox2 KD decreases Six2-mediated metastatic outgrowth, we next examined whether this relationship may be relevant to human breast cancer. We first performed RT-qPCR on multiple human TNBC cell lines to determine if the SIX2–SOX2 axis is present in additional human TNBC cells. Indeed, we find that in both BT549 and SUM159 cells (TNBC cell lines with high endogenous SIX2 expression; ref. 9), transient KD of SIX2 leads to a corresponding decrease in SOX2 (Supplementary Fig. S4). These data suggest that similar to what we observe in the MDA-MB-231 cell line model, SIX2 also regulates SOX2 in other human TNBC models. We next performed correlation analysis of SIX2 and SOX2 mRNA expression obtained from cbioportal, Oncomine, and R2 Cancer Databases and demonstrate that SIX2 and SOX2 positively correlate in breast cancer datasets encompassing all breast cancer subtypes (Supplementary Fig. S5A); a relationship that is also observed when specifically examining TNBC datasets (Fig. 7A). These data suggest that the relationship between SIX2 and SOX2 in our cell line models is recapitulated in human breast cancer.

We then generated a Six2-gene expression signature compiled of the top 15 significantly upregulated genes in response to Six2 OE after converting mouse gene symbols to their human equivalent (Supplementary Fig. S5B). Of interest, the top 15 Six2-regulated genes (based on highest expression) are all associated with roles in normal and CSC function and/or tumor progression (40–51). Importantly, high expression of the Six2-gene signature significantly correlates with a decrease in distant metastasis-free and relapse-free survival across all breast cancer patient samples as well as in TNBC-specific samples (Fig. 7B and C). Additionally, the Six2-gene signature significantly correlates with a decrease in time to distant metastasis, distant metastasis-free, and recurrence free survival across all breast cancer (we were unable to query TNBC specifically in this dataset; Fig. 7D). These data demonstrate that the relationship between SIX2 and metastatic burden is not limited to our cell line models, and that the SIX2 transcriptional program is predictive of poor outcomes in human breast cancer.

Previous data from our group, using two different triple negative mouse mammary carcinoma xenograft models, established Six2 as critical for late stage metastasis (9). In particular, the fact that Six2 results in outgrowth of 4T07 cells, which are known to reach the lungs after orthotopic injection into the mammary gland, but to not efficiently colonize the lungs (15), demonstrates a role for Six2 in the colonization step of the metastatic cascade. Further, data within this manuscript demonstrate that SIX2 also regulates late-stage metastasis in human TNBC.

Although our data suggest that a main function of SIX2 is to mediate outgrowth of cells that have newly arrived at a secondary site, given its described role in mediating mesenchymal phenotypes and anoikis resistance (9, 14), SIX2 may in fact affect several additional aspects of late-stage metastasis, although this remains to be formally tested. If SIX2 regulates multiple steps of late stage metastasis (such as survival in the bloodstream and exit from the vasculature in addition to colonization at secondary sites), this finding would further emphasize the importance of developing novel means to target its function, as significant data suggest that metastases can seed additional metastatic sites (52). In the future, it will be important to determine whether SIX2 is required for maintenance of metastatic lesions. This question is critical to address, as patients that already present with significant metastatic burden would benefit from a therapy that inhibits the continued growth of cells at the secondary site. Nonetheless, even if SIX2 is only critical for the establishment and initial stages of metastatic outgrowth, many patients with breast cancer do not present with overt metastases, and thus inhibiting this step could be valuable. However, as SIX2 is a transcription factor, its direct inhibition may prove difficult. Thus, understanding the mechanism by which SIX2 mediates late stage metastasis may enable the identification of novel downstream targets or pathways that can potentially be targeted in the future.

Previous studies have shown that Six2 maintains the kidney progenitor pool during nephrogenesis, in part through impinging on stem cell-associated pathways (7). Based on our RNA-seq analysis, we find that similarly, in the context of breast cancer, SIX2 regulates a genetic stem-cell like program. We demonstrate that Six2 OE increases, whereas SIX2/Six2 KD decreases, the expression of multiple genes that are enriched in normal stem cells. These data suggest that SIX2 may regulate stem cell-like properties across multiple tissue types and across normal and disease states. Interestingly, although it has mostly been studied in the context of nephrogenesis (7), previous studies from our laboratory have shown that during mammary gland development, Six2 expression is increased in the mammary gland of Six1 KO mice (53). Because Six1 OE in the mammary gland results in the expansion of mammary stem/progenitor cells (8), yet Six1 loss does not affect mammary gland development (53), it is possible that upregulation of Six2 compensates for Six1 loss to maintain mammary stem cell function. Similar to SIX1, SIX2 is rarely expressed in adult tissues, but it is re-expressed in the context of breast cancer (9). The role of SIX2 in nephrogenesis, coupled with its re-expression in human breast cancers, suggests that SIX2 may be a key regulator of both mammary stemness and mammary tumorigenesis and/or metastasis. However, the functional role of SIX2 in the normal mammary stem cell compartment and in mammary gland development has yet to functionally explored. Analysis of mRNA expression from isolated mouse fetal mammary stem and adult mouse differentiated mammary epithelial cells (23) shows that Six2 mRNA is more highly expressed in the fetal mammary stem cell population compared with the adult mammary cell population, suggesting that Six2 may play a role in fetal mammary stem cells (Supplementary Fig. S6). Thus, investigation into whether Six2 contributes to normal mammary development may be warranted.

Previous studies using ChIP-seq analysis in the developing kidney identified Gdnf, Eya1, and Six2 itself as direct transcriptional Six2 targets (7, 38). Interestingly, our RNA-seq analysis in 4T07 cells identifies Eya1 as one of the top upregulated genes in response to Six2 OE, suggesting conservation of specific downstream targets of Six2 between normal development and cancer. As outlined above, we have identified Sox2 as a novel direct transcriptional target of Six2 in breast cancer cells that is critical for Six2-mediated metastasis. We further demonstrate that Six2 directly regulates Sox2 transcription by binding to, and activating, the Srr2 enhancer region of the Sox2 gene. Thus, it will be of interest to determine whether Six2 similarly regulates Sox2 in the context of normal embryonic development. Although it is important to continue to understand how SOX2, as a key stem cell factor, mediates SIX2-regulated late-stage metastasis, unfortunately, SOX2, like SIX2, is a transcription factor and thus not an ideal drug target. Moving forward, it will be important to identify more targetable downstream effectors of SIX2 that are critical in mediating its ability to promote metastatic burden.

Overall, our studies indicate that SIX2 plays an important role in promoting cancer stem cell–associated properties through a novel SIX2/SOX2 axis that is critical for metastatic colonization. Importantly, the relationship between SIX2 and SOX2 is relevant in human breast cancer, and high expression of a Six2-mediated gene signature significantly correlates with poor metastatic prognosis. Although it is well established that targeting cancer stem cells at the primary tumor site can increase therapeutic efficacy (54), our findings support growing evidence that these strategies should be considered when targeting metastatic tumor cells, in which there may be different regulators of tumor outgrowth that provide an advantage in a new microenvironment (3). Our results support this emerging concept by providing critical knowledge regarding how stem cell properties in newly arriving tumor cells promote late-stage metastasis.

No potential conflicts of interest were disclosed.

Conception and design: M.U.J. Oliphant, J. DeGregori, H.L. Ford

Development of methodology: M.U.J. Oliphant, V. Zaberezhnyy, J.C. Costello, J.M. Espinosa, H.L. Ford

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.U.J. Oliphant, K.R. Johnson, J.C. Costello, H.L. Ford

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.U.J. Oliphant, M.Y. Vincent, M.D. Galbraith, A. Pandey, P. Rudra, J.C. Costello, D. Ghosh, J.M. Espinosa, H.L. Ford

Writing, review, and/or revision of the manuscript: M.U.J. Oliphant, M.Y. Vincent, M.D. Galbraith, J.C. Costello, D. Ghosh, J.C. Costello, H.L. Ford

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.U.J. Oliphant, M.Y. Vincent, V. Zaberezhnyy, J.M. Espinosa

Study supervision: M.U.J. Oliphant, H.L. Ford

Other (first author of the manuscript and performed the bulk of the experiments.): M.U.J. Oliphant

Other (graphical abstract): M.Y. Vincent

We thank Yasushi Okada for providing the pGF-mSox2-ONL and pCNL-N1_3xNLS plasmids for use in the SRR2 luciferase reporter assays. We also thank the University of Colorado Cancer Center Animal Imaging (P30CA046934), Flow Cytometry (P30CA046934), Functional Genomics (P30CA046934), and Genomics Shared Resources (P30CA046934) for aiding in experimental analysis used within this work. This work was supported by R01CA180175 (to J. DeGregori), R01CA095277 (to H.L. Ford), The Cancer League of Colorado - 25A5976 (to H.L. Ford), R01CA117907 (to J.M. Espinosa), the NCI F31 Ruth L. Kirschstein National Service Research Award 5F31CA210622-02 (to M.U.J. Oliphant), and the NCI F99/K00 Fellow Transition Award F99CA223023 (to M.U.J. Oliphant).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.