FOXA1 is Prognostic of Triple Negative Breast Cancers by Transcriptionally Suppressing SOD2 and IL6.

Having markers feasible for breast cancer subtyping, especially for triple negative breast cancer identification is crucial for improving the treatment outcome of such cancers. Here we explore the role of FOXA1 in characterizing triple negative breast cancers and the driving mechanisms. Through in vitro examination of the expression pattern at both transcriptional and translational levels, patient relapse-free survival analysis, immunohistochemistry staining and prediction power assessment using clinical samples, as well as functional studies, we systematically compared the role of FOXA1 in identifying triple negative and luminal type of breast cancers and explored the mechanisms driving such functionalities. We report that FOXA1 under-expression can lead to increased malignancy and cancer stemness, and is a subtyping marker identifying triple negative breast cancers rather than the luminal subtype by transcriptionally suppressing the expression of SOD2 and IL6. We are the first to systematically address the significance of FOXA1 in triple negative breast cancer identification as a biomarker and elucidate the mechanism at the molecular level, through a sequential bioinformatics analysis and experimental validations both in vitro and in clinics. Our discoveries compliment the current biomarker modalities once verified using larger clinical cohorts and improve the precision on characterizing breast cancer heterogeneity.

Introduction

As the most common cancer in women worldwide 1 , breast cancer is not a single disease. Perou et al. firstly divided breast cancers into four intrinsic subtypes, i.e., luminal, HER2 positive, basal, and normal-like tumors 2. Sorlie et al. have sub-divided luminal tumors into A and B subtypes 3. A lot more molecular subtypes have been identified later on, with triple negative breast cancers (TNBCs) being the most heterogeneous and morphologically diverse encompassing at least basal, claudin-low, metaplastic breast cancers and the interferon-rich subtypes 4. Different molecular subtypes have distinct clinical implications, with TNBC being the most aggressive that still lacks efficient targeted therapies 5. How to appropriately identify TNBCs and understand the mechanism driving its malignancy is of great importance in improving the diagnosis and therapeutics of TNBCs.

We have previously identified FOXA1 as a potential marker for breast cancer subtyping through a series of bioinformatics analysis6, 7. By in vitro screening of FOXA1 and its correlated genes using 10 breast cancer cell lines covering four subtypes at both gene and protein expression levels followed by a series of computational verifications, functional studies and clinical validations, we propose that low FOXA1 expression is associated with TNBCs, and it functions as a transcriptional suppressor of SOD2 and IL6 to contribute to the invasive and stem-like features of TNBCs. By systematically comparing the performance of FOXA1 in characterizing TNBC and luminal tumors, we propose it as a marker highly associated with TNBC, which contradicts with the canonical conception that FOXA1 is representative of ER and associated with luminal type of cancers 8 , and elucidate the driving mechanism in vitro. The outputs of this study enrich the toolbox for the diagnosis and therapeutics of TNBCs.

Materials and methods

Experimental materials and protocols

Breast cancer cell lines

A total of 10 breast cancer cell lines, obtained from American Type Tissue Culture Collection (Manassas, VA, USA), were used in this study (Supplementary ). These cell lines were authenticated through monitoring of cell morphology, and were passaged in our laboratory for fewer than 6 months after receipt or resuscitation.

Quantitative RT-PCR

Total RNA was extracted using RNA simple Total RNA Kit (TIANGEN, China), and complementary DNA (cDNA) was synthesized from total RNA using the PrimeScript PT reagent Kit with gDNA Eraser (TaKaRa, China) following the manufacturer's protocol. The primers designed for each gene are listed in Supplementary . Real time PCR was conducted on the ABI PRISM 7500 Quantitative PCR system (Life Technologies, USA) using the SYBR Green Real time PCR master mix (QPK-201, Toyobo). Each sample was examined in triplicates. The relative gene expression levels were normalized to GAPDH (internal control) using the 2-△△ Ct method.

Western blot

