Gene Expression Meta-Analysis of Potential Metastatic Breast Cancer Markers.

BACKGROUND: Breast cancer metastasis is a highly prevalent cause of death for European females. DNA microarray analysis has established that primary tumors, which remain localized, differ in gene expression from those that metastasize. Crossanalysis of these studies allow to revile the differences that may be used as predictive in the disease prognosis and therapy.
OBJECTIVE: The aim of the project was to validate suggested prognostic and therapeutic markers using meta-analysis of data on gene expression in metastatic and primary breast cancer tumors.
METHOD: Data on relative gene expression values from 12 studies on primary breast cancer and breast cancer metastasis were retrieved from Genevestigator (Nebion) database. The results of the data meta-analysis were compared with results of literature mining for suggested metastatic breast cancer markers and vectors and consistency of their reported differential expression.
RESULTS: Our analysis suggested that transcriptional expression of the COX2 gene is significantly downregulated in metastatic tissue compared to normal breast tissue, but is not downregulated in primary tumors compared with normal breast tissue and may be used as a differential marker in metastatic breast cancer diagnostics. RRM2 gene expression decreases in metastases when compared to primary breast cancer and could be suggested as a marker to trace breast cancer evolution. Our study also supports MMP1, VCAM1, FZD3, VEGFC, FOXM1 and MUC1 as breast cancer onset markers, as these genes demonstrate significant differential expression in breast neoplasms compared with normal breast tissue.
CONCLUSION: COX2 and RRM2 are suggested to be prominent markers for breast cancer metastasis. The crosstalk between upstream regulators of genes differentially expressed in primary breast tumors and metastasis also suggests pathways involving p53, ER1, ERB-B2, TNF and WNT, as the most promising regulators that may be considered for new complex drug therapeutic interventions in breast cancer metastatic progression. Copyright© Bentham Science Publishers; For any queries, please email at epub@benthamscience.org.

Introduction

The most common cause of cancer death in European females is breast cancer; the majority of breast cancer mortalities arise from metastatic spread of the cancer and not from the cancer at the primary site [1, 2]. DNA microarray studies have established that tumors which remain localized and primary tumors which metastasize differ in gene expression [3]. Hence gene-expression profiling can aid the identification of prognostic markers for patient outcomes and could potentially aid the identification of new therapeutic targets. Accurately predicting the risk of metastasis could also enrich patient quality of life, sparing those without metastatic disease from aggressive chemotherapy which invokes a range of short and long term side effects [4].

A key application of gene expression profile analysis is to identify a small number of gene signatures that classify cancer phenotypes in relation to their prognosis [5, 6]. Genes can be linked to mediating breast cancer metastasis to different areas of the body and to specific stages of metastasis. Genes such as IL3RA2, VCAM1 and MMP2 are currently associated with mediating aggressive breast cancer metastasis to the lung [7], whereas ST6GLANAC5 is a specific mediator of breast cancer metastasis to the brain [8] and cytokeratin-19 has been found as a putative marker of stem cells from the breast [9]. Tumor suppressors such as Nm23, KAI1 and BRMS1 are linked to the prevention of detachment of tumor cells from primary tumors, while KISS1 and MKK4 are linked to reduced growth at secondary sites [10]. A subgroup of metastasis-inducing proteins (MIPs), namely secreted phosphoprotein 1 (SPP1) [11], S100 calcium binding proteins A4 (S100A4) [12, 13] and S100P [14] and anterior gradient 2 (AGR2) [15] have also been shown to be overexpressed in patients exhibiting sporadic and metastatic breast cancer and are associated with reduced patient survival [14, 16-19]. Hence the up-regulation and down-regulation of identified gene signatures can be assessed in terms of prognosis.

