FOXQ1 is Differentially Expressed Across Breast Cancer Subtypes with Low Expression Associated with Poor Overall Survival

Introduction

Breast cancer (BC) is well known as a highly heterogeneous disease comprising distinct molecular subtypes.^1,2 Each BC subgroup differs from each other in terms of biological characteristics, risk factors, treatment responses, and patient survival outcomes.^1,2 This variability in BC behaviours and characteristics presents major clinical challenges and implications with regard to prognosis and BC management.^3,4 Therefore, much effort has been made to stratify heterogeneous BC subtypes, in order to increase our knowledge of the pathobiology of BC as well as to discover new treatments.^5–8 Typically, BC tumors have been divided into subtypes according to immunohistochemical (IHC) markers of estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2), with the combination of differential gene expression.⁹ Hence, there are three major subtypes of BC, and these subtypes are luminal (ER positive), triple-negative breast cancer (TNBC), and HER2 BC.⁹ Clinicians can determine what therapeutic options are most effective by using well-defined standard IHC markers to identify BC tumor subtype. However, care is now increasingly shifting towards looking at genetic differences to develop more personalized medicine approaches.¹⁰ The Prediction Analysis of Microarray 50 (PAM50) predictor is a multigene signature that stratified traditional pathological classification of BC (TNBC, luminal, HER2) into intrinsic molecular subtypes—luminal A, luminal B, HER2 enriched, normal-like and basal-like (which commonly is TNBC) BC—based on multigene differential expression.^11,12 With this move to precision medicine for cancer treatment, there has been an increasing need to identify new potential biomarkers across the different subtypes of BC.

Forkhead box (FOX) family members are highly conserved transcription factors that have a major role in a plethora of biological functions. These family members have a common DNA-binding domain (DBD) composed of 100 amino acids that is known as the forkhead box or “winged helix” domain.¹³ Many FOX transcription factors have been linked to tumorigenesis, where oncogenic and/or tumor-suppressive roles of FOX genes have been suggested in several cancers.^14,15 Members of the FOXA, FOXC, FOXM, FOXO, and FOXP subclasses of FOX proteins, in particular, have direct effects on the initiation, maintenance, progression, and drug resistance of cancers.¹⁶ For example, overexpression of FOXM1 accelerates the proliferation and progression of prostate cancers in mouse models.¹⁷ The widely studied FOXO proteins are key regulators of tumour suppression, as the simultaneous deletion of FOXO1, FOXO3, and FOXO4 alleles in somatic cells invokes thymic lymphomas and systemic haemangiomas in mouse models.¹⁸ Moreover, FOXF2, FOXM1, FOXA1, and FOXC1 genes were suggested to play critical roles in BC initiation, proliferation, migration, and invasion.^15,19–21 FOXQ1 is another member of the FOX protein superfamily of transcription factors that regulate the expression of genes necessary for embryonic development, cell proliferation, differentiation, and apoptosis.²² This single-exon gene is located on chromosomal region 6p25.3 and its encoded protein is characterized by a distinctive evolutionarily conserved DBD.²³ Functionally, FOXQ1 plays a role in a diverse range of important biological processes such as angiogenesis, epithelial differentiation, smooth muscle differentiation, mucin secretion, and natural killer cell effector function activation.^13,22,24,25 Moreover, FOXQ1 has been associated with the aggressive phenotype of cancers such as breast, colorectal, ovarian, and pancreatic cancers, since altered FOXQ1 expression has been previously detected in human tumor specimens.^26–31 Increasing numbers of studies have examined the role of FOXQ1 in tumor progression and it has been suggested that suppressing FOXQ1 expression may decrease invasion and migration in TNBC cell lines.^32,33 FOXQ1 has also been suggested as potentially being a driving force in the heterogeneity of BC.⁶ Thus, additional knowledge of the role of FOXQ1 in BC could have the potential to improve the diagnosis and treatment of BC tumors. In particular, the role and expression of FOXQ1 in luminal and HER2 BC needs to be determined.

In this study, we investigated FOXQ1 mRNA expression and copy number variation (CNV) across BC patient subtypes and cell lines. Furthermore, we assessed the prognostic significance of FOXQ1 for BC patients. Our results revealed for the first time that FOXQ1 mRNA is differentially expressed in TNBC, luminal, and HER2 BC patients and cell lines. We found that FOXQ1 expression is significantly lower in luminal and HER2-positive tumors and that low expression of FOXQ1 predicts poor overall survival in BC patients. We also found that FOXQ1 has more copies in tumors from TNBC patients compared to normal breast tissue samples. Interestingly, FOXQ1 mRNA expression is independent of its CNV in BC cell lines. Our findings highlight that low expression of FOXQ1 mRNA is associated with significantly poorer survival for different classes of BC.

Materials and Methods

Tissue Samples

BC patients’ samples were obtained with approval of the University of Alberta Health Research Ethics Board (Pro00018758) with written informed consent, in accordance with the Declaration of Helsinki. The tumor samples from BC patients were collected at surgery and frozen in liquid nitrogen within 20 min of devitalization for further experiments. Normal human breast tissue was obtained from breast reduction surgery. BC subtypes have been defined using immunohistochemistry (IHC) markers; ER, Estrogen receptor; PR, Progesterone; and HER2, Human Epidermal Receptor (Table 1).

Table 1 BC Patient Demographics

Overall Survival Analysis