Cells were lysed with ice-cold RIPA buffer with 0.5 mM Phenylmethanesulphonyl fluoride (PMSF). The protein content was determined according to Bradford's method. After separation using SDS-PAGE gel, proteins were transferred to PVDF membrane and probed with human antibodies against FOXA1 (1:5000, Abcam, UK), SOD2 (1:5000, Abcam), MYC (1:1000, Proteintech, US), HER2 (1:1000, Proteintech), IL6 (1:1000, Proteintech) and GAPDH (1:4000, Proteintech).

Cell proliferation and apoptosis

Cell viability was assessed by a Cell Counting Kit 8 (CCK-8, Japan) according to the manufacturer's protocol. The apoptosis rates were identified using an Annexin V-FITC and PI Apoptosis Detection Kit (BD Biosciences, USA). Cells were detected using the FACS scan flow cytometer, and the data was analyzed using the Flowjo software.

Gene knockdown by siRNA transfection

Breast cancer cells were transfected with human-specific FOXA1, SOD2, IL6 or MYC siRNA, GAPDH siRNA (for optimization), and non-silencing siRNA (negative control siRNA) (Gene Pharma, China) using the siRNA-mate transfection agent (Gene Pharma). In addition, to avoid off-target effects of siRNAs, we used another siRNA sequence for each gene. The sequences of siRNAs for FOXA1, SOD2, and Myc are shown in Supplementary

Gene up-regulation by CRISPR editing

FOXA1 was overexpressed using CRISPR/ dCas9 Synergistic Activation Mediator (SAM) system following protocols described previously 9. Three sgRNAs targeting FOXA1 (sequences in Supplementary ) were concatenated and cloned into one plasmid (Synbio Technologies, China) followed by co-transfection with the dCas9 Synergistic Activation Mediator Lentivector (Applied Biological Materials Inc, Canada) into BT474 using Lipofectamine 2000 (Invitrogen, USA). Positive cells were selected using G418 disulfate salt (300ug/ml) and Puromycin (0. 25ug/ml).

Cell migration detection by transwell

Transfected and non-transfected cells were incubated for 48 hours under normoxic and anaerobic conditions, respectively. Cell medium was added on the lower layer of 24-well culture plate and the chambers were placed in the medium. Cells were collected following pancreatic digestion, re-suspended and added to the chambers (2×105/well). The culture media inside the chambers were discarded after 20 hours, and cells were washed by PBS (phosphate buffered saline). Migrated cells under the chambers were fixed by methanol followed by staining with 0.1% crystal violet solution.

ALDEFLUOR assay and separation of the ALDH positive population by FACS

ALDEFLUOR assays were performed according to the manufacturer's instructions (Stem Cell Technologies, Durham, NC, USA). In brief, 2.5×105 cells were suspended in 500 μL ALDEFLUOR assay buffer containing 5 μL/mL ALDEFLUORTM substrate and incubated for 30 minutes at 37 °C in darkness. As a negative control, cells were stained under identical conditions in the presence of the specific ALDH inhibitor diethylaminobenzaldehyde (DEAB). After 30 minutes, cells were centrifuged, the supernatant was removed and the remaining pellet was suspended in ice-cold ALDEFLUORTM assay buffer and kept on ice. Cells were immediately assayed with FACS Calibur (Becton Dickinson Biosciences, Franklin Lakes, NJ, USA) using DEAB controls as baselines to gate ALDH+ and ALDH- cell populations.

Mammosphere formation assay

Mammosphere formation assays were performed to determine the sphere-forming activity of cancer stem cells (CSCs) as previously described10. Briefly, single-cell suspensions prepared from human SKBR3 cells (with or without being supplemented with IL6) were cultured at 2000 to 5000 cells/mL per well in 24-well ultra-low attachment plates (Corning) using serum-free DMEM/F-12 medium supplemented with 20 ng/mL basic FGF, 20 ng/mL EGF, 4 μg/mL insulin, 4 μg/mL heparin, 0.5 μg/mL hydrocortisone, 0.4% BSA and B27 (Invitrogen). Culture medium was replaced every other day with 50% fresh medium. Tumor spheres were counted and photographed after 7 days of culture. Cells forming tumor spheres were harvested and cultured as single clones to examine their ability of forming secondary tumor spheres following the same procedures.

Chromatin immunoprecipitation assay