It is a major target in breast cancer research to determine the genetic mechanisms that underpin the metastatic processes, which include: tumor cell intravasation, cell survival during circulation, extravasation into new tissues and successfully inhibited growth at a secondary site [20]. The initial step facilitating the transmission of tumors from primary to secondary sites is thought to involve epithelial-mesenchymal transmission (EMT); a process that allows epithelial cells to exhibit enhanced motility and invasion [21-24]. Many genes influence EMT, overexpression of the FOXM1 gene has been demonstrated to stimulate EMT-like changes in cells and EMT has also been shown to increase the regulation of MMP1 [25, 26]. The next step, intravasation, involves cancer cells forming circulating tumor cells (CTC’s) in the blood or lymphatics and cancer dissemination occurs either via angiogenesis or lymphangiogenesis [27-29]. The final steps involve extravasation to different parts of the body, such as the lung, bone and liver and adaptation and proliferation within new tissues [30, 31]. Cancer cells which have escaped from the primary lesion prior to its removal at surgery can remain dormant for 10 years or more before manifestation [32].

Meta-analysis of the existing gene expression data aims to facilitate the identification of a range of prognostic biomarkers for metastatic breast cancer within this project. It is proposed that a fundamental reason why genes within prognostic signatures are so unreliable is due to the difficulty in distinguishing their role as ‘passenger’ or ‘driver’ genes within the metastatic phenotype [33] and that expression of passengers may vary greatly across different cases due to many factors, such as noise in transcriptional regulation, which consequently results in differential expression [34]. Constant and reliable differential expression though may indicate a role of a gene as a driver, and the defined marker gene candidates should be experimentally validated. There is currently an abundance of online databases holding microarray and RNAseq data, which can be used in a context-query driven manner and can provide information on gene expression, giving information on up-regulation and down-regulation of genes within different tissue types and perturbations. In this study we used a number of applications for meta-analysis of gene expression, alongside interactions data to validate a general significance of already-suggested metastatic markers. We have analyzed several defined datasets representing different contrasts of gene expression in metastatic breast cancer compared to non-metastatic breast cancer and normal tissue, to define metastasis-specific cell functions distinguishing metastasis from the non-metastatic breast neoplasms.

Methods

Gene list compilation: 52 genes were compiled from literature, including genes up-regulated in metastasis, genes down-regulated in metastasis and genes that mediate metastasis to specific sites such as the lung and the bone. The combinations of key words ‘breast metastasis & gene’, ‘gene expression’, ‘microarray’, ‘PCR’, ‘up-regulated’ and ‘down-regulated’ were used to search for relevant publications in the NCBI database. 52 genes were selected as they had been previously reported as differentially expressed in metastasic tumors when compared to primary tumors.

