Abstract
Changes in iron regulation characterize the malignant state. However, the pathways that effect these changes and their specific impact on prognosis remain poorly understood. We capitalized on publicly available microarray datasets comprising 674 breast cancer cases to systematically investigate how expression of genes related to iron metabolism is linked to breast cancer prognosis. Of 61 genes involved in iron regulation, 49% were statistically significantly associated with distant metastasis-free survival. Cases were divided into test and training cohorts, and the supervised principal component method was used to stratify cases into risk groups. Optimal risk stratification was achieved with a model comprising 16 genes, which we term the iron regulatory gene signature (IRGS). Multivariable analysis revealed that the IRGS contributes information not captured by conventional prognostic indicators (HR = 1.61; 95% confidence interval: 1.16–2.24; P = 0.004). The IRGS successfully stratified homogeneously treated patients, including ER+ patients treated with tamoxifen monotherapy, both with (P = 0.006) and without (P = 0.03) lymph node metastases. To test whether multiple pathways were embedded within the IRGS, we evaluated the performance of two gene dyads with known roles in iron biology in ER+ patients treated with tamoxifen monotherapy (n = 371). For both dyads, gene combinations that minimized intracellular iron content [anti-import: TFRCLow/HFEHigh; or pro-export: SLC40A1 (ferroportin)High/HAMPLow] were associated with favorable prognosis (P < 0.005). Although the clinical utility of the IRGS will require further evaluation, its ability to both identify high-risk patients within traditionally low-risk groups and low-risk patients within high-risk groups has the potential to affect therapeutic decision making. Cancer Res; 71(21); 6728–37. ©2011 AACR.
Introduction
A long history of experimental and clinical literature attests to the importance of iron in tumor cell growth (1). For example, high levels of dietary iron have been linked epidemiologically to the increased development of tumors in humans (2, 3). In animal models, high dietary iron increased the rate of tumor formation (4, 5); conversely, a low iron diet led to less rapid growth of tumor xenografts (6, 7) and reduced spontaneous mammary tumors in rats (7). This has led to the exploration of agents that restrict iron availability, such as iron chelators and antitransferrin receptor antibodies, as antitumor agents (8, 9).
Despite the wealth of descriptive data linking iron and cancer, pathways that underlie the apparent demand for iron by tumors remain largely uncharacterized. We recently uncovered a role for a pathway that mediates iron efflux in breast cancer growth and metastasis (10). This pathway is mediated by ferroportin and hepcidin. Ferroportin is an iron efflux pump and hepcidin is a peptide hormone that binds to ferroportin and triggers its degradation (11). We observed that ferroportin expression is reduced in breast cancer cells relative to normal mammary epithelial cells. Low ferroportin expression correlated with high levels of metabolically available iron and increased growth of tumor xenografts. Strikingly, the expression of ferroportin in concert with the expression of hepcidin independently predicted for metastasis-free survival of women after definitive primary treatment of their breast cancer in multiple independent breast cancer cohorts (10).
Although our previous observations linked genes mediating iron efflux to breast cancer outcomes, we did not attempt to identify an optimal iron homeostasis gene signature related to breast cancer prognosis. Rather, we focused on providing an explanation for the experimental observation that ferroportin and hepcidin were altered in breast cancer. In this article, we cast a wider net, using a formal bioinformatic analysis to assess the predictive value of all genes with readily identifiable roles in iron metabolism. The goal was to both identify an optimal iron gene signature and to discover whether one or more iron regulatory pathways were embedded in this signature. Here, we identify an iron regulatory gene signature (IRGS) that predicts outcome in breast cancer patients. Simple combinations of genes within this signature reveal different pathways of iron regulation that converge on a similar breast cancer phenotype.
Materials and Methods
Breast cancer datasets and microarray data processing
Expression profiles of iron-associated genes were analyzed in breast tumors using a multi-institutional “microarray meta-cohort” (MMC) comprising 6 distinct breast cancer cohorts (17–20) totaling 759 primary breast cancer cases. Datasets and associated clinical annotations were downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/). All tumor samples were analyzed from surgically resected frozen tissue and profiled on both the Affymetrix U133A and U133B series microarrays, or the U133 PLUS 2.0 GeneChip, according to standard Affymetrix protocols (17–20). These datasets were selected, in particular, for maximal coverage of iron-associated genes (some of which are found exclusively on the U133B or U133 PLUS 2.0 platforms, but not on the U133A GeneChip). Of 759 cases, 674 were annotated for distant metastasis-free survival (DMFS) with a minimum event-free follow up of 5-years. Relevant cohort details are summarized in Supplementary Table S1. For all cases, raw data (CEL files) were quantile-normalized within cohorts by computing Robust Multichip Average (RMA) expression values [Bioconductor “Affy Package (R)”] as previously described (12). Batch effects between cohorts were estimated by mixed model ANOVA and corrected using the Partek Genomics Suite Batch Remover program. The complete matrix of RMA-normalized, batch-corrected microarray expression data and corresponding clinicopathologic tumor annotations can be found in Supplementary Datafiles S1 and S2, respectively. Intrinsic molecular subtypes of breast cancer were assigned to samples using the Single Sample Predictor (SSP) algorithm as described (13) and utilized by Fan and colleagues (14). Affymetrix probe sets were matched to the genes comprising the SSP centroids. Prior to batch correction, the signal intensities for each gene were mean centered, and Spearman correlation was used to find the SSP centroid most closely associated with each tumor sample. In cases where a correlation greater than 0.1 was not achieved with at least one centroid, a subtype was not assigned to that sample.
Gene survival analysis and construction of prognostic classifiers
The 674 cases of the MMC were randomized into 2 groups: a “training” and a “testing” cohort (n = 337 in each), having similar Kaplan–Meier estimated survival curves (i.e., with nonsignificant differences) and near equal proportions of cases (<10% variation) with respect to cohort, ER status, nodal status, and treatment group. For the Fisher's exact and permutation testing analysis, Affymetrix probes sets were as given: (i) annotated according to the latest Affymetrix' probe set annotation (11-15-2009), (ii) filtered to exclude those lacking gene names/symbols (n = 35,515 remaining), (iii) analyzed for DMFS associations by Cox proportional hazards regression (Partek Genomics Suite), (iv) ranked by Cox likelihood ratio test P value, and (v) filtered to removed redundancy (i.e., probe sets targeting the same gene were filtered such that the one with the most significant P value remained; n = 18,428 remaining). P values associated with Fisher's exact test and permutation testing (100,000 iterations) were generated using custom JavaScript. To construct prognostic models, the supervised principal component method of Bair and Tibshirani (15) was used as implemented in the BRB-ArrayTools package developed by Dr. Richard Simon (National Cancer Institute) and the BRB-ArrayTools Development Team. Briefly, the supervised principal component method was employed for the following 2 reasons: (i) it is one of few well-validated methods that incorporates survival time information in the model training process, and (ii) the algorithm generates a scalable prognostic index for each test case. All models were generated using doubly nested 10-fold cross-validation. The prognostic index for the IRGS can be computed by the formula ∑iwi xi + 7.952256, where wi and xi are the weight and logged gene expression for the i-th gene. A new sample is then predicted as high risk if its prognostic index is larger than 0.265029, medium risk if the index is within −0.361502 and 0.265029, and low risk if the index is smaller than or equal to −0.361502. The IRGS gene weights can be found in Supplementary Datafile S3. All Kaplan–Meier and multivariable analyses were conducted in Sigma Plot 11.0.
Results
Genes that regulate iron biology exhibit prognostic associations in breast cancer
To investigate the prognostic impact of iron regulatory genes in breast cancer, we assembled a large retrospective cohort of breast cancer cases previously profiled on Affymetrix microarrays. The cohort consists of 759 breast tumors representing 6 independent, previously published expression microarray datasets (Supplementary Table S1). Of the 759 cases, 674 were annotated for DMFS with a minimum event-free follow up of 5 years. To construct and confirm prognostic models based on iron-associated genes, we randomized cases to 2 groups with similar population characteristics: a training and a test cohort, each comprising 337 cases (Supplementary Table S2). Next, we assembled a comprehensive list of genes with known functions in regulating iron biology. This list was derived from Gene Ontology (GO) categories related to iron metabolism (16, 17) and review of the literature; it comprised 63 genes and overlapped largely with the iron network devised by Hower and colleagues (18). Sixty-one of these genes could be mapped to one or more corresponding microarray probe sets found on the Affymetrix U133A or U133B GeneChips. Using the training cohort, we examined statistical associations between expression patterns of the iron regulatory genes and patient DMFS. Strikingly, we found that 49% of the genes were significantly associated with DMFS by Cox regression (P < 0.05; likelihood ratio test) as measured by one or more probe sets (Supplementary Table S3). To determine the likelihood that such an observation would occur by chance alone, we analyzed the training cohort as well as the entire combined cohort by Fisher's exact test and permutation testing (Table 1). Using 2 P value thresholds for assigning significance to DMFS-associated genes (P < 0.05 or P < 0.01), we observed by Fisher's exact test a statistically significant enrichment for DMFS-associated genes among the 61 iron regulatory genes (as compared with the remaining population of genes represented on the microarray). This statistical significance was observed at both P value thresholds and in both cohorts. Next, we carried out permutation testing whereby in each cohort, we randomly selected a 61-gene set 100,000 times, and at each iteration, counted the number of DMFS-associated genes at each P value threshold. In both cohorts, and at both thresholds, we again observed a statistically significant enrichment for DMFS-associated genes among the iron regulatory genes (P ≤ 0.02 in all cases). These data indicate that iron regulatory genes, as a group, are statistically unique in their tendency to be associated with DMFS, perhaps reflecting a pathologic role for iron regulation in the clinical behavior and progression of breast cancer.
. | Training cohort (n = 337) . | Combined cohort (n = 674) . | ||
---|---|---|---|---|
. | P < 0.05 . | P < 0.01 . | P < 0.05 . | P < 0.01 . |
No of iron genes (of n = 61) significant: | 30 | 16 | 34 | 26 |
No of all genes (of n = 18,428) significant: | 5,718 | 3,012 | 7,778 | 4,963 |
Fisher's exact Pa | 0.003 | 0.03 | 0.02 | 0.006 |
Permutation Pb | 0.0009 | 0.02 | 0.01 | 0.003 |
. | Training cohort (n = 337) . | Combined cohort (n = 674) . | ||
---|---|---|---|---|
. | P < 0.05 . | P < 0.01 . | P < 0.05 . | P < 0.01 . |
No of iron genes (of n = 61) significant: | 30 | 16 | 34 | 26 |
No of all genes (of n = 18,428) significant: | 5,718 | 3,012 | 7,778 | 4,963 |
Fisher's exact Pa | 0.003 | 0.03 | 0.02 | 0.006 |
Permutation Pb | 0.0009 | 0.02 | 0.01 | 0.003 |
aTwo-tailed.
b100,000 iterations.
Iron regulatory gene signatures are predictive of breast cancer recurrence
Using the supervised principal component method of Bair and Tibshirani (15), we investigated the prognostic potential of iron regulatory genes in statistical models predictive of breast cancer recurrence. The training cohort was used to generate, in total, 9 prognostic models.
To construct the models, we considered 3 threshold significance levels (α, alpha) for feature selection (α = 0.01, α = 0.001, or α = 0.0001). The threshold significance level refers to the statistical significance with which gene expression patterns associate with DMFS of breast cancer patients. At α = 0.01, α = 0.001, and α = 0.0001, the number of genes selected for model inclusion after cross-validation were 16 genes (19 probe sets), 7 genes (8 probe sets), and 4 genes (5 probe sets), respectively (Note that some genes are represented by more than one microarray probe set).
At each alpha threshold significance level, we also considered 1, 2, or 3 principal components for assigning gene weights within each prognostic model (see Materials and Methods). Principal components refer to linear combinations of iron regulatory genes that explain the variance–covariance structure underlying the gene expression patterns. The combination of 3 alpha threshold significance levels and 3 principal components yielded 9 prognostic models (Supplementary Table S4).
Each model generates a prognostic index for each tumor sample, and this index reflects the relative likelihood of metastasis-free survival. During model construction, we specified 3 risk groups: low risk, intermediate risk, and high risk on the basis of prognostic index cut points set at the 33rd and 66th percentiles of the training cohort. After each model was generated, we applied it to the test cohort for independent verification of model performance.
In all training and testing scenarios, the prognostic models stratified cases into predicted risk groups with significantly different survival rates (P < 0.005; log-rank test; Supplementary Table S4), indicating that the various different combinations of iron regulatory genes have robust prognostic potential regardless of model parameters. Maximal risk stratification in the training cohort, as determined by the HRs in Supplementary Table S4, was achieved with models comprising 16 genes (models #2 and #3; α = 0.01); between these, the model with 2 principal components was selected for further analysis (model #2; α = 0.01, 2 principal components), because this model exhibited the greatest risk group stratification at the clinically relevant endpoint of 10 years (data not shown). The genes comprising this model and their univariable associations with DMFS are shown in Table 2. As shown in Fig. 1, Kaplan–Meier survival plots show a robust association between this collection of 16 iron regulatory genes and distant metastasis-free survival.
Probe set IDs . | Symbol . | Name . | Pa . | HRb . | 95% CI . |
---|---|---|---|---|---|
222453_at; 217889_s_at | CYBRD1 | Cytochrome b reductase 1 | 1.83E-07 | 0.60 | 0.5–0.73 |
205542_at | STEAP1 | Six transmembrane epithelial antigen of the prostate 1 | 4.21E-06 | 0.59 | 0.47–0.74 |
225871_at | STEAP2 | Six transmembrane epithelial antigen of the prostate 2 | 2.02E-05 | 0.60 | 0.48–0.76 |
206087_x_at | HFE | Hemochromatosis | 3.05E-04 | 0.34 | 0.19–0.61 |
229839_at; 235849_at | SCARA5 | Scavenger receptor class A, member 5 (putative) | 4.02E-04 | 0.44 | 0.28–0.69 |
202018_s_at | LTF | Lactotransferrin | 4.16E-04 | 0.84 | 0.76–0.93 |
240686_x_at | TFRC | Transferrin receptor (p90, CD71) | 6.16E-04 | 3.54 | 1.72–7.3 |
223044_at; 233123_at | SLC40A1 | Solute carrier family 40 (iron-regulated transporter), member 1 (Ferroportin) | 7.00E-04 | 0.76 | 0.64–0.89 |
209075_s_at | ISCU | Iron-sulfur cluster scaffold homolog | 7.74E-04 | 0.41 | 0.24–0.69 |
218392_x_at | SFXN1 | Sideroflexin 1 | 8.33E-04 | 2.02 | 1.34–3.06 |
200878_at | EPAS1 | Endothelial PAS domain protein 1 | 1.07E-03 | 0.57 | 0.41–0.8 |
226179_at | SLC25A37 | Solute carrier family 25, member 37 | 2.13E-03 | 0.55 | 0.37–0.80 |
209735_at | ABCG2 | ATP-binding cassette, subfamily G (WHITE), member 2 | 3.82E-03 | 0.45 | 0.26–0.77 |
241999_at | SFXN5 | Sideroflexin 5 | 5.41E-03 | 4.00 | 1.51–10.6 |
226648_at | HIF1AN | Hypoxia inducible factor 1, alpha subunit inhibitor | 6.53E-03 | 0.51 | 0.32–0.83 |
218487_at | ALAD | Aminolevulinate dehydratase | 8.68E-03 | 0.49 | 0.29–0.84 |
Probe set IDs . | Symbol . | Name . | Pa . | HRb . | 95% CI . |
---|---|---|---|---|---|
222453_at; 217889_s_at | CYBRD1 | Cytochrome b reductase 1 | 1.83E-07 | 0.60 | 0.5–0.73 |
205542_at | STEAP1 | Six transmembrane epithelial antigen of the prostate 1 | 4.21E-06 | 0.59 | 0.47–0.74 |
225871_at | STEAP2 | Six transmembrane epithelial antigen of the prostate 2 | 2.02E-05 | 0.60 | 0.48–0.76 |
206087_x_at | HFE | Hemochromatosis | 3.05E-04 | 0.34 | 0.19–0.61 |
229839_at; 235849_at | SCARA5 | Scavenger receptor class A, member 5 (putative) | 4.02E-04 | 0.44 | 0.28–0.69 |
202018_s_at | LTF | Lactotransferrin | 4.16E-04 | 0.84 | 0.76–0.93 |
240686_x_at | TFRC | Transferrin receptor (p90, CD71) | 6.16E-04 | 3.54 | 1.72–7.3 |
223044_at; 233123_at | SLC40A1 | Solute carrier family 40 (iron-regulated transporter), member 1 (Ferroportin) | 7.00E-04 | 0.76 | 0.64–0.89 |
209075_s_at | ISCU | Iron-sulfur cluster scaffold homolog | 7.74E-04 | 0.41 | 0.24–0.69 |
218392_x_at | SFXN1 | Sideroflexin 1 | 8.33E-04 | 2.02 | 1.34–3.06 |
200878_at | EPAS1 | Endothelial PAS domain protein 1 | 1.07E-03 | 0.57 | 0.41–0.8 |
226179_at | SLC25A37 | Solute carrier family 25, member 37 | 2.13E-03 | 0.55 | 0.37–0.80 |
209735_at | ABCG2 | ATP-binding cassette, subfamily G (WHITE), member 2 | 3.82E-03 | 0.45 | 0.26–0.77 |
241999_at | SFXN5 | Sideroflexin 5 | 5.41E-03 | 4.00 | 1.51–10.6 |
226648_at | HIF1AN | Hypoxia inducible factor 1, alpha subunit inhibitor | 6.53E-03 | 0.51 | 0.32–0.83 |
218487_at | ALAD | Aminolevulinate dehydratase | 8.68E-03 | 0.49 | 0.29–0.84 |
Abbreviation: CI, confidence interval.
aLikelihood ratio test P.
bValues <1.0 indicate that expression is positively correlated with good survival.
To evaluate the potential clinical impact of the IRGS, we compared its prognostic attributes with those of conventional markers of breast cancer recurrence in a multivariable analysis involving the test cohort (Table 3). In the presence of nodal status, tumor size, patient age, histologic grade, and ER status, the IRGS remained statistically significant (P = 0.005) indicating that the IRGS contributes additive prognostic information not captured by these conventional markers.
. | Univariable analysis . | Multivariable analysis . | ||||
---|---|---|---|---|---|---|
Covariates . | Cox P valueb . | HR . | 95% CI . | Cox P valueb . | HR . | 95% CI . |
IRGS (low-, mid-, high-risk) | <0.001 | 1.85 | 1.37–2.51 | 0.005 | 1.59 | 1.15–2.21 |
ER status (0, 1) | 0.944 | 0.98 | 0.49–1.96 | 0.68 | 1.16 | 0.56–2.45 |
LN status (0, 1) | 0.005 | 1.93 | 1.22–3.05 | 0.121 | 1.48 | 0.90–2.42 |
Tumor size (cm; <2, 2-3, >3) | <0.001 | 1.91 | 1.37–2.66 | 0.019 | 1.56 | 1.08–2.27 |
Patient age (<40, 40–65, >65) | 0.647 | 1.10 | 0.72–1.69 | 0.373 | 1.22 | 0.79–1.89 |
Histologic grade (1, 2, 3) | <0.001 | 1.75 | 1.26–2.43 | 0.207 | 1.29 | 0.87–1.91 |
. | Univariable analysis . | Multivariable analysis . | ||||
---|---|---|---|---|---|---|
Covariates . | Cox P valueb . | HR . | 95% CI . | Cox P valueb . | HR . | 95% CI . |
IRGS (low-, mid-, high-risk) | <0.001 | 1.85 | 1.37–2.51 | 0.005 | 1.59 | 1.15–2.21 |
ER status (0, 1) | 0.944 | 0.98 | 0.49–1.96 | 0.68 | 1.16 | 0.56–2.45 |
LN status (0, 1) | 0.005 | 1.93 | 1.22–3.05 | 0.121 | 1.48 | 0.90–2.42 |
Tumor size (cm; <2, 2-3, >3) | <0.001 | 1.91 | 1.37–2.66 | 0.019 | 1.56 | 1.08–2.27 |
Patient age (<40, 40–65, >65) | 0.647 | 1.10 | 0.72–1.69 | 0.373 | 1.22 | 0.79–1.89 |
Histologic grade (1, 2, 3) | <0.001 | 1.75 | 1.26–2.43 | 0.207 | 1.29 | 0.87–1.91 |
a292 cases had complete clinical annotation for LN status, tumor size, patient age, and histologic grade.
bLikelihood ratio test P value.
The iron regulatory gene signature is prognostic of DMFS in breast tumor subtypes and treatment groups
To further understand the role of the IRGS in breast cancer, we evaluated its prognostic significance in specific molecular subtypes of breast cancer and patient treatment groups. Of note, we limited this analysis to subsets of the independent test cohort, to avoid overly optimistic interpretations that may result from inclusion of the training data.
First, we separately examined IRGS performance in the 295 cases positive for estrogen receptor (ER+) and the 40 cases lacking estrogen receptor (ER−). In the ER+ cases (Fig. 2A), the 3 risk groups significantly stratified patients by DMFS (P < 0.001; log-rank test). The significance of this stratification was largely driven by the particularly good outcome of the predicted low-risk cases as compared with the intermediate- and high-risk groups, which showed similar, and comparably poorer, survival rates. In the smaller ER− population (Fig. 2B), the IRGS showed a consistent trend toward the correct classification of low-, intermediate-, and high-risk groups, but did not achieve statistical significance. Notably, however, we observed that the prognostic index assigned by the IRGS to the ER−cases, as a continuous variable, was in fact statistically significantly associated with DMFS by Cox regression analysis (P = 0.035; likelihood ratio test), indicating that alternate IRGS cut points may be necessary for assigning risk to ER− cases.
We next examined the prognostic relevance of the IRGS in the intrinsic breast cancer subtypes (19, 20). First, we examined the basal subtype, which tends to comprise ER− tumors with poor outcomes (13, 19, 20). Consistent with this poor outcome association, the majority of cases (67%) were assigned by the IRGS to the predicted high-risk group (Fig. 2C). Similar to that observed of the ER− population, the IRGS showed a substantial but nonsignificant risk stratification of the basal cases. The continuous prognostic index, however, did not reach significance by Cox regression (P = 0.1).
Next, we considered the subtypes that largely comprise ER+ breast cancer, namely, the luminal A (LumA), luminal B (LumB), and Normal-like (NL) subtypes. Consistent with previous observations that LumA and NL subtypes exhibit more favorable survival outcomes (13), the IRGS classified the majority of LumA (54%; Fig. 2D) and NL subtypes (55%; Fig. 2E) into the low-risk group. However, within the LumA subtype, the IRGS predicted intermediate- and high-risk cases that showed significantly poorer survival (P = 0.03; log-rank test), indicating that the IRGS can further risk stratify LumA disease. In a similar fashion, the IRGS further risk stratified the NL subtype (P = 0.01; log-rank test) predicting a small fraction of high-risk cases (4%) that showed an increased rate of distant metastasis. The LumB subtype has been historically associated with poor survival outcomes, and concordantly in our patient population, the majority of LumB cases were classified by the IRGS as high risk (69%; Fig. 2F). However, both the predicted low- and intermediate-risk groups showed similarly poor outcomes, indicating that the IRGS likely does not have prognostic utility within this tumor subtype. An additional subtype, the HER2+-like subtype, has also been described; however, too few cases were present (n = 24) to assess performance of the IRGS. However, we did note that all but one of the HER2+-like cases were classified by the IRGS as high or intermediate risk, consistent with previous observations that this subtype is associated with poor outcome (Supplementary Datafile S2). Taken together, these observations suggest that the IRGS, while recapitulating some of the prognostic features of the molecular subtypes, may provide valuable additive prognostic information to the LumA and NL subtypes, and potentially ER- negative breast cancer, although further validation in larger breast cancer populations is necessary.
As the potential effects of different treatments were not accounted for in the previous analyses, we examined the prognostic attributes of the IRGS in uniformly treated patients. First, we considered the subset of 104 patients (in the test cohort) who, following surgery, received no adjuvant systemic therapy (Fig. 3A). The IRGS risk groups stratified these patients with statistical significance (P = 0.05; log-rank test), indicating that the IRGS has, to some extent, a purely prognostic component uncoupled from adjuvant therapy prediction. The largest, most uniformly treated patient subgroup in the test cohort comprised ER+ patients [negative or positive for lymph node (LN) involvement] treated in the adjuvant setting with tamoxifen monotherapy (n = 185). Determining treatment for these patients is a particular clinical challenge, as the desire to treat these patients aggressively with combination tamoxifen and adjuvant chemotherapy is counterbalanced by the small gain in therapeutic benefit imparted by added chemotherapy and the severity of the adverse side effects caused by chemotherapy. We evaluated the IRGS in 2 subsets of this treatment group: the ER+, LN− subset (n = 99) and the ER+, LN+ subset (n = 86). As shown in Fig. 3B and C, the IRGS predicted for a low-risk group that exhibits significantly better DMFS than the predicted intermediate- and high-risk groups, and this observation is consistent in both the LN− and LN+ populations. As the intermediate- and high-risk groups showed no differences in actual DMFS rates, we considered the relevance of the IRGS as a binary classifier (i.e., 0 = low risk and 1 = high risk/intermediate risk) in a multivariable analysis of the tamoxifen-treated population. Although 4 covariates (IRGS, LN status, tumor size, and histologic grade) were significantly associated with DMFS by univariate analysis, only the IRGS and tumor size remained significant in the multivariable model (P = 0.007 and P = 0.035, respectively) (Supplementary Table S5). Together, these observations show that the IRGS contributes important clinical value in predicting subsets of ER+ patients (even those with LN+ disease) that will show excellent long-term distant metastasis-free survival if they are treated with tamoxifen and spared adjuvant chemotherapy.
Iron export and iron import gene dyads are complementing prognostic factors in breast cancer
To better understand the transcriptional dynamics of the prognostic IRGS genes, we investigated their correlation structure in the combined breast cancer cohort by hierarchical clustering (Supplementary Fig. S1). Surprisingly, we found that the transcriptional patterns of the IRGS genes are largely diverse, displaying an average Pearson correlation of −0.1. This suggests that the IRGS genes may represent multiple regulatory pathways, each of which affects iron homeostasis in an independent way.
To test this hypothesis, we examined the effect of transferrin receptor 1 (TFRC) and hereditary hemochromatosis (HFE) on prognosis. The rationale for this approach was derived from our previous observation that expression of ferroportin and hepcidin, 2 genes whose products work together to regulate iron export, affect prognosis in breast cancer (10). Specifically, we had previously found that high levels of the iron efflux pump ferroportin (SLC40A1), which leads to low levels of intracellular iron, were associated with favorable prognosis. In patients who expressed high levels of SLC40A1, concomitant expression of low levels of hepcidin (HAMP), a protein that degrades ferroportin, further improved prognosis. Because these results suggest that decreased intracellular iron is associated with a favorable prognosis, we reasoned that other genes that decrease intracellular iron might similarly affect prognosis. We therefore tested whether TFRC and HAMP, 2 genes whose products work together to regulate iron import, might represent a complementary regulatory pathway embedded in the IRGS gene set. Cellular uptake of iron is predominantly driven by endocytosis of iron-loaded transferrin bound to TFRC. The HFE protein negatively regulates TFRC-mediated iron uptake (21–23). Thus, whereas TFRC acts to promote iron import, HFE acts to block it, through mechanisms that are still under investigation (24).
To study the prognostic interplay between TFRC and HFE, as well as their associations with SLC40A1 and HAMP, we studied the expression–survival associations of these genes in the full combined cohort with a focus on the ER+ patients uniformly treated with tamoxifen monotherapy (n = 371). All tumors were assigned the binary annotation of “low” or “high” expression for a given gene on the basis of whether the signal intensity fell below or above the population mean. A pro-export phenotype was assigned to tumors having high SLC40A1 and low HAMP concomitantly, whereas an anti-export phenotype was assigned to those having both low SLC40A1 and high HAMP. In a similar vein, a pro-import phenotype was assigned to tumors showing both high TFRC levels and low HFE levels, whereas an anti-import phenotype was assigned to those having concomitant low TFRC and high HFE. Shown in Fig. 4 are distant metastasis-free survival estimates of breast cancer cases categorized according to iron export (Fig. 4A) and iron import (Fig. 4B) phenotypes. Consistent with the hypothesis that cellular iron content is a determinant of breast cancer behavior, both pro-export and anti-import phenotypes predicted for reduced metastasis rates, whereas anti-export and pro-import phenotypes predicted for significantly increased metastasis rates. Next, to investigate the prognostic relationship between the iron export and import phenotypes, each dyad was compared in a multivariable model for its prognostic contributions (Supplementary Table S6). Both dyads remained highly significant in the multivariable model (P < 0.005) indicating that they each contribute additively (in a noncolinear fashion) to the prediction of metastatic recurrence. Indeed, although both pro-export and anti-import phenotypes predicted for low metastatic risk, they largely comprised different patients (Fig. 4C), suggesting these phenotypes represent independent pathways to a final common endpoint. Together, these data present the possibility that different regulatory modulators of cellular iron content may directly, and through distinct mechanisms, impact the clinical progression of breast cancer.
Discussion
Datasets of tumor expression profiles represent a rich resource for hypothesis testing and may provide insights into novel biological properties of human tumors. Previous work has successfully used such datasets to identify groups of functionally unrelated genes that collectively associate with breast cancer risk (19, 20). In this work, we used a different strategy, capitalizing on publicly available microarray datasets derived from breast cancer patient cohorts to test the hypothesis that perturbations of iron metabolism affect breast cancer risk. Specifically, we tested (i) whether expression of genes related to iron metabolism are collectively linked to breast cancer prognosis, (ii) whether an optimal IRGS could be identified, (iii) whether specific pathways involving different aspects of iron management could be identified within this signature, and (iv) whether this signature exhibited potential clinical utility.
We observed a statistically significant association between almost 50% of 61 genes involved in iron metabolism and breast cancer prognosis, indicating a remarkably robust association between expression of genes related to iron metabolism and breast cancer prognosis (Table 1). Maximal stratification into low-, intermediate-, and high-risk groups was achieved with an IRGS comprising 16 of these genes, which therefore represents our optimal signature. The IRGS does not seem to simply recapitulate information provided in other molecular classifiers of breast cancer and provides additional useful information to discriminate among patients. For example, we found that the IRGS could further stratify LumA and normal-like tumors (Fig. 2D and E) into high-, intermediate-, and low-risk groups. Similarly, multivariable analysis revealed that the IRGS contributed information in addition to that provided by the conventional markers of nodal status, tumor size, patient age, histologic grade, or ER status.
The IRGS also contains embedded information about molecular pathways of iron metabolism that can be used to probe pathways that are perturbed in cancer. Our previous results identified ferroportin and hepcidin as 2 genes important in one such pathway (10). The products of these genes can be considered an iron efflux dyad: ferroportin is an iron efflux pump whose stability is controlled by hepcidin. We observed that a “low intracellular iron” phenotype conferred by high ferroportin and low hepcidin was associated with good prognosis, whereas a “high intracellular iron” phenotype conferred by low ferroportin and high hepcidin was associated with poorer prognosis (10).
In this study, we interrogated the IRGS to test whether pathways that mediate iron import might represent a second pathway that similarly modifies prognosis. We selected TFRC, the major iron importer in most mammalian cells, and HFE, a protein that negatively regulates TFRC-mediated iron uptake, to represent an iron import dyad. In agreement with results obtained with the iron export dyad of ferroportin and hepcidin, a “low intracellular iron” phenotype was associated with good prognosis, and a “high intracellular iron” phenotype was associated with poorer prognosis (Fig. 4). Thus, high levels of TFRC in conjunction with low levels of its negative modulator HFE were associated with poorer DMFS than low levels of TFRC and high levels of HFE. These results are consistent with the literature: an increase in levels of the TFRC protein in breast and other cancers has been known since the 1980s (25), as well as the dependency of cancer cell proliferation on TFRC-mediated iron uptake (26–28). In fact, the transferrin receptor is frequently used as a targeting ligand in the design of anticancer drugs (29). Overall, our results are concordant with and begin to provide molecular specificity to the historical view of iron as an element that favors both malignant transformation and tumor growth (30).
Our data may be useful for future hypothesis generation and testing. Not all gene associations that we observed conform to a simple picture in which increased iron content is associated with poorer prognosis: we found that expression of some genes whose products have been ascribed functions related to iron import (e.g., CYBRD1) were associated with improved rather than decreased survival (Table 2). Interestingly, the association of CYBRD1 with good prognosis is supported by a recent study that sought to identify genes that distinguish primary breast tumors and their matched nodal metastases. This study identified 4 of the 16 IRGS genes including CYBRD1 (the others are LTF, STEAP1, and STEAP2) in a shortlist of 431 named genes as having statistically significantly decreased transcript levels in metastases compared with the primary tumors (31). Such findings may suggest that products of these “iron” genes exhibit multiple functions, or that they play alternative roles in breast tissue when compared with tissues involved in systemic iron management (e.g., the duodenum, liver, etc), where the roles of many of these genes have been elucidated.
Some gene associations we observed confirm expectations based on the literature. For example, we found that HIF1A, which is frequently a negative prognostic indicator (32), was associated with reduced DMFS (Supplementary Table 3). Conversely, increased expression of the HIF1alpha inhibitor (HIF1AN) was associated with improved DMFS. It has been reported that breast cancer in younger women represents a unique biological entity, with increased expression of genes related to hypoxia and other signaling pathways (33). The relatively small number of young patients in our cohort did not allow us to discriminate potential special features of IRGS gene expression as a function of age and, in particular, whether expression of HIF1A and other genes related to hypoxia were linked to patient age. We also did not have access to information on anemia or dietary iron in these subjects. The contribution of factors such as age, diet, and anemia to the IRGS will require further study.
Assessing the expression of the 16 genes comprising the IRGS in breast cancer patients may have clinical utility. We observed that the IRGS could be successfully used to predict the outcome of ER+ patients (Fig. 2A) and patients exhibiting favorable molecular subtypes [LumA and normal-like (Fig. 2D and E)]. Whereas the IRGS showed only a near-significant trend toward risk stratification in ER− patients by Kaplan–Meier analysis (Fig. 2B), Cox regression showed a significant association with DMFS when the IRGS prognostic index was used as a continuous variable. This suggests that future study of the IRGS in risk stratification of ER− patients will require a rescaling of the survival group assignments as they pertain to the IRGS prognostic index. When we divided patients into groups that had received homogeneous treatment, the IRGS was able to stratify all patients, including those who had received no treatment and those treated with tamoxifen monotherapy, independent of lymph node status (Fig. 3). These results suggest that the IRGS may be useful in 2 different clinical settings. First, it may allow identification of potentially high-risk patients within the traditional low-risk LumA and NL subtypes, which comprise a sizable fraction of luminal breast cancer cases. Second, it may spare patients unneeded chemotherapy by allowing the identification of subsets of ER+ patients (even those with LN+ disease) who will show excellent DMFS when treated with tamoxifen monotherapy. That the IRGS significantly risk-stratified LumA but not LumB cases, suggests that an integrated prognostic model combining aspects of the IRGS and luminal subtypes [or the subtype-derived Risk of Relapse (ROR) score; refs. 34, 35] could provide additive prognostic power to outcome prediction.
Overall, this study shows a strong link between genes that govern iron metabolism and breast cancer, and suggests new tools to guide breast cancer prognosis. Over the longer term, it may also help uncover metabolic differences that distinguish normal and malignant breast cancer cells that can be used to therapeutic advantage.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Grant Support
The work was supported in part by grants R37 DK42412 (F.M. Torti) and R01DK071892 (S.V. Torti).
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.