Chromatin immunoprecipitation (ChIP) assay was performed according to the manufacturer's protocol (Beyotime, China) with slight modifications. Chromatin solutions were sonicated and incubated with a monoclonal goat anti-human FOXA1 antibody (0.02 μg/μL; Abcam) or control IgG overnight at 4℃. DNA-protein cross-links were reversed and chromatin DNA was purified and subjected to PCR analyses (primers are in Supplementary ). After amplification, PCR products were resolved using 3% agarose gel and visualized by ethidium bromide staining.

Luciferase reporter assay

The pGL3 basic plasmids with or without adding the SOD2 and IL6 promoter sequence were co-transfected with the plasmid expressing Renilla luciferase (internal control) to BT474 cells. Thirty-six hours after transfection, cells were harvested and the luciferase activities were measured using a Luciferase Reporter Gene Assay Kit (Beyotime, China).

Clinical sample collection, IHC staining and statistical analysis

A total of 82 human primary breast cancer tissues were collected at Affiliated Hospital of Jiangnan University from years 2008 to 2012 with informed consent (Supplementary ), and this project was approved by the Clinical Research Ethics Committees of Affiliated Hospital of Jiangnan University.

IHC staining was performed on 4 μm sections from formalin fixed paraffin embedded breast cancer tissues using anti-FOXA1 antibody (1:100, Abcam) and anti-SOD2 antibody (1:200, Abcam). FOXA1 nuclear expression was scored based on the staining intensity and the positive percentage of tumor cells as previously reported with some modifications 11. Detailed evaluation method is described in Supplementary . For SOD2, IHC staining was observed in the cytoplasm of tumor cells, and a scale of 0 to 3 was used to score relative expression intensity. Chi-square test was conducted to examine the significance of FOXA1 or SOD2 IHC status in differentiating tumors of different subtypes.

Data and statistical analysis

Public dataset

METABRIC, TCGA (level 3), GSE24450 and E-MTAB-181 data were used for statistical analysis (Supplementary ).

Statistical analysis

Gene expression was stratified by breast cancer subtypes to evaluate the subtyping performance of potential markers using clinical samples. Student T test was used to assess the statistical significance of dichotomizing patients into TNBC and non-TNBC cohorts.