Kaplan-Meier curves were generated using the online analysis tools OncoLnc (http://www.oncolnc.org)³⁴ and Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/).³⁵ OncoLnc and GEPIA use patient survival data and level 3 RNA sequencing expression data of BC tumor samples from the Cancer Genome Atlas (TCGA) invasive breast carcinoma project (BRCA) database (https://portal.gdc.cancer.gov/projects/TCGA-BRCA). Patient samples (n=1006) and (n=1055) from TCGA-BRCA were used to generate the Kaplan-Meier (KM) curves in OncoLnc and GEPIA online tools, respectively. Patient samples in OncoLnc and GEPIA were divided into high and low groups based on the median and the quartile of FOXQ1 expression in BRCA database. The median cut-off divided patients into two groups, high FOXQ1 expression group are above the 50th percentile, and low FOXQ1 expression group are below the 50th percentile. For the quartile cut-off analysis, patients were categorized into two groups, high FOXQ1 expression group are above the third quartile (higher than 75th percentile), and low FOXQ1 expression group are below the first quartile (lower than 25th percentile). More details on the mRNA expression data and their normalization can be found in (http://www.oncolnc.org).³⁴ The long-rank test was used for the hypothesis evaluation. In regard to the censorship method, we used right censoring method in which the patients who did not experience the event (death) during the study duration and had no available further follow-up data were censored. The mean age at diagnosis of the BC patients in TCGA-BRCA database is 58.4± 13.2 years and the median follow-up time is 27.7 months where the median follow-up time for the overall survival time to event is 41.8 months and the median follow-up time for the overall survival time to censor is 25.0 months.³⁶

Kaplan-Meier survival curves were generated of FOXQ1 expression in BC subtypes using UALCAN online tool (http://ualcan.path.uab.edu/index.html).³⁷ TCGA-BRCA patient samples were divided into two groups, high and low FOXQ1 expression groups. Patients with FOXQ1 transcript per million (TPM) above the third quartile (higher than 75th percentile) were categorized as high FOXQ1 expression group, while patients with FOXQ1 TPM below the third quartile (lower than 75th percentile) were categorized as low/medium FOXQ1 expression group. Subsequently, patients within each group were further stratified into TNBC, luminal, and HER2 BC subtypes. Log rank test was used to generate the p value to test statistical significance of survival correlation between BC subtypes and FOXQ1 expression. UCSC Xena website: https://xenabrowser.net/³⁸ was used to generate Kaplan-Meier curves for the TCGA-BRCA data with a cut-off time after 10 years. Log rank test was used to generate the p value to test statistical significance of survival correlation.

SPSS version 25 (IBM Corp. Released 2017. IBM SPSS Statistics for Windows, Version 25.0. Armonk, NY: IBM Corp.) was used to run the Cox’s regression hazard model. p-value <0.05 was used for statistical significance. We used the TCGA-BRCA data (https://portal.gdc.cancer.gov/projects/TCGA-BRCA) to conduct the univariate overall survival of Cox’s regression analysis for FOXQ1 expression. One thousand six TCGA-BRCA patients data with high and low FOXQ1 expression with respect to the median and quartile cut-off were downloaded from OncoLnc (http://www.oncolnc.org).³⁴ OncoLnc includes only patients who have all the necessary clinical data for each cancer such as age, sex, and grade as well as gene expression. More details on TCGA_BRCA patient data used for FOXQ1 Cox regression analysis can be found in Supplementary Table 1. Unadjusted and adjusted Cox’s regression hazard model was used to determine the association of FOXQ1, FOXF2, and FOXM1 with overall survival, more details can be found in Supplementary Table 2.

Finally, we also used TCGA-BRCA data (https://portal.gdc.cancer.gov/projects/TCGA-BRCA) to conduct Cox’s regression analysis of FOXQ1 expression adjusted to BC subtypes. One thousand six BC patient data with high and low FOXQ1 expression with respect to the median were downloaded. Then, we identified each BC patient subtype as TNBC, luminal, and HER2 as previously described in Ref. 39. More details on the latter analysis can be found in Supplementary Table 3.

RT-qPCR

RT-qPCR was conducted according to the detailed protocol as previously described.⁴⁰ Total RNA was extracted from BC patients’ tumor samples and human breast normal tissue using the RNAqueous kit (Ambion, Streetsville, ON, Canada) according to the manufacturer’s instruction. cDNA was synthesised from 1 µg RNA using Superscript II reverse transcriptase.⁴¹ RT-qPCR assays were performed with a QuantiTect SYBR Green PCR kit (Applied Biosystems Foster City, CA, United States) and analyzed on a LightCycler 96 Real-Time PCR System (Roche Life Science, Penzberg, Germany), at least three times with each reaction in triplicate. mRNA levels were normalized to TFRC (selected as an internal control gene) through the ΔΔCt method and changes in mRNA levels were described in fold change compared to the control samples. Primers for FOXQ1 and TFRC were designed using the Primer3 software (http://primer3.ut.ee/); FOXQ1: F: 5′-CGGAGATCAACGAGTACCTCA-3′, R: 5′-CAGTCGTTGAGCGAAAGGTT-3′; TFRC F: 5′-AACAACAGATTTCGGGAATGC-3′, R: 5′-CGTAGGGAGAGAGGAAGTGATA-3′. For statistical analyses, we first compared FOXQ1 expression in all four groups; normal breast tissue, TNBC, luminal, and HER2, using one-way ANOVA. Then, one-tailed unpaired Student’s t-tests and Bonferroni correction were used for multiple comparisons to assess statistical differences.

Copy Number Variation

FOXQ1 gene copy number in BC patients was quantified through Real-time qPCR, as previously described⁴² to determine if the altered dosage of FOXQ1 gene affected expression levels in BC tumors. Genomic DNA was extracted from human breast normal tissue and from TNBC, luminal, and HER2 BC patient tumor samples using the EZ‐10 Spin Column Genomic DNA Minipreps Kit (Bio Basic Inc., Markham, ON, Canada) according to the manufacturer’s instructions. Real-time qPCR was then conducted on the genomic DNA of all samples to measure FOXQ1 dosage. PCR reactions for normal breast tissue, TNBC, luminal, and HER2 BC patient tumor samples (Table 1) were analyzed in triplicate. Normal breast tissue samples were used as a control to normalize FOXQ1 dosage. We used the ΔΔCT method as our quantification strategy, with GJA5 selected as an internal control gene. Average CT values of triplicates were calculated for each sample. ΔCT for each sample was then calculated by subtracting the average CT number of FOXQ1 from that of GJA5. FOXQ1: F: 5′-CGGAGATCAACGAGTACCTCA-3′, R: 5′-CAGTCGTTGAGCGAAAGGTT-3′; GJA5 F: 5′-AGTTCCCAGCCAATAGACAGC-3′, R: 5′-AAGGCTGAGTAGAGGGAGGAG-3′. We used a two-tailed unpaired Student’s t-tests for comparisons of each BC subtype with normal breast tissues to assess statistical differences.

K-Means Clustering

K means clustering analysis was performed to investigate and compare FOXQ1 CNV and its mRNA expression. In order to conduct K-means clustering, we first identified the ideal number of clusters. Our initial approach was to compare the biological ideal number of clusters, which is 3 clusters, based on cell line annotation and classification; TNBC, luminal, and HER2 (supervised K-means clustering). For the unsupervised ideal number of clusters, we used the elbow method⁴³ to determine the ideal statistical number of clusters for K-means with the scikit-learn’s method of the MinMaxScaler module used with a sum of squares distance and fit between 1 and 15 possible clusters.⁴³ We then devised Python scripts to identify the cluster sets by applying standard supervised K-means clustering using both number of clusters (supervised and unsupervised) and FOXQ1 mRNA expression and copy number variation (CNV) as seeds as described.⁴⁴ The mRNA and CNV data of FOXQ1 were obtained from Cancer Cell Line Encyclopedia online (CCLE) database.⁴⁵ K-means clustering performed 1000 iterations using the Euclidean distance and using seeds between 1 and 3 clusters and 1 and 4 clusters. Thirty-six BC cell lines were used for the analysis. Cell lines were grouped as TNBC, luminal, and HER2 BC cell lines as described.^46,47 CCLE CNV data of BC cell lines were obtained using genome-wide human Affymetrix SNP Array 6.0.⁴⁵ The mRNA data on the CCLE website were produced by Affymetrix Human Genome U133 Plus 2.0 Array followed by Robust Multiarray Averaging (RMA) method of microarray normalization and were uploaded as log(2) gene expression signal.⁴⁵ The log(2) FOXQ1 expression data for each of the 36 BC cell lines from CCLE were obtained and then scaled to values between 0 and 1 using the equation;

Zi = {xi-min(x)}/{max(x)-min(x)}

where x = (x₁, … …,x_n) producing z_i as the scaled value for FOXQ1 mRNA expression in each of the 36 BC cell lines. CNV raw data created by Affymetrix (CEL file) were normalized to copy number estimates using a GenePattern pipeline.⁴⁵ 0.5 was used as a cut-off, where cell lines with values above the cut-off were considered as relative high expression of FOXQ1 cell lines, while cell lines with FOXQ1 expression below the cut-off were considered as relative low expression of FOXQ1 cell lines.

Breast Cancer Patient Publicly Available Dataset

RNAseq data for FOXQ1 expression in breast cancer patients from the TCGA-BRCA (https://portal.gdc.cancer.gov/projects/TCGA-BRCA) were obtained from UCSC Xena website: https://xenabrowser.net/³⁸ to generate Boxplots analysis across breast cancer subtypes. One-way ANOVA was applied to generate the p value to test statistical significance. Microarray data of FOXQ1 expression in breast cancer patients from the Van’t Veer⁴⁸ study were accessed using UCSC Xena website: https://xenabrowser.net/.³⁸ Welch’s t-test was applied to generate p values to test statistical significance.

FOXQ1 mRNA expression correlation with each of PDGFRα and β, Twist1, Zeb2, and DACH1 mRNA expression was analyzed in the TCGA-BRCA database using the co-expression analysis tool on cBioPortal.^96,97 Spearman and Pearson were used for correlation analysis. FOXQ1, PDGFRα and β, Twist1, Zeb2, and DACH1 log2 mRNA values of normalized RSEM RNAseq V2⁴⁹ in 994 patient and normal samples were obtained from the PanCanAtlas 2018 of the TCGA-BRCA https://gdc.cancer.gov/about-data/publications/pancanatlas.⁵⁰

Results

FOXQ1 mRNA Expression Varies Between BC Patient Subtypes

To evaluate FOXQ1 expression across the different subtypes of BC, RT-qPCR assays were conducted on TNBC, luminal, and HER2 BC patient samples and on normal breast tissue samples (Table 1). One-way ANOVA analysis showed that FOXQ1 mRNA expression varied between BC patients’ subtypes, F(30,52) = 3.8, p=0.026. Furthermore, we found that FOXQ1 expression was significantly lower in luminal BC when it was compared to normal breast tissue (Figure 1A). The mRNA level of FOXQ1 was lower in HER2 BC compared to normal breast tissue (p=0.02) (Figure 1A). Interestingly, FOXQ1 mRNA levels were lower in luminal and HER2 BC compared to TNBC (p=0.017 and p=0.029, respectively) (Figure 1A). We further confirmed our results in an expanded panel of BC patients (n = 522) using a publicly available dataset from the Cancer Genome Atlas (Figure 1A). Consistent with our RT-qPCR, we found that the expression of FOXQ1 mRNA varied across BC patient subtypes (Figure 1A, F = 24.64, p<0.0001). To substantiate our findings of significant low expression of FOXQ1 in luminal BC, we further investigated FOXQ1 expression in a public dataset of BC patients after they were stratified by ER status into ER positive and ER negative. Consistent with our RT-qPCR, we found that the expression of FOXQ1 was significantly lower in ER-positive BC patients compared to ER-negative BC patients (Figure 1B). These results indicate that FOXQ1 mRNA levels are differentially expressed among BC patients and that FOXQ1 could potentially have different roles based on the BC subtype.

Figure 1 FOXQ1 is differentially expressed across breast cancer subtypes and is under-expressed in luminal and HER2 breast cancer patient samples. (A) Left, qPCR experiments were conducted to measure FOXQ1 mRNA levels in normal breast tissue (Ctr, n=6) acquired from reduction mammoplasties, TNBC (n=6), luminal (n=6), and HER2 (n=6) BC patient samples. One-way ANOVA analysis was used to asses statistical difference in all groups *F(30,52) = 3.8, p=0.026 followed by unpaired t-tests and Bonferroni corrections were used for multiple comparisons. α= 0.05, adjusted α= 0.0125. The ΔΔCT method was used for analysis and TFRC was used as a reference gene. Error bars represent standard error of the mean (SEM). (A) Right, box plots showing FOXQ1 expression across BC intrinsic subtypes (PAM50). FOXQ1 RNAseq data in BC subtypes from the TCGA-BRCA data95 were used. Luminal B (orange), n = 127, Minimum = 0.00, Q1 = 2.0, Median = 3.29, Q3 = 4.40, Maximum = 6.26; Luminal A (red), n = 231, Minimum = 0.527, Q1 = 3.31, Median = 4.32, Q3 = 5.26, Maximum = 7.08; HER2-enriched (brown), n = 58, Minimum = 1.17, Q1 = 3.58, Median = 4.50, Q3 = 5.44, Maximum = 6.94; Normal-like (pink), n = 8, Minimum = 3.13, Q1 = 4.25, Median = 4.58, Q3 = 5.35, Maximum = 5.49; Basal-like (purple), n = 98, Minimum = 0.00, Q1 = 3.54, Median = 5.25, Q3 = 6.68, Maximum = 9.65. One-way ANOVA was used for statistical significance (F = 21.25, p<0.0001). FOXQ1 expression is expressed in RSEM (RNA-seq by expectation-Maximization) format. (B) Left, box plot of FOXQ1 RNAseq data in BC patients stratified by ER status. ER-positive (red), n = 601, Minimum = 0.00, Q1 = 2.82, Median = 4.12, Q3 = 5.05, Maximum = 8.13; ER-intermediate (orange), n = 2, Minimum = 4.24, Q1 = 4.24, Median = 4.50, Q3 = 4.50, Maximum = 4.50; ER-negative (blue), n = 179, Minimum = 0.00, Q1 = 3.58, Median = 4.68, Q3 = 6.26, Maximum = 9.65. One-way ANOVA was used for statistical significance (F = 24.64, p<0.0001). (B) Right, box plot of microarray data of FOXQ1 expression in BC patients stratified by ER status from the Van’t Veer et al48 study. ER-positive (red), n = 71, Minimum = −1.46, Q1 = −0.67, Median = −0.20, Q3 = 0.08, Maximum = 1.06; ER-negative (blue), n = 46, Minimum = −1.22, Q1 = −0.38, Median = 0.5, Q3 = 1.06, Maximum = 2.87. Welch’s t-test was used for statistical significance (***p<0.0001). Whiskers represent minimum and maximum values. The TCGA-BRCA95 and the Van’t Veer et al48 data were accessed using UCSC XENA cancer genome browser: https://xenabrowser.net/.

Low Expression of FOXQ1 is Associated with Poorer Overall Survival in Breast Cancer Patients

Kaplan-Meier curves were performed on the TCGA-BRCA dataset using the OncoLnc,³⁴ GEPIA,³⁵ and UALCAN³⁷ online tools to investigate whether FOXQ1 mRNA expression is prognostic for life expectancy. BC patients were divided into two risk groups with high and low FOXQ1 expression levels with respect to the quartile and the median used as the cut-off values. BC patients with low FOXQ1 mRNA expression had significantly shorter overall survival time compared to those with high FOXQ1 mRNA expression (Figure 2, p=0.00473 and p=0.0479) and (Supplementary Figure 1, p=0.0087 and p=0.042), and after 10 years follow-up (Supplementary Figure 2).

Figure 2 Low FOXQ1 expression predicts poor overall survival in BC patient. Kaplan-Meier (KM) analysis identified low and high-risk BC patient groups based upon significant differences of FOXQ1 mRNA in BRCA-TCGA database. BC patients were divided with high and low FOXQ1 expression levels with respect to the quartile and the median used as the cut-off value. This graph was generated using the bioinformatics online tool OncoLnc (http://www.oncolnc.org)34 and then modified to include the number of patients at risk at 0, 5 and 10 years.

In addition, a univariate Cox’s regression analysis was conducted for overall survival on the TCGA-BRCA patient data with respect to the expression of FOXQ1. The results of the univariate Cox’s regression analysis in Table 2 indicate that BC patients who have high FOXQ1 expression (≥50th percentile) were significantly associated with 29% lower risk of death as compared to BC patients who have low FOXQ1 expression (p=0.048). Further analysis was conducted, comparing the highest quartile and the lowest quartile of FOXQ1 expression. The results indicate that FOXQ1 high expression (≥75th percentile) was significantly associated with about 51% lower risk of BC patient’s death as compared to the low FOXQ1 expression (below 25th percentile) with a p-value of 0.006 (Table 2).

Table 2 Hazard Ratios from the Univariate Cox’s Regression Analysis for the TCGA-BRCA Database

Furthermore, using UALCAN online tool, we examined the correlation of FOXQ1 expression and BC subtypes by dividing patients into two groups; high and low/medium FOXQ1 expression cohorts (Figure 3). Kaplan-Meier curves indicate that patients with low/medium FOXQ1 expression had significantly shorter overall survival time than those with high expression of FOXQ1 (Figure 3, p=0.011). Subsequently, using the TCGA-BRCA data, we conducted Cox’s regression analysis for FOXQ1 expression stratified by BC subtype as well as for FOXQ1 expression adjusted for BC subtype (Table 3). For HER2 subtype, high FOXQ1 expression was significantly associated with about 73% lower risk of death as compared to low FOXQ1 expression (Table 3, p=0.024). The same trend was observed for luminal and TNBC subtypes indicating that high FOXQ1 expression showed a reduced risk of death by 19% and 32%, respectively. However, these associations were not statistically significant (Table 3).

Figure 3 Low FOXQ1 expression predicts poor overall survival in BC patient subtypes. Kaplan-Meier (KM) analysis of the effect of high and low/medium FOXQ1 expression on overall survival of HER2, luminal, and TNBC patients shows a cumulative significance of p=0.011. This graph was generated and modified using the bioinformatics online tool UALCAN online tool (http://ualcan.path.uab.edu/index.html)37.

Table 3 Cox’s Regression Analysis for FOXQ1 Expression Based on Median Cut-off Stratified by BC Subtypes and Adjusted for BC Subtypes

When adjusting for BC subtype, high FOXQ1 expression showed a 32% lower risk of death (Table 3, p=0.036). The unadjusted analysis from Table 2 and adjusted analysis from Table 3 indicate that FOXQ1 is an independent predictor of overall survival in BC patients. Together, these results demonstrate that low expression of FOXQ1 predicts poor overall survival in patients with BC.

Relation of FOXQ1 Copy Number in BC Subtypes to mRNA Levels

It has recently been suggested that CNV correlates with gene expression in BC cell lines and patient tissue.⁵¹ To determine FOXQ1 CNV in the subtypes of BC patients, qPCR assays were conducted on TNBC, luminal, and HER2 BC patient DNA samples (Table 1) and on normal breast tissue DNA samples. The number of copies of FOXQ1 was significantly higher in TNBC patients compared to normal breast tissue p=0.03 (Figure 4). Although FOXQ1 appeared to have a trend for more copies in HER2 BC patients, this did not reach statistical significance (Figure 4, p=0.08).

Figure 4 FOXQ1 has more copies in TNBC compared to control samples. qPCR experiments were conducted on genomic DNA to measure FOXQ1 dosage in TNBC, luminal, and HER2 BC patient tumor tissues. The ΔCT method was used for analysis and GJA5 was used as a reference gene. FOXQ1 dosage was calculated using 2 [ΔCT sample – ΔCT control]. Control (Ctr) samples are normal breast tissue samples acquired from reduction mammoplasties. Ctr (n=6), TNBC (n=6), luminal (n=6), and HER2 (n=6). Two-tailed unpaired Student’s t-tests were used for comparisons of each BC subtype with normal breast tissue to assess statistical differences. Error bars represent standard error of the mean (SEM). *p<0.05.

Unsupervised and supervised K-means clustering analyses (Table 4 and Figure 5) were performed to identify clusters of BC cell lines as well as to investigate any associations between FOXQ1 mRNA expression and gene copy number. Unsupervised K-means clustering identified 4 clusters among BC cell lines for FOXQ1 mRNA compared to their copy number (Table 4 and Figure 5A).

Table 4 BC Cell Line Names and Types That Were Used for Standard K-Means Clustering

Figure 5 FOXQ1 expression is independent of its CNV in BC cell lines. (A) unsupervised and (B) supervised K-means clustering analyses show different clusters of BC cell lines that have similar ranges of CNV but different FOXQ1 expression (the orange cluster vs the red cluster). Unsupervised clustering identified an extra cluster in BC cell lines compared to supervised clustering (the olive cluster, lower left).

We found some cell lines clusters that have similar copy number (orange cluster vs red cluster, Figure 5A) but different mRNA expression (Figure 5A). We also found similar results after we conducted supervised K-means clustering, where some cell line clusters had high and/or low FOXQ1 mRNA expression but similar ranges of variation in relative copy number (Figure 5B). For the supervised K-means clustering, we chose 3 as the number of clusters to correspond to the 3 different subtypes of BC; TNBC, luminal, and HER2 (see K-means clustering in method’s section). Interestingly, unsupervised and supervised K-means clustering analyses suggest that the steady state level of FOXQ1 mRNA expression appears to be independent of gene copy number in BC cell lines. Thus, FOXQ1 mRNA expression levels are not correlated with FOXQ1 CNV in BC subtypes. However, consistent with our results where FOXQ1 expression was lower in luminal and HER2 BC patient tumors (Figure 1), we found that FOXQ1 expression was low in all luminal and HER2 cell lines (Figure 5A and B). Similar to our observation of increased expression of FOXQ1 in TNBC patient tumors compared to luminal and HER2 (Figure 1), we found that FOXQ1 has elevated expression in some TNBC cell lines as compared to luminal and HER2 BC cell lines (Figure 5A and B). Interestingly, we found numerous TNBC cell lines cluster with other luminal and/or HER2 cell lines and have low FOXQ1 expression. This difference of FOXQ1 mRNA expression within the TNBC classification might be attributed to the distinct subgroups (basal A and basal B) that exists within TNBC cell lines.⁴⁷

FOXQ1 Expression is an Independent Predictor of Overall Survival in BC Patients When Adjusted to FOXF2 and FOXM1 Expression

We also investigated the overall survival of FOXQ1 adjusted for FOXF2 and FOXM1 expression. We evaluated FOXF2 and FOXM1 genes due to their critical role in BC initiation, proliferation, migration, and invasion.^19,20,52

TCGA-BRCA patients were divided into low and high expressing FOXF2, FOXM1, and FOXQ1 groups using the medians as cut-off values, followed by a univariate and multivariate Cox’s regression overall survival analysis. FOXM1 high expression was associated with 1.5 times higher risk of BC patient death as compared to FOXM1 low expression values (p=0.013; Table 5). Adjusted Cox’s regression model with FOXQ1, FOXM1, and FOXF2 expression levels indicated that FOXQ1 and FOXM1 were significantly associated with overall survival (Table 5; p=0.05 and p=0.011, respectively). The results indicate that FOXQ1 and FOXM1 were independent predictors of overall survival in BC patients. Having high expression of FOXQ1 was associated with 32% lower risk for overall survival when adjusted for FOXF2 and FOXM1 (p=0.05; Table 5). Higher expression of FOXM1 was associated with 1.6 times higher risk of death when adjusted for FOXQ1 and FOXF2 expressions, this association was statistically significant (p=0.011; Table 5). Our results show that FOXQ1 expression correlates inversely with overall survival in BC patients and that FOXQ1 expression is an independent predictor of BC overall survival when adjusted to the expression of FOXF2 and FOXM1 in BC patients.

Table 5 Cox’s Regression Analysis for FOXQ1, FOXF2, and FOXM1 Expression Levels from the TCGA-BRCA dataset Based on Median Cut-off

Discussion

FOXQ1 plays an important role in BC tumor re-initiation,⁵³ stemness and chemoresistance,⁵⁴ epithelial–mesenchymal transition (EMT), invasion, and metastasis.²⁸ Increasing evidence supports the hypothesis that FOXQ1 is a biomarker for cancer prognosis and diagnosis. Indeed, high levels of FOXQ1 expression can predict poor overall survival in hepatocellular carcinoma, gastric, colorectal, pancreatic, and non-small cell lung cancers.^29,55–57 However, the prognostic value of FOXQ1 in BC is not well known. In particular, the expression of FOXQ1 in other subtypes of BC, specifically luminal and HER2 BC, is not well studied. Thus, a molecular understanding of the mechanism of altered expression of FOXQ1 warrants further investigation. In this study, we report for the first time that FOXQ1 mRNA is differentially expressed across different types of BC patients and cell lines. In contrast to the earlier studies of a limited number of BC cell lines,^28,30,32,58 our studies of an expanded panel of thirty-six BC cell lines show that FOXQ1 has low expression in many TNBC cell lines and in all luminal, HER2 BC patients and cell lines. We also show that this low expression of FOXQ1 is associated with poor prognosis in BC subtypes and that FOXQ1 expression is an independent predictor of overall survival in BC patients.

In order to understand FOXQ1 expression across BC subtypes, we investigated FOXQ1 expression in BC tumor tissues (Table 1). FOXQ1 expression across BC subtypes was analyzed in TNBC, luminal, and HER2 tumor tissues (Figure 1A). Two previous papers^6,30 studied FOXQ1 expression data of BC tumors from the TCGA-BRCA database and found an overexpression of FOXQ1 in TNBC compared with luminal BC³⁰ and with HER2 BC.⁶ Similarly, it was suggested that FOXQ1 is overexpressed in BC tumor tissue compared with normal adjacent tissue;³² however, it was not clear what BC subtypes were used for this analysis. Consistent with these studies, we found that FOXQ1 mRNA levels were higher in TNBC tumor tissues when compared with luminal BC tumor tissues (p=0.017) (Figure 1A). As well, FOXQ1mRNA levels were higher in TNBC compared with HER2 BC (p=0.029) (Figure 1A). Importantly, we reveal for the first time that FOXQ1 is expressed at lower levels in luminal BC (p=0.008) and in HER2 BC (p=0.029) when each subtype was compared with unmatched normal breast tissue (Figure 1A). A limitation of our RT-qPCR analysis was the lack of available samples of matched normal breast tissue. However, our RT-qPCR analysis results were substantiated by our analyses of publicly available data from the Cancer Genome Atlas (Figure 1A and B). Interestingly, we found lower expression of FOXQ1 in luminal B than in luminal A (Figure 1A). Luminal B has been suggested to have higher levels of proliferation and cell cycle markers genes, but lower luminal markers compared to luminal A.^59–61 In addition, luminal A has a good prognosis compared to luminal B, where the latter has a higher recurrence rate, lower overall survival rate, and more invasive and proliferative cancers.^62,63 Whether FOXQ1 has a subtype-specific tumor-suppressive role in luminal BC, in particular in luminal B similar to its neighbour gene, FOXC1,⁶⁴ has yet to be determined. Nevertheless, our results reveal that significant differences in FOXQ1 expression occur in BC subtypes, with high FOXQ1 expression in TNBC and low FOXQ1 expression in luminal and HER2.

We further investigated the impact of FOXQ1 expression on BC patients’ overall survival. Using the OncoLnc,³⁴ and GEPIA³⁵ online tools, overall survival curves identified two risk groups with high and low FOXQ1 mRNA expression levels (Figure 2 and Supplementary Figures 1 and 2). We found that BC patients with low FOXQ1 expression had significantly shorter overall survival than BC patients with high FOXQ1 expression (Figure 2, and Supplementary Figures 1 and 2). Univariate Cox regression analysis for FOXQ1 expression, where BC patients who have high FOXQ1 expression were found to be significantly associated with lower risk of death as compared to BC patients who have low FOXQ1 expression (Table 2). Kaplan-Meier curves also showed that TNBC, luminal, and HER2 BC patients with low FOXQ1 expression have significantly poorer overall survival compared to TNBC, luminal, and HER2 BC patients with high FOXQ1 expression (Figure 3). When adjusting for BC subtype, high FOXQ1 expression showed significantly lower risk of death (Table 3) indicating that FOXQ1 is an independent predictor of overall survival in BC patients. This suggests that low expression of FOXQ1 in BC patients is significantly correlated with poor prognosis.

Importantly, we found FOXQ1 has lower expression in HER2 BC patient tissues (Figure 1A) and that HER2 BC patients with low FOXQ1 expression had significantly shorter overall survival compared to HER2 BC patients with high FOXQ1 expression (Table 3, Figure 3). This is an interesting result giving that HER2 BC, in comparison to other BC subtypes, is normally known to respond well to numerous targeted therapies that favourably impact patients overall survival.⁶⁵ On the other hand, TNBC lacks targeted therapy and is considered to have a clinically aggressive phenotype compared to other BC subtypes.⁶⁶ Although TNBC patient tissue had higher FOXQ1 mRNA expression compared to luminal and HER2 (Figure 1A), those TNBC patients with low expression of FOXQ1 appeared to have poorer overall survival compared to those with high FOXQ1 expression (Table 3 and Figure 3). This suggests that, clinically and case-by-case, BC patients in the TNBC cohort who have low FOXQ1 expression could have poorer clinical outcomes. Further studies that focus on stratifying the TNBC subtype further using gene signatures and cancer stages could reveal better understating of FOXQ1 expression in TNBC.

Intriguingly, while it has been reported that high FOXQ1 expression predicts poor overall survival in hepatocellular carcinoma, gastric, colorectal, pancreatic, and non-small cell lung cancers,^29,55–57 we found that high expression of FOXQ1 favourably impacts overall survival in BC (Tables 2 and 3), suggesting a dual roles of FOXQ1 across human cancers. Dual role of FOXQ1 and other FOX genes in cancers has been reported before. For instance, it was suggested that the ectopic expression of FOXQ1 promotes tumour growth in colorectal cancer by inhibiting apoptosis.²⁹ At the same time, FOXQ1 over-expression did not impact BC cell survival and/or proliferation in the highly metastatic 4T1 mouse cell line and it induced cell death in epithelial BC cell lines, MCF7 (labeled as luminal BC cell line) and BT-20 (labeled as TNBC cell line)²⁸ supporting a tumor-suppressive role of FOXQ1. High expression of FOXF2 enhanced EMT, migration and invasion of lung cancer cells⁶⁷ and in TNBC,⁶⁸ which directly contrasts other reports where low expression of FOXF2-induced EMT and was associated with poor overall survival in BC.¹⁹ Furthermore, high expression of FOXA1 was associated with poor overall survival in prostate cancer,⁶⁹ while high expression of the same protein favourably impacted BC prognosis.^70,71 In addition, FOXA1 activated luminal genes and repressed basal genes in ER-positive and ER-negative BC cells, respectively,⁷² indicating different roles of FOXA1 in different BC subtypes. This further highlights the complex and differing roles of FOX proteins in cancer. Our data also contrast the findings of Qiao et al who suggested that overall survival was significantly poorer in BC patients with high FOXQ1 expression.⁵⁸ Qiao et al⁵⁸ used the van de Vijver cohort to generate Kaplan-Meier plots of their overall survival analysis but the number of patients within that cohort who had high and low levels of FOXQ1 was not indicated.⁵⁸ Moreover, Van de Vijver et al⁷³ included microarray data of 295 BC tumors, some of which were lymph-node/ER-positive and some of which were lymph-node/ER-negative. OncoLnc,³⁴ GEPIA,³⁵ and UALCAN,³⁷ on the other hand, use up-to-date RNA sequencing data of the TCGA-BRCA database. In our study, we applied OncoLnc and GEPIA (respectively) to obtain Kaplan-Meier curves on BC tumors (n=1006 and n=1055), with high (n=503 and n=524), and low (n=503 and n=531) levels of FOXQ1 expression (Figure 2 and Supplementary Figure 1). The TCGA-BRCA data were also used to perform the adjusted and unadjusted Cox regression (see methods and Supplementary Tables 1, 2, and 3). The differences in numbers of BC tumor tissue samples, BC subtypes, and assays used could explain why we obtained different clinical outcomes of FOXQ1 expression compared to Qiao et al⁵⁸.

We also measured FOXQ1 CNV in BC samples (Table 1) to determine if altered copy number could underlie the differences in mRNA expression levels in TNBC, luminal, and HER2 tumors (Figure 4). Although there are numerous studies that have identified an association between copy number alteration and altered gene expression, how CNVs alter gene expression in BC is not well understood.^51,74,75 We found that FOXQ1 copy number is significantly amplified in TNBC compared with normal breast tissue (Figure 4). In contrast, no significant changes in FOXQ1 copy number were found in luminal and HER2 BC compared to control samples (Figure 4). To explore the possible ramifications of these findings further, we performed K-means clustering to investigate and compare FOXQ1 CNV and its mRNA expression. Our K-means clustering analyses took two approaches; unsupervised (Figure 5A) and supervised (Figure 5B). For the supervised K-means clustering, we determined the biological ideal number of clusters based on BC cell lines classification, TNBC, luminal, and HER2. For the unsupervised K-means, we used the elbow method⁴³ to identify the ideal statistical number of clusters. For both methods, the numbers of clusters and the FOXQ1 expression and CNV were used as seeds. We found that some clusters have high copy number with low mRNA expression, while other clusters have high mRNA expression with no CNV (Figure 5A and B). Intriguingly, this suggests that FOXQ1 expression is independent from FOXQ1 copy number in BC cell lines.

Several interesting patterns, however, emerge from our K-means clustering analysis. Our supervised analysis investigated if FOXQ1 expression and CNV correlate within a specific BC subtype more than others. Intriguingly, while it was expected that TNBC cell lines would cluster together based on FOXQ1 expression and CNV (Figure 5B), we found some TNBC cell lines also cluster with HER2 and with HER2 and luminal BC cell lines (red and black clusters, Figure 5B). Furthermore, our unsupervised analysis identified 4 clusters (Table 4 and Figure 5A), an extra cluster compared to the supervised method (Table 4 and Figure 5B). This suggests that FOXQ1 expression and CNV could potentially identify a sub-population of BC within TNBC, HER2, and luminal BC cell lines. This supports the findings of Yang et al, which suggest that FOXQ1 expression could drive the heterogeneity of BC subtypes.⁶ Consistent with our findings of low levels of FOXQ1 expression in luminal and HER2 BC patient tumors (Figure 1), we found FOXQ1 mRNA to be low in all luminal and HER2 BC cell lines (Figure 5A and B). Interestingly, while other studies reported FOXQ1 overexpression in TNBC cell lines,^28,32 our results show that FOXQ1 has low expression in many TNBC cell lines (Figure 5A and B), but FOXQ1 is overexpressed in some TNBC cell lines (Figure 5A and B). The difference in FOXQ1 overexpression across TNBC cell lines is striking. It has been previously suggested that, on the basis of molecular features and morphology and invasion potential, TNBC might be further subdivided into basal A and basal B types.⁴⁷ Basal A has a less differentiated and more epithelial-like morphology, whereas basal B has a more mesenchymal-like morphology which is more invasive.⁴⁷ Basal B has additionally been characterized as exhibiting more stem-cell-like characteristics (“stemness”), as cells of this subgroup have been found to possess the CD44⁺/CD24^−/low phenotype⁴⁷ normally associated with mammary cancer stem cells (CSCs).^76,77 Moreover, a role of FOXQ1 in stemness,^54,78 and EMT²⁸ has been reported. This could provide an explanation of why FOXQ1 is overexpressed in certain TNBC cell lines rather than other FOXQ1 TNBC cell lines (Figure 5A and B). Further classification of TNBC cell lines into basal A and B could provide a clearer picture of the different consequences of FOXQ1 expression between these two TNBC subtypes.

Very few direct targets of FOXQ1 are known. These genes may be positively regulated by FOXQ1 such as PDGFRα and β, Twist1 and Zeb2⁵⁴ or negatively regulated by FOXQ1 such as DACH1 gene.³⁰ Consistent with these studies, we found in preliminary studies of BC gene expression databases a significant positive correlation between FOXQ1 and PDGFRα or β, Twist1 and Zeb2 and a significant negative correlation between FOXQ1 and DACH1 (Supplementary Figure 3). While it has been suggested that FOXQ1 regulates these genes that play roles in BC stemness and EMT, there is no specific transcriptional signature that has been proposed for FOXQ1 in either development or cancer. The underlying mechanisms altering FOXQ1 expression in BC subtypes, as well remain elusive. This is not surprising given the lack of knowledge on the regulatory machinery controlling the expression of transcription factors such as FOXQ1. The actual number of genes regulated by FOXQ1 is unknown but could number in the thousands based upon similar analyses of the direct targets of the related FOXC1 transcription factor.⁷⁹ Thus, while a “transcriptional addiction in cancer” has been proposed,⁸⁰ it remains a challenge to precisely identify in cancer how and when transcription factors are themselves regulated and activated, and then to understand the target genes regulated by the activity of these transcription factors in normal and disease states.⁸¹ Epigenetics, gene function redundancy, tissue specificity, predicting enhancer-promoter regions, and the role of coactivators are a few among many other challenges to studying transcription factors in cancer.⁸¹ Future studies that reveal the regulation of FOXQ1 expression, together with the identification of direct transcriptional targets of FOXQ1, might reveal how FOXQ1 expression is controlled and the consequences of its dysregulation of this regulation in BC.

Our research, and that from other groups,^29,55–58 clearly demonstrate that FOXQ1 is another member of the FOX class of transcription factors with a key emerging role in cancer.^14,15 Previous studies have investigated the role of FOXQ1 expression, regulation, and function in relation to the expression of FOXF2 and FOXM1 genes, function and roles in BC and other diseases. Our findings that FOXQ1 is an independent predictor of overall survival in BC when adjusted to FOXM1 and FOXF2 expressions highlight the prognostic significance of FOXQ1 expression in BC patients (Table 5). FOXQ1 and FOXM1 expression levels have previously been shown to have prognostic value in colorectal cancer patients.⁵⁶ Moreover, a recent study suggested that FOXQ1 expression is negatively regulated by FOXF2 in BC cells.⁸² FOXQ1, FOXF2, and FOXC1 are all located within a 300 kb region of chromosome 6, physically linking in close proximity three FOX transcription factors with substantial roles in BC,^15,19,83,84 development,^{13,24,25,85–87} and in several diseases.^88–94 Our results suggest that FOXQ1 expression in BC patients could have significant prognostic value for survival from BC.

Conclusion

In our study, we show that FOXQ1 is differentially expressed in BC subtypes. Low expression of FOXQ1 predicts poor overall survival in BC and is an independent prognostic predictor. Moreover, FOXQ1 CNV alterations are not associated with its expression in BC. Understanding the consequences of FOXQ1 expression in BC subtypes could have clinical value.