Meta-analysis of Gene Expression data was performed using the Genevestigator Biomedical V4 microarray database (https://genevestigator.com/gv/biomed.jsp) [34, 35]. The software was used to identify genes that are specifically up-or down-regulated in response to a set of perturbations [34]. The conditions tool was used to identify gene sets that were differentially expressed across different breast cancer neoplasms. The Genevestigator Similarity Search tool provides co-expressed gene relationships calculated from array data [35]. More information on how the Genevestigator software can be operated is described in extensive detail in the Genevestigator user manual (https://genevestigator.com/userdocs/manual/). The top 10 significant, differentially expressed genes from the literature list associated with breast metastasis were selected for further analysis. A number of metastatic tissue datasets and corresponding independent gene expression experiments were limited to 3 compared to 9 of independent primary breast cancer gene expression analysis datasets. Moreover, according to the data source notation, lymph nodes were used as sources of the metastatic tissue, though in the case of primary breast tumors the samples were taken directly from breast. It was the only available option to compare primary and metastatic breast cancer.

The Ingenuity Pathway Analysis (IPA) (Ingenuity Systems, Inc., Redwood City, CA, USA) (www.qiagen.com/ingenuity) [36] tool was used to identify the regulators upstream of the gene set and to reconstruct the relational functional networks. A set of gene IDs were entered into the software for a Core Analysis. The tool produces graphical networks from the uploaded gene dataset and the molecules associated with this dataset in IPA Knowledge base, the Upstream Regulator Analysis can be used to suggest transcriptional factors ranked by the z-scores reflecting the number of associated genes in the dataset.

Results

Functional Analysis of Literature-Derived List of Genes Associated with Breast Cancer Metastasis

Genes were identified as relevant to metastasis via literature mining and analyzed using the Genevestigator software. The perturbations tool was used to visualize the expression of genes over a variety of conditions including the expression of genes across different breast cancer types (Genevestigator database). 10 top differentially expressed genes which were consistently reported in the literature as up or down regulated in metastatic cancer progression (Table ) were selected for further analysis. Although all of these genes showed slight differential expression across different breast cancer samples, p-values for each gene within a set of different perturbations were not significant. Insignificance could occur due to the way samples were collected; metastasis samples were extracted from the lymph nodes of patients whilst non-metastatic samples were taken from primary tumors of breast ductal carcinomas with no sign of distant metastases, inconsistencies could also arise due to a limited number of metastatic samples [37]. An extensive literature search was performed to validate a potential functional input of selected genes with detected insignificant p-values (Table ).

Instances where Genevestigator results were significant (p<0.05) were automatically indicated by (✓). Instances where p >0.05 were consistently confirmed by literature (✓), or have consistently contradicted to published data (X). Consistency here, by a rule for an initial selection of these genes from a literature gene set, means that an opposite vector of differential expression for a gene hasn’t been reported.

Due to the fact that metastatic samples used for the original microarray experiments were derived from lymph nodes and that contamination by immune cells can occur, we could not rely on the promising marker contrasts where increment in gene differential expression being characteristic only for metastatic neoplasms (i.e. in the case of VCAM1). From our analysis, VCAM1 appears to be the best potential marker for metastatic breast cancer and hence should be validated via analysis of the gene expression in metastatic tissue samples, where contamination with immune cells has been avoided. We can also dismiss gene expression changes that occurred only in primary neoplasms when contrasted against normal tissues, as these represent markers for primary breast cancer. The FZD3 gene however, may be an exception, its expression is thought to be insignificant but is still high in metastatic tumors when compared to primary tumors and this tendency is supported by literature (Table ). Gene expression contrasts observed between metastatic tissues and primary breast tumors (DEPDC1, NUSAP1, FOXM1, MUC1) demonstrate good prognostic markers for metastatic cancer development.

Our analysis suggested that the expression of the COX2 gene is significantly downregulated in metastatic tissue compared to both primary tumors and normal tissue, but not in primary tumors compared with normal breast tissue and may be used as a differential marker in metastatic cancer diagnostics.

RRM2 expression decreases in metastatic breast cancer progression and can be suggested as a marker to monitor breast cancer signatures. Our study also supports MMP1, VCAM1, FZD3, VEGFC, FOXM1 and MUC1 as breast cancer onset markers as these genes demonstrate significant differential expression in breast neoplasms when compared with normal breast tissue (Table ).

It is important to note that differential gene expression vectors across breast cancer neoplasms were sometimes contradictory when comparing the results of the Genevestigator analysis and sources of literature. A factor which could largely contribute to discrepancies in the results are the sample types used within different studies. For metastatic samples, the Genevestigator analysis used samples solely from the lymph nodes, whereas in other instances, samples were taken from surgically resected breast tumors [36] while some studies used a combination of freshly frozen human breast tumor samples and transgenic mice models [40]. Differences in cell lines used for each study may contribute to inconsistencies in results.

Extending the Gene List by Co-Expressed Genes

The list of genes strongly and differentially expressed in metastasis (Table ) was used as a seed to identify other relevant genes potentially associated with metastatic breast cancer.

Thus, a co-expression analysis was performed for each of the top 10 differentially expressed genes listed in Table . A co-expression value was calculated as a Pearson’s correlation coefficient and co-expressed genes which possessed a Pearson’s correlation coefficient of >0.8 were included in further analysis (Table ). The aim of identifying other genes, from microarray data, that exhibit strong co-differential expression with members of our initial gene list (Table ) was to facilitate prediction of the upstream regulators associated with metastatic process.

Functional Analysis of the Metastasis-Specific Gene Set Extended by their Co-Expression Partners

The merged list of genes retrieved from the literature and co-expression analysis (Table & Table ) was uploaded in IPA for functional interrogation. Table presents the top-ranked results of analysis of upstream regulators based on their activation z scores calculated by IPA predictive algorithm. The list of potential regulators may be biased by our initial selection of genes. We however discuss a relevance, based on literature analysis, of predicted factors to metastatic breast cancer in our Discussion section.

Fig. ( reconstructs its predicted regulatory interactions and shows regulatory interactions potentially leading to activation of a positive metastatic breast cancer prognostic marker, VCAM1, and two positive prognostic markers for breast cancer, which do not show specific links to metastatic transition in our analysis, MMP1 and MUC1. Interestingly, only RRM2, which shows decreased expression in metastasis in comparison to primary breast cancer, is downregulated by TNF and p53, according to results based on IPA knowledge-base data. However, TNF up-regulates VCAM1. ERB2 up-regulates primary breast cancer markers and rather pro-metastatic TNF and p53 genes. The oestrogen receptor may also modulate expression of p53. The core predicted regulators form an interacting circuit that may regulate primary breast cancer markers and lead to differential expression of genes involved in the metastatic transition.

Fig. (1)

Functional interactions between top breast cancer metastasis markers (bright blue, red) and primary breast cancer prognostic markers (pastel blue red) with the IPA-predicted upstream regulators. Highlighted in purple-regulatory circuit that may lead to development of metastatic transition. Arrows indicate activation, blocked lines-inhibition. Bright line colours correspond to regulatory interactions that potentially lead to differential expression of the metastatic markers. Dashed lines indicate predicted direct interactions between proteins, depicting protein-protein binding and regulatory actions. Validated transcriptional regulation is shown by solid lines.

Discussion

Gene expression profiling is widely used to define potential disease biomarkers, however individual experiments and small cohorts may lead to predictions that are not sustained on a larger scale. An integration of available data in meta-analysis gives an opportunity to validate the consistency of relationships between gene expression and physiological conditions, and our analysis aimed for the validation of such relationships between predicted gene markers and metastatic breast cancer.

The availability of proper datasets for metastatic cancer tissue limits our understanding of the specificity of the detected associations. We assume, for instance, that there is a high degree of probability that immune cells could contaminate the metastatic tissue derived from lymph nodes. This restricts our ability to suggest genes, which are highly prevalent in metastasis when compared to primary breast cancer. Thus, we would rather focus on genes which are strongly downregulated in metastatic tissue. However, even the sustainable predictions supported by the literature, need to be validated in a clearly designed experimental setting. Our analysis suggests that expression of COX2 is significantly downregulated in metastatic tissues, but not in primary tumors when compared with normal breast tissue and expression of RRM2 also decreases during the evolution of primary breast cancer into the metastatic one. Therefore, these genes may be suggested for disease monitoring and diagnostics. Our study supports MMP1, VCAM1, FZD3, VEGFC, FOXM1 and MUC1 as breast cancer onset markers, as these genes have a significant differential expression in breast cancer neoplasms, compared with normal breast tissue (based on data from 9 microarray experiments/Genevestigator), however it does not show confident evidence of their differential expression in metastatic tissue compared to primary tumors (based on data from 3 microarray experiments/ Genevestigator).

Several potential regulators and functional pathways of the differential gene expression in breast metastasis were suggested by our analysis. CDKN1A [66-68] is identified as the top-ranked regulator of MMP1 and FOXM1 (Table ), two genes where up-regulation is associated with breast cancer metastasis and primary tumors. CDKN1A may act downstream of c-ERBB-2 [68], and c-ERBB-2 also increases the expression of MMP1 [69, 70]. In agreement with the results of our upstream analysis performed via IPA, c-ERBB-2 was shown to regulate expression of VEGFC and MUC1 [71, 72] and associated with poor prognosis in breast cancer [73].

TNF alpha is predicted to up-regulate the expression of several top breast cancer and metastasis markers (Table ). TNF has indeed been shown to increase the expression of MMP1 via up-regulation of IL4 [74] and to induce the upregulation of MUC1 [75] and VEGFC [76]. Interestingly, RRM2 is also on this list of TNF-regulated genes, as well as on the list of genes regulated by p53 function (Table ). The protein product of the RRM2 gene is known to inhibit WNT signaling [77] and may have a direct role in WNT signaling regulation during the metastatic process. One of the components of WNT signaling, The Frizzled protein, FZD3, is up-regulated in breast cancer metastasis [78] via the WNT pathway and was also defined in our analysis as a potential prognostic marker. Crosstalk between TNF and WNT signaling can be suggested for future experimental analysis as one of the promising interactions driving or strongly reflecting the metastatic process. The multi-step process of tumor development and expansion driven by accumulation of new mutations may explain a transition from TNF-dependent primary tumor to less proliferating metastatic cells due to altered p53 function and dysregulated LGALS3 and CDKN1A downstream pathways (Fig. ).

The product of the Lectin Galactoside-Binding soluble 3 (LGALS3) gene is involved in the regulation of MUC1 and in up-regulation of VCAM1 [79, 80]; it also activates VEGFC expression, which is linked to tumor angiogenesis [81]. Expression of VEGF factors and VCAM1 is also controlled by NFKB [81, 82] stimulating leukocyte recruitment and leukocyte infiltration of a tumor [83-85].

The Forkhead box M1 (FOXM1) gene was linked to the up-regulation of several genes within the gene set (Table ). It was also shown to increase the expression of estrogen receptor alpha (ER1) [86] predicted to be upstream of a number of genes within the gene set. Estrogen has been shown to increase the expression of MMP1 [87], but also to suppress the expression of VEGFC [88, 89], that correlates with the vector of differential expression of these genes in breast cancer contrasted to normal breast tissue and metastatic tissue contrasted to primary breast tumor (Table ).

A gene expression profiling study similar to ours analyzed microarray data from 189 invasive breast carcinomas combined with data from three published gene expression datasets. The study identified genes in parallel with this study such as FOXM1 and BIRC5 as genes that are overexpressed in grade 3 tumors providing strength for our findings [90].

Conclusion

The number of genes that were suggested in publications as markers for breast cancer metastasis (Table ) did not show consistent differential expression in association with metastatic conditions, their differential expression vectors mentioned in the literature were also contradictory. Without dismissing the potential biomarker role of the genes that did not show confident differential signal in our meta-analysis, we advocate a re-evaluation of their differential association with metastatic breast cancer and the usefulness of their expression level measurement in diagnostics. On the other hand, further analysis and validation is required for the gene-markers and interactions suggested in this study. Crosstalk between the regulatory pathways and transcription factors (Fig. ) is likely to define a complex landscape for realization of the metastatic cell program in different gene and environmental contexts and should be taken into account for complex therapies and a personalized medicine approach.

Table 1

Genes Identified by literature mining and their corresponding mRNA expression.

Study	MMP1	VCAM1	FZD3	VEGFC	COX2	DEPDC1	NUSAP1	RRM2	FOXM1	MUC1
Metastatic vs. normal breast tissue	+log2=0.15fold=1.09pval=0.620	+log2=0.51fold=1.59pval=0.348	+log2=2.71fold=5.9pval=0.006	-log2=-2.33fold=-0.09pval=0.012	-log2=-1.06fold=-2.09pval=<0.001	+log2=0.24fold=1.18pval=0.146	+log2=0.43fold=1.42pval=0.237	+log2=1.32fold=1.84pval=0.398	+log2=0.13fold=1.00pval=0.806	+log2=2.50fold=5.78pval=0.059
Metastatic vs. primary breast cancer	-log2=-2.73fold=-14.61pval=0.098	+log2=1.72fold=2.93pval=0.048	+log2=0.41fold=1.11pval=0.562	-log2=-0.09fold=-1.13pval=0.901	-log2=-1.04fold=-2.08pval=<0.001	-log2=-0.53fold=-1.66pval=0.332	-log2=-1.02fold=-3.22pval=0.368	-log2=-1.84fold=-3.36pval=0.038	-log2=-2.35fold=-9.59pval=0.076	-log2=-0.68fold=-2.47pval=0.617
Primary breast cancer vs. normal breast tissue	+log2=2.87fold=15.91pval=0.017	-log2=-1.21fold=-1.84pval=0.039	+log2=2.30fold=5.36pval=<0.001	-log2=-2.24fold=-4.04pval=<0.001	-log2=-0.02fold=-1.00pval=0.892	+log2=0.77fold=1.96pval=0.056	+log2=1.45fold=4.56pval=0.076	+log2=3.16fold=6.18pval=<0.001	+log2=2.48fold=9.63pval=0.012	+log2=3.18fold=14.25pval=0.004

Log(2)- ratios, fold change values and p values. Log(2)-ratios (log2), fold change values (fold) and p-values (pval) from Genevestigator analysis are shown. Up-regulation of genes is indicated by (+), down-regulation of genes is indicated by (-).

Table 2

Validated Genevestigator results.

	MMP1	VCAM1	FZD3	VEGFC	COX2	DEPDC1	NUSAP1	RRM2	FOXM1	MUC1
Metastaticvs. normal breast tissue	✓ + [38]	✓ + [39]	SIGNIFICANT✓+	SIGNIFICANT✓-	SIGNIFICANT✓-	✓ + [40]	✓ + [40]	✓ + [41] [42]	✓ + [43]	✓ + [40]
Metastatic vs. primary breast cancer	✗ +[44 [7][45 [46]	SIGNIFICANT✓+	✓ + [47]	✗ + [48] [49] [50] [51] [52]	SIGNIFICANT✓-	✗ + [53] [40]	✗ + [54] [40] [55]	SIGNIFICANT✓-	✗ + [43] [56]	✗ + [57][58
Primary breast vs. normal breast tissue	SIGNIFICANT✓+	SIGNIFICANT✓-	SIGNIFICANT✓+	SIGNIFICANT✓-	✗ + [59] [60] [61] [62]	✓ + [63] [40]	✓ + [40] [64][65	SIGNIFICANT✓+	SIGNIFICANT✓+	SIGNIFICANT✓+

Table 3

Genes co-expressed with the genes in the literature-derived dataset (Table 2) obtained by the Genevestigator co-expression analysis.

Gene Co-Expressed with a Gene in the Dataset	Gene (Literature Analysis Dataset)	Pearson’s Correlation Coefficient
GPX8	VEGFC	0.85
FST	VEGFC	0.84
LOX	VEGFC	0.84
PXDN	VEGFC	0.84
EHD2	VEGFC	0.82
HNRNPM	COX2	0.93
ATP6	COX2	0.93
DCAF6	COX2	0.91
ND2	COX2	0.87
ND3	COX2	0.86
GTSE1	DEPDC1	0.82
HJURP	DEPDC1 & FOXM1	0.8 & 0.92
KIF2C	DEPDC1	0.8
MKI67	NUSAP1	0.77
TPX2	FOXM1	0.9
DLGAP5	FOXM1	0.9
AURKB	NUSAP1	0.82
CCNA2	NUSAP1	0.81
UBE2C	NUSAP1	0.79
KIF4A	FOXM1	0.9
ELF3	MUC1	0.86
AGR2	MUC1	0.83
PIGR	MUC1	0.83 & 0.84
TMC4	MUC1	0.82 & 0.88
RASEF	MUC1	0.8 & 0.84
TMC5	AGR2	0.89
SLC44A4	AGR2	0.88
KRT19	MUC1	0.79
TSPAN1	AGR2	0.87
C9orf152	AGR2	0.89
ST6GALNAC1	AGR2	0.88
LOC100505989	AGR2	0.86
KIAA0101	RRM2	0.93
TOP2A	RRM2 & NUSAP1	0.93 & 0.82
ZWINT	RRM2	0.93
DTL	RRM2	0.93
CCNB2	RRM2 & FOXM1	0.92 & 0.88
DLGAP5	FOXM1	0.9
Gene Co-Expressed with a Gene in the Dataset	Gene (Literature Analysis Dataset)	Pearson’s Correlation Coefficient
HMMR	DEPDC1	0.84
MELK	RRM2	0.91
BIRC5	DEDPC1 & FOXM1 & NUSAP1	0.79 & 0.91 & 0.78
ASPM	FOXM1	0.87
NUF2	RRM2	0.89

Pearson’s Correlation Coefficient is shown for each pair of gene profiles.

Table 4

IPA Upstream regulators associated with the genes selected from literature analysis and co-expressed genes.

Upstream Regulator	Molecule Type	P-Value	Target Molecules in Integrated Dataset
CDKN1A	Kinase	6.31E-09	BIRC5,CCNA2,FOXM1,KIAA0101,KRT18,MMP1,TOP2A (7)
AR	ligand-dependent nuclear receptor	7.80E-09	BIRC5,CCNA2,KIF2C,MUC1,NUSAP1,PGR,PIGR,TPX2,UBE2C,ZWINT (10)
ERBB2	Kinase	1.38E-08	ASPM,BIRC5,CCNA2,CCNB2,MKI67,MMP1,MUC1,TOP2A,VEGFC,ZWINT (10)
FOXO1	transcription regulator	3.86E-08	ASPM,BIRC5,CCNB2,DEPDC1,DLGAP5,MMP1,NUSAP1 (7)
TNF	Cytokine	2.30E-07	ADIPOQ,BIRC5,ELF3,FST,HLA-DRB4,MMP1,MT-CO2,MUC1,PIGR,PPP2R1B,RRM2,VCAM1,VEGFC (13)
FOXM1	transcription regulator	1.70E-06	AURKB,BIRC5,CCNA2,FOXM1,GTSE1 (5)
estrogen receptor	Group	5.25E-06	KRT18,KRT19,MMP1,MT-CO2,MUC1,PGR,PRLR,VEGFC (8)
ESR1	ligand-dependent nuclear receptor	7.09E-07	ASPM,BIRC5,CCNA2,FST,KRT19,MKI67,MMP1,PDZK1,PGR,PRLR,RAMP3 (11)
ESR2	ligand-dependent nuclear receptor	5.67E-04	BIRC5,MMP1,PGR (3)
LGALS3	Other	1.70E-06	CCNB2,KRT18,KRT19,MUC1,VCAM1 (5)
TP53	transcription regulator	3.06E-04	AURKB,BIRC5,CCNA2,HMMR, KIAA0101,RRM2,TOP2A,TPX2,UBE2C (9)

P-value measures whether there is a statistically significant overlap between dataset genes and the downstream targets of a regulator. The number of different proteins that each upstream regulator has been reported to be associated with is indicated in brackets. The top genes differentially expressed in metastasis (Table ) are highlighted in red. In bold are TFs that may be responsible for regulation of suggested prognostic marker, RRM2.