KM plotter 12 (http://kmplot.com/analysis/) was used to screen genes with significant association with patient survival (Supplementary ). KM plotter binarizes patients by the median expression of a gene. Five years relapse free survival (RFS) was conducted using 3951 patients (the number of available patients for conducting RFS analysis). The hazard ratio and p value from log rank test were computed for each gene.

Pearson correlation analysis was conducted to evaluate gene correlations, with p values from Pearson's Chi-squared test being used to access the statistical significance.

The receiver operating characteristic (ROC) curve was drawn to access the precision and recall of the performance of a marker as a classifier, where higher area under the curve (AUC, defined as the definite integral of the curve) represents a better marker performance.

The fitness of the markers in modeling TNBC or luminal breast cancers as independent cohorts was evaluated by fitting marker gene expression and clinical subtyping information into a linear model. The adjusted R2 and p values computed from Chi-squared test were used to access the statistical fitness of the model.

Results

FOXA1 identifies triple negative cells and exhibits an opposite expression pattern with SOD2 and IL6 in breast cancer cells

Through conducting Pearson correlation analysis using METABRIC, TCGA, GSE24450 (r>0.5, p<0.05) and Patients' 5-year RFS analysis using KM plotter (which includes 3951 patients), we found 18 genes transcriptionally associated with FOXA1 and being prognostic of patient clinical outcome with statistical significance (Supplementary ). Additionally, IL6 could induce the percentage of breast CSCs and mammosphere formation (Figures ) and is a verified transcriptional target of IL6DBP that is among the candidates 13.

We found that 9 out of the 19 candidates were capable of differentiating breast cancer subtypes with statistical significance at the transcriptional level (Figure , Supplementary ), and three out of them (FOXA1, SOD2, and IL6) could be verified at the protein expression level (Figure . Such patterns were also observed at the transcriptional level using E-MTAB-181 data set (Figure ).

It was interesting that FOXA1 exhibited an opposite pattern with SOD2 and IL6, both of which characterize triple negative cells. That is, while SOD2 (p=0.0132 for mRNA expression, p=0.0118 for protein expression) and IL6 (p=0.0277 for mRNA expression, p=0.005 for protein expression) were over-represented in the triple negative subtype, FOXA1 (p=0.0005 for mRNA expression, p=5.21e-6 for protein expression) was under-expressed in this type of tumor cells (Figure ).

An exception was observed in SKBR3, a HER2 positive cell line, which exhibited low FOXA1 level and slightly over- expressed SOD2 and IL6 (Figure ). BRCA1 mutation status did not alter the expression of FOXA1 nor IL6 with statistical significance (Figures ), but further up-regulated SOD2 expression (p=0.0335 for gene expression, p=0.0433 for protein expression as compared with the rest TNBC cells, Figures ).

FOXA1 is a subtyping marker associated with TNBCs rather than luminal tumors using clinical samples

IHC staining of 82 collected breast cancer samples including 29 TNBCs demonstrated that FOXA1 was significantly down-regulated in TNBCs (Figure ), and its subtyping potential was substantially higher in distinguishing TNBCs rather than luminal cancers from the other subtypes (p=3.52e-13 versus p=2.22e-03 for identifying TNBCs and luminal cancers, respectively, Table ). In particular, 82.72% of TNBCs exhibited low FOXA1 expression (below score 2), whereas 68.75% of luminal cancers expressed high levels of FOXA1 (score≥2). The precision and recall of FOXA1 as a subtyping marker was evaluated using AUC, which were 0.8307 and 0.7477, respectively, for discriminating the TNBC and luminal subtypes from the rest breast cancers (Figure ). We obtained similar results using public datasets METABRIC (AUC= for 0.9743 TNBC versus AUC= for 0.9085 luminal tumor identification, Figure ) and TCGA (AUC=0.9635 for TNBC versus AUC=0.9208 for luminal tumor identification, Figure ).

Patients' 5 years RFS analysis from 3951 breast cancer samples deposited in KM plotter suggested that low mRNA expression of FOXA1 was significantly associated with poor patient outcome in clinics (p=1.7e-11, Figure ). However, such an association diminished when samples were stratified by molecular subtypes (Supplementary ), further consolidating our claim that FOXA1 is a molecular marker for categorizing tumors of different subtypes but not for differentiating tumors within the same subtype. However, we could not completely exclude the possibility that the insignificance is caused by the limited number of patients available in certain subtype cohorts.

By grouping 502 breast cancer patients from TCGA into luminal, HER2 positive and triple negative subtypes, we found that the expression of FOXA1 in triple negative cancers was significantly lower than that in the other subtypes (p=1.19e-97, Figure ). More significant results were observed using METABRIC data (p=5.69e-74, Figure ).

Using METABRIC and TCGA data which have relatively large sample sizes (i.e., METABRIC encompasses 1217 patients and TCGA has 502 samples), we examined the fitness of using FOXA1, alone or in combination with the canonical marker panel (ER and HER2) in modeling TNBC subtype as an independent cohort. FOXA1 had a comparable fitness with the combined marker panel (adjusted R2=0.6844 and p=2.98e-308 for FOXA1, adjusted R2=0.6909 and p=8.38e-311 for combined panel, Table ) that was considerably higher than using ER and HER2 (adjusted R2=0.5146 and p=8.82e-193, Table ) from the analysis of the METABRIC data. Similar results were obtained from TCGA data. That is, the adjusted R2 was 0.6023 (p=2.35e-101), 0.5535 (p=4.97e-91), 0.5381 (p=4.99e-86) for the combined panel, FOXA1, and ER-HER2 combination respectively. It was interesting that as FOXA1 outperformed canonical ER-HER2 in discriminating TNBCs, ER-HER2 was better than FOXA1 in capturing the luminal feature of breast cancers. For example, the adjusted R2 were 0.6604 (METABRIC) and 0.6205 (TCGA) using ER-HER2, which became 0.4494 (METABRIC) and 0.4708 (TCGA) using FOXA1. Therefore, while the canonical ER-HER2 classifier mostly captures the properties of luminal cancers, low FOXA1 expression is a potential independent prognostic marker of TNBCs, which can be either integrated into the ER-HER2 panel for improved accuracy in breast cancer subtyping or used to substitute ER-HER2 for TNBC identification with enhanced precision and reduced complexity.

SOD2 is a weak subtyping marker associated with TNBCs and IL6 does not have diagnostic value using clinical samples

Both METABRIC and TCGA data showed that SOD2 expression can significantly distinguish TNBCs (p=4.25e-52 for METABRIC, p=1.38e-07 for TCGA, Supplementary ), despite the reduced discrimination power as compared with FOXA1. High SOD2 gene expression was significantly prognostic of poor patient outcome in the 5 years RFS analysis using KM plotter (p=5.50e-10, Supplementary ). However, only a slightly higher expression of SOD2 was observed in TNBCs as compared with samples of the other subtypes using our 82-patient cohort (Supplementary , Supplementary ), suggesting that SOD2 is a weak subtyping maker associated with TNBCs.

No varied expression of IL6 across cancer subtypes was obtained using neither METABRIC nor TCGA data (Supplementary ). No distinctive effect of IL6 on 5 years' patient RFS from KM plotter (Supplementary ) was observed. However, IL6DBP over-expression was associated with poor patients 5 years' RFS (p=8.6e-08, Supplementary ). Thereby, while IL6 could distinguish TNBCs in vitro (Figures ) and was transcriptionally regulated by IL6DBP that was associated with patient outcome (Supplementary ), its diagnostic value was not visible in clinics given the complex tumor microenvironment and diversified roles of chemokines in living tissues.

FOXA1 inhibits cancer cell proliferation and migration, SOD2 and IL6 play the opposite roles

Knocking down FOXA1 decreased apoptosis (Figure p=0.0013 for MCF7; p=0.0240 for BT474) and accelerated cell proliferation (Figure p=0.0016 for MCF7, Supplementary ; p=0.0064 for BT474, Supplementary ). Oppositely, SOD2 knockdown promoted apoptosis (Figure p=0.0023 for SUM149PT; p=0.025 for MDAMB436) and inhibited cell growth (Figure p=0.0002 for SUM149PT, Supplementary ; p=0.0006 for MDAMB436, Supplementary ); and silencing IL6 led to inhibited cell growth (Figure p=0.0003 for SUM149PT, Supplementary ; p=0.0350 for MDAMB436, Supplementary ).

Enhanced motility (p=0.0014 for MCF7, p=0.0020 for BT474) and increased proportion of stem cells (p=0.0001 for MCF7, p=0.0004 for BT474) were observed in FOXA1-silenced MCF7 and BT474 cells (Figures ). On the contrary, knocking down SOD2 inhibited cell migration (p=0.0013 for SUM149PT, p=0.0002 for MDAMB436) and decreased the amount of CSCs (p=0.0001 for SUM149PT, p=0.0002 for MDAMB436,) in these luminal cells (Figures ). Similar behavior was observed for cells lacking IL6 expression, i.e., cell migration was significantly inhibited (p=0.0020 for SUM149PT, p=0.0003 for MDAMB436) and the percentage of stem cells largely decreases (p=0.0009 for SUM149PT, p=0.0006 for MDAMB436) in IL6-silenced SUM149PT and MDAMB436 cells (Figures ).

FOXA1 is a transcription factor of SOD2 and IL6

Transcription factor (TF) binding site prediction was conducted using JASPAR (http://jaspar.genereg.net/) 14 to investigate whether FOXA1 could regulate SOD2 and IL6 as a TF. FOXA1 was shown to be the most probable TF of SOD2 and IL6 among all TF binding motifs in vertebrates (the top scores of SOD2 and IL6 are 12.195 and 10.2481, respectively).

The ChIP assay performed using BT474 cells (cells expressing high level of FOXA1) revealed that FOXA1 could bind the promoter region of SOD2 and IL6 (Figure )15. Silencing FOXA1 (p=4.17e-07) significantly enhanced SOD2 and IL6 expression at both transcriptional (p=9.82e-05 for SOD2, p=0.0059 for IL6, Figure ) and translational (Figure ) levels; over-expressing FOXA1 using CRISPR method (p=6.27e-05) showed the reverse pattern, both reaching statistical significance at the mRNA (p=4.82e-05 for SOD2, p=0.0013 for IL6, Figure ) and protein expression levels (Figure ).

Luciferase reporter assay further validated our claim that FOXA1 physically binds the promoter region of SOD2 and IL6 and suppresses their expression (Figure ). That is, by adding the promoter sequences of SOD2 or IL6 to the basic plasmid, the expression of these genes enhanced (p=6.92e-07 and p=5.86e-06, respectively, for SOD2 and IL6, respectively), which were further increased with statistical significance (p=0.0001 for SOD2, p=1.75e-05 for IL6) when FOXA1 was down-regulated (p=0.0002).

We further analyzed the correlations between the expression of FOXA1 and SOD2, FOXA1 and IL6, and SOD2 and IL6 using clinical data from METABRIC, TCGA and GSE24450. The gene expression of FOXA1 showed a negative correlation with SOD2 and IL6 (average correlation coefficient is -0.61 and -0.32 for FOXA1-SOD2 and FOXA1-IL6, respectively), and SOD2 and IL6 were transcriptionally positively correlated using all these public datasets (average correlation coefficient is 0.3, Supplementary ).

The Pearson correlation score between FOXA1 and SOD2 from our 82-sample cohort was -0.4129 (p=1.16E-4), which dropped to -0.3562 (p=1.02E-3) once the data was discretized by score 2. Though some information was lost due to data discretization, this adds further evidence to the negative regulatory relationship between FOXA1 and SOD2.

FOXA1 influences TNBC features by transcriptionally regulating SOD2 and IL6

FOXA1, SOD2 and IL6 were effectively knocked down in MCF7 (p=0.0026 for FOXA1, p=0.0070 for SOD2, p=0.0225 for IL6). The proliferation of FOXA1-low cells was significantly reduced (p= 0.0020) and the apoptosis of FOXA1-low cells was significantly increased (p=0.0040) after knocking down SOD2 and IL6 despite the fact that FOXA1 was under-expressed (Figure ).

Discussion

FOXA1 is a subtyping marker for TNBC identification

FOXA1 is known indispensable for the expression of 50% of ER-related genes 8, 16 and implicated in various cancers including breast, prostate, lung, thyroid and esophageal squamous cell carcinomas 17. There is a strikingly high Pearson correlation between ER/FOXA1 and apoptosis- related genes 18, suggesting its tumor suppressive role. However, the functionalities of FOXA1 have largely been associated with that of ER and luminal features, very rare has been reported to associate FOXA1 with TNBCs 19-21. We found from our in vitro and clinical results (using both public datasets and our collected sample cohort) that FOXA1 significantly outperformed ER-HER2 (canonically applied in clinics for breast cancer subtyping) in TNBC identification and could be used alone to achieve enhanced accuracy and reduced complexity than the conventional diagnostic panel. Alternatively, FOXA1 can significantly enrich information captured by the ER-HER2 panel, which can be jointly applied for breast cancer subtyping with considerably enhanced diagnostic accuracy. These findings are not mutually exclusive from the previous reports that FOXA1 is a luminal differentiation factor 22, but provide a more accurate definition of FOXA1 towards its subtyping functionalities, i.e., excluding HER2 positive breast cancers from non- luminal tumors.

The exceptional observation of low FOXA1 expression in a HER2 positive cell line (SKBR3) that harbors an MYC amplification warrants us that the diagnostic accuracy of FOXA1 may be influenced by MYC mutation. This suggests a combinatorial use of FOXA1 and MYC in the diagnostic practice, i.e., favorable clinical outcome is expected for patients carrying high FOXA1 and low/normal MYC expression but not for cases otherwise.

FOXA1 contributes to TNBC phenotypic features through suppressing SOD2 and IL6

The discovery of FOXA1 as a prognostic marker of TNBCs uncovers its mechanistic functionalities in driving the invasiveness and stem-like features of TNBCs. Through TF predictions followed by ChIP and luciferase reporter assay, we identified FOXA1 as a transcriptional suppressor of SOD2 and IL6. As FOXA1 interferes with the recruitment of NF-κB to the promoter regions of IL6 15 and SOD223, it is possible that FOXA1 suppresses NF-κB signaling by competing with NF-κB for the binding sites of downstream effectors such as IL6 and SOD2, leading to halted tumor migration and reduced cancer cell stemness (Figure ).

SOD2 is a known anti-apoptotic factor diluting the detrimental effects of reactive oxygen species (ROS) to cancer cells. Via suppressing SOD2 expression, FOXA1 is likely to function as a safe guard controlling cell life-death switch. This is in accordance with our observation that cells underwent faster proliferation with less apoptosis when FOXA1 was knocked down, and the reverse pattern was observed once SOD2 was silenced (Figure ). SOD2 expression was over-represented in TNBCs and, in particular, in tumors with concomitant BRCA1 mutation (namely the basal subtype, Figure ). That is, with enhanced ability of bypassing apoptotic signals (as empowered by SOD2 over-expression) and reduced capacity in repairing DNA damage (as enabled by BRCA1 deficiency24), cancer cells undergo accelerated evolutionary alterations that ultimately transit them towards a more malignant state, e.g., being invasive 25.

IL6 encodes a cytokine that is implicated in CSC maintenance and progenitor-enriched mammosphere formation 26-28 (Figures ). Increased stemness could ultimately lead to increased invasiveness of cancer cells. Yang et al. found that the IL6/CXCR2 axis played a critical role in the metastasis of oral squamous cell carcinoma 29, and IL-6 activation promoted cervical tumorigenesis through autocrine and paracrine pathways in tumor microenvironment 30. Thus, FOXA1 may halt cancer metastatic transition through suppressing IL6 that contributes to the CSC properties of TNBCs.

Even FOXA1 transcriptionally suppresses SOD2 and IL6 expression, and modulating FOXA1 and SOD2/IL6 leads to opposite cell behaviors, we cannot claim that FOXA1 influences TNBC features by transcriptionally regulating SOD2 and IL6. Given that FOXA1 is a known proliferation marker31, we examined cell proliferation and apoptosis in the recovery test. Concomitantly knocking down SOD2 and IL6 in FOXA1-low MCF7 cells restores cells with the luminal features (Figure ), suggesting that FOXA1 likely exerts its roles via controlling the expression of SOD2 and IL6, which are the direct executioners regulating the phenotypic switch between triple negative and luminal features.

MYC suppresses FOXA1 expression in HER2 positive cells - Explanation of the exceptional in vitro observation

An exceptionally low FOXA1 expression and slightly increased levels of SOD2 and IL6 were observed in the SKBR3 cell line (Figures ) that is HER2 positive and harbors an MYC amplification32. This suggests strong signaling interventions among FOXA1, MYC and HER2 that partially interfere with SOD2 and IL6. We hypothesized a triangle negative feedback loop (). That is, FOXA1 binds the promoter region of HER2 and elevates its expression, increased level of HER2 leads to decreased MYC expression, and MYC binds the promoter region of FOXA1 and suppresses its transcription, formulating a negative feedback loop. This concords with the report that FOXA1 regulates the transcription of a panel of genes in HER2 signaling, including reversely regulated MYC 33, and that MYC is a downstream effector of HER2 34. It is predicted from Jasper (http://jaspar.genereg.net/) 14 that FOXA1 is a potential TF of HER2 (the top binding affinity score is 10.721), and MYC potentially regulates FOXA1 expression (the highest binding affinity score is 14.020).

The results showed that silencing MYC led to increased FOXA1 and HER2 expression (), and knocking down FOXA1 repressed HER2 but increased MYC expression () in SKBR3, which has a higher expression of MYC than MDAMB453, another HER2 positive cell line ().

The rest genes imply other candidate markers and targets for breast cancer management

The rest genes (Supplementary ), which are transcriptionally associated with FOXA1 and prognostic of patient clinical outcome with statistical significance, also implicate candidate markers for breast cancer subtyping and worth further exploration.

FBP1 loss is required for cancer cells to acquire the basal-like phenotype 35. The tumor suppressive role of FBP1 has been reported in various cancers such as non-small-cell lung cancer, a hepatocellular carcinoma, colon cancer, and pancreatic cancer 36-38. Our results showed that FBP1 could differentiate luminal tumors with a superior performance than ER (p=0.0001 for FBP1, p=0.0021 for ER), suggesting its subtyping potential.

AGR2 and MYO5C were over-expressed in the luminal A subtype with statistical significance (p=0.0286 for AGR2, p=0.0001 for MYO5C). AGR2 is an estrogen-responsive gene that is positively correlated with ER at the transcriptional level 39. It promoted tumor growth and metastasis 40, 41, and was suggested as a potential drug target or biomarker for various cancers including colorectal 42, gastric 43 and breast cancers 40, 44. MYO5C functions in the trafficking of integral membrane proteins to melanomas 45. Our results implicate their tumorigenic roles that worth further investigations.

LRBA and WWP1 were over represented in luminal B cells (p=0.0001 for both LRBA and WWP1). LRBA promoted cancer cell growth and was associated with p53 and Rb mutations 46. WWP1 was considered an oncogenic factor in prostate 47, breast and oral cancers 47-49. Further explorations on these genes may uncover the discrepancies between luminal A and luminal B types of breast cancers and the mechanism driving such a differentiation. Also, if confirmed as markers specific to the luminal B type of cancers (provided with more luminal B cell lines and clinical evidence), these may represent novel therapeutic targets besides their subtyping potential.

Conclusion

This study proposes FOXA1 as an independent subtyping marker for TNBCs identification through bioinformatics analysis followed by in vitro and clinical validations. Also, we propose that FOXA1 likely halts the triple negative feature of cancer cells by transcriptionally suppressing SOD2 that helps cells bypass DNA-damage-induced apoptosis and inhibiting IL6 that enables cells with stem-like features and invasive nature. That is, by impeding cells at two critical transitions towards carcinogenesis, i.e., life- death control and metastatic switch, FOXA1 plays a tumor suppressive role and is under-expressed in TNBCs. We are the first to propose FOXA1 as a TNBC identification marker and elucidate the potential mechanism, which are of the guiding significance in the precise control of TNBCs.

Supplementary figures and tables.

Click here for additional data file.

Table 1

Results showing the fitness of FOXA1, ER-HER2 and their combined panel in modeling triple negative or luminal breast cancers as independent cohorts.

Variables	Data	TNBC vs. non-TNBC			luminal vs. non-luminal
		Combined	FOXA1	ER-HER2	Combined	FOXA1	ER-HER2
		ER, HER2, FOXA1	FOXA1	ER, HER2	ER, HER2, FOXA1	FOXA1	ER, HER2
Adjusted R2	METABRIC	0.6909	0.6844	0.5149	0.6838	0.4494	0.6604
p (model)		8.376e-311	2.98e-308	8.82e-193	7.96e-305	1.22e-160	3.09e-287
p (ER)		1.62e-07	-	5.59e-178	2.50e-99	-	4.07e-286
p (HER2)		5.17e-03	-	3.60e-53	6.44e-08	-	7.48e-01
p (FOXA1)		1.08e-121	2.98e-308	-	5.56e-21	1.22e-160	-
Adjusted R2	TCGA	0.6023	0.5535	0.5381	0.6294	0.4708	0.6205
p (model)		2.35e-101	4.97e-91	4.99e-86	4.38e-109	2.67e-72	1.26e-107
p (ER)		1.86e-11	-	6.93e-73	1.92e-40	-	3.43e-106
p (HER2)		2.43E-09	-	6.95e-23	3.14e-02	-	4.66e-05
p (FOXA1)		2.13E-18	4.97e-91	-	3.11e-04	2.67e-72	-

The results were produced by fitting data (METABRIC and TCGA data) to a linear model. Adjusted R2 and the p values of the model as well as that of each predictor were used to evaluate the fitness of each model.

Table 2

IHC staining results of FOXA1 in 82 breast cancer tissue samples.

Subtype	N	FOXA1 expression		p value
		Score 0-1	Score 2-3
TNBC	29	24	5	3.52e-13
non-TNBC	53	11	47
Luminal	32	10	22	2.22e-03
non-Luminal	50	25	25