Downregulation of CYB5D2 is associated with breast cancer progression

We report here that CYB5D2 is associated with tumor suppression function in breast cancer (BC). CYB5D2 expression was significantly reduced in tamoxifen resistant MCF7 cells and in MCF7 cell-derived xenografts treated with TAM. CYB5D2 overexpression induced apoptosis in MCF7 cells; CYB5D2 knockdown enhanced MCF7 cell proliferation. Using the TCGA and Curtis datasets within the Oncomine database, CYB5D2 mRNA expression was downregulated in primary BCs vs breast tissues and HER2-positive or triple negative BCs vs estrogen receptor (ER)-positive BCs. Using the TCGA and Metabric datasets (n = 817 and n = 2509) within cBioPortal, 660 and 4891 differentially expressed genes (DEGs) in relation to CYB5D2 were identified. These DEGs were enriched in pathways governing cell cycle progression, progesterone-derived oocyte maturation, oocyte-meiosis, estrogen-mediated S-phase entry, and DNA metabolism. CYB5D2 downregulation decreased overall survival (OS, p = 0.0408). A CYB5D2-derived 21-gene signature was constructed and robustly correlated with OS shortening (p = 5.72e-12), and independently predicted BC deaths (HR = 1.28; 95% CI 1.08–1.52; p = 0.004) once adjusting for known clinical factors. CYB5D2 reductions displayed relationship with mutations in PIK3CA, GATA3, MAP3K1, CDH1, TP53 and RB1. Impressively, 85% (560/659) of TP53 mutations occurred in the 21-gene signature-positive BC. Collectively, we provide the first evidence that CYB5D2 is a candidate tumor suppressor of BC.

cancer; its enforced expression reduced the in vitro invasion and in vivo lung metastasis of HeLa cells 11 . The CYB5D2 gene resides at 17p13.2; 17p13. 2-13.3 is lost in 50% of breast cancer 29 , indicating that CYB5D2 may be a novel tumor suppressor of BC.
In support of this possibility, we report here a significant reduction of CYB5D2 expression following the progression of tamoxifen resistance both in vitro and in vivo as well as in more than 3000 primary BCs. CYB5D2 downregulation is correlated with mutations in PIK3CA, GATA3, MAP3K1 and TP53 as well as reductions in overall survival (OS) of breast cancer.

Methods
Tissue culture and the development of tamoxifen-resistant cells. MCF7 cells were purchased from the American Type Culture Collection (ATCC) 30 . Cell were cultured in DMEM accompanied with 10% fetal bovine serum and 1% Penicillin-Streptomycin (Life Technologies, Burlington, ON). Tamoxifen resistant MCF7 (TAM-R) cells were developed by continuous culture of MCF7 cells in phenol-red-free DMEM media in the presence of 1 μM of 4-hydroxyl-tamoxifen (Sigma Aldrich, Oakville, ON) for 12 months 30 . The TAM resistance status was confirmed.
TUNEL apoptotic detection. MCF7 cells were seeded in chamber slide and transfected transiently with either GFP or GFP-CYB5D2 for 48 hours. TUNEL procedures were then carried out with a TUNEL kit (Abcam) following the manufacturer's instructions.
Knockdown of CYB5D2 and proliferation assay. MCF7 cells were transfected using lentivirus-based (control: CTRL) shCTRL or shCYB5D2 (a pool of three individual knockdown constructs; Santa Cruz); the knockdown was confirmed. MCF7 shCTRL and MCF7 shCYB5D2 cells (5 × 10 4 cells/well) were seeded in a 6-well tissue culture plate with cell numbers counted every three days for nine days.
Determination of TAM-derived cytotoxicity. Cells (10 5 cells/well) were first seeded into a 6 well plate with phenol-red free DMEM media, and cultured for 2 days prior to treatment with either 3 μM TAM or DMSO (1:1000) in serum-free media for 48 hours. Cells were then cultured in compete medium without TAM for 96 hours 31 , followed by fixation in a solution containing 2% formaldehyde and 0.2% glutaraldehyde for 20 minutes prior to addition of a crystal violet solution (0.5% Crystal violet, 20% methanol, 150 mM NaCl) for 30 minutes. The plates were then washed in water and allowed to dry, after which images were obtained. The staining was then released using 2 mL of 33% acetic acid for quantification by measuring absorbance at 550 nm 30 . Treatment of xenograft tumor with TAM. Four to five week old ovariectomized nude mice had an 0.72 mg estrogen pelle inserted. 3 × 10 6 MCF7 cells were implanted into the flank of each mouse after which animals were randomly divided into two groups. Half received a 5 mg tamoxifen pellet and the other half served as controls. Animals were maintained for 28 weeks or until endpoint was reached (1000 mm3). Tumor volumes were determined using a caliper according to the standard formula: L × W 2 × 0.52, where L and W are the longest and shortest diameters, respectively. All animal work was performed according to protocols approved by the McMaster University Animal Research Ethics Board 30 . β-actin was used as an internal control. CYB5D2 mRNA levels were normalized to those of β-actin. Experiments were repeated three times; means ± SD (standard deviation) were graphed. Statistical analysis was performed using Student's t-test (2-tails). (B) Western blot analysis of CYB5D2 protein expression in MCF7 and TAM-R cells. Experiments were repeated three times; typical images are provided (please see Fig S14 for the non-cropped Western blot). CYB5D2 protein expression was normalized to Actin; means ± SD are graphed; *p < 0.05 determined by 2-tailed Student's t-est. (C,D) MCF7 cell-derived xenograft tumors were generated in NOD/SCID mice, followed by treatment with and without TAM. CYB5D2 expression in treated (+TAM) and untreated (−TAM) xenograft tumors (n = 5 for each group) was determined by real time PCR (C) and IHC (D). Typical IHC images are presented. **p < 0.01 in comparison to untreated tumors (2-tailed Student's t-test). (E,F) CYB5D2 mRNA expression data were extracted from the datasets of TCGA (E) 34 and Curtis (F) 6 . Mean ± SD for the indicated BC subtypes are www.nature.com/scientificreports www.nature.com/scientificreports/ 5′-ACCGAGCGCGGCTACAG-3′; R: 5′-CTTAATGTCACGCACGATTTCC-3′). Reactions were performed using the ABI 7500 Fast Real-Time PCR System (Applied Biosystems, Burlington, ON) in the presence of SYBR-green. All samples were run in triplicate 30 . Measurement of CYB5D2 mRNA levels. The mRNA expression data of CYB5D2 were extracted from TCGA 34 , Curtis 6 , Finak 35 , and Karnoub 36 within Oncomine (Compendia Bioscience, Ann Arbor, MI). The genomic mutational data of TP53, PIK3CA, GATA3, MAP3K1, and other genes were retrieved from the Curtis and TCGA-Cell datasets within cBioPortal (http://www.cbioportal.org/) 37,38 . A variety of statistical methods were used to examine CYB5D2 expression and its association with OS (see the section of Statistical analysis). Data related to CYB5D2-associated co-alteration of mutations and gene expression were extracted from the Metabric (n = 2509) and TCGA-Cell (n = 817) datasets within the cBioPortal database 39 . A signature consisting of CYB5D2 reduction and a set of genomic mutations was derived using the Cox regression model. Furthermore, a signature consisting of 21 genes was established using differentially expressed genes (DEGs) related to CYB5D2 downregulation; this was achieved by inputting these DEGs into the Cox model to select for their contributions to hazard ratio (HR) by either forward addition or backward elimination using SPSS Statistics version 23 40 .
Pathway enrichment analysis. The R packages of Reactome and GAGE as well as Ingenuity Pathway Analysis (IPA) were selected to analyze pathways that were enriched in the differentially expressed genes (DEGs) with the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.
Statistical analysis. Statistical analysis was performed using Student's t-test. Kaplan-Meier surviving curves, log-rank test, receiver-operating characteristic (ROC) curve, univariate and multivariate Cox proportional hazards regression analyses 40 were performed using the R survival package and SPSS Statistics version 23. p < 0.05 is considered statistically significant 39 .

Downregulation of CYB5D2 during BC tumorigenesis. Our recent identification of CYB5D2
as a tumor suppressor in cervical cancer 11 as well as its genomic location at 17p13.2; this region is frequently lost in breast cancer 29 ; led us to study the possible involvement of CYB5D2 in BC. We have recently derived a Tamoxifen-resistant (TAM-R) line from MCF7 cells (Fig. S1), and detected reductions of CYB5D2 mRNA and protein expression in TAM-R cells (Fig. 1A,B). The downregulation was also demonstrated in MCF7 cell-derived xenograft tumors treated with TAM compared to xenografts produced in untreated mice by either real-time PCR (Fig. 1C) or IHC (Fig. 1D). Addition of TAM decreased MCF7 empty vector (EV) cell-derived xenograft tumor growth 31 . By taking advantage of large datasets of BC gene expression available from the Oncomine database, we showed a significant decrease in CYB5D2 mRNA in BC compared to normal breast tissues in two large patient cohorts (Fig. 1E,F) and two small BC populations (Fig. S2). A further decrease in CYB5D2 expression was demonstrated in aggressive versus less aggressive BC sub-types: ER-versus ER+, PR-versus PR+, and HER2+or TN versus ER+in 2 large BC cohorts, the TCGA 34 and Curtis datasets 6 ( Fig. 1E,F). In both cohorts, CYB5D2 downregulations separate BC from breast tissues with area under curve (AUC) 0.712 and 0.696, respectively (Fig. 1G,H). Furthermore, CYB5D2 downregulation is associated with a rapid course of overall survival (OS) reduction in BC, including either ER+ or PR+BC ( Fig. S3A-C).
We subsequently characterized CYB5D2 downregulation-associated decreases in OS with the Curtis data (n = 1980) (cBioPortal) 6,41 . The dataset was retrieved from cBioPortal; tumors/patients were subsequently separated into a group without CYB5D2 reduction and a group with CYB5D2 downregulation at the levels below 1 standard deviation (SD, −1SD), −1.5 SD, or −2SD from the reference population mean which was derived from tumors that are diploid for the intended gene or the tumor population (http://www.cbioportal.org/faq.jsp). In comparison to BCs without the CYB5D2 downregulations, those with the downregulations are associated with significant reductions in OS within the initial follow-up period of 120 months (Fig. S4).
To investigate the impact of CYB5D2 in BC tumorigenesis, we have made a major effort to stably express CYB5D2 in MCF7 cells. This approach seemed feasible, as CYB5D2 was stably expressed in HeLa cells 11 . However, despite multiple tries by two individuals, a MCF7 cell line stably expressing CYB5D2 could not be established, suggesting that CYB5D2 potently inhibits MCF7 cell viability or proliferation. To study this possibility, we transiently expressed a CYB5D2-GFP (green fluorescence protein) fusion protein or GFP in MCF7 cells, and noticed a large number of rounded-up CYB5D2-GFP cells compared to GFP cells ( Fig. 2A,B), indicative of possible apoptosis in MCF7 cells expressing CYB5D2-GFP. Indeed, we detected 29.2% and 53.4% of the CYB5D2-GFP cells as TUNEL-positive 24 hours and 48 hours following transient transfection, respectively (Figs 2C, S5). To examine the effects of CYB5D2 in MCF7 cell growth, CYB5D2 was knocked-down in MCF7 cells (Fig. 2D, inset); knockdown of CYB5D2 significantly enhanced MCF7 cell proliferation (Fig. 2D). The ER + MCF7 cells require ER signaling to survival. Of note, CYB5D2 significantly reduced ER-derived transcription activity based on the luciferase activity driven by an ER enhancer reconstituted SV40 promoter (Fig. 2E). Taken together, we provide the first in vitro, in vivo, as well as clinical evidence suggesting functional downregulations of CYB5D2 during BC tumorigenesis.
graphed. *p < 0.05 by an unpaired, two-tailed, welch-corrected t-test; *p < 0.05 in comparison to normal breast tissues (Breast); $ p < 0.05 in comparison to the respective ER+ and PR+ tumors; and # p < 0.05 in comparison to ER+ breast tumors. (G,H) Receiver-operating characteristic (ROC) curves of normal breast tissues versus primary breast tumors were derived from the indicated datasets. AUC: area under the curve. (2019) 9:6624 | https://doi.org/10.1038/s41598-019-43006-y www.nature.com/scientificreports www.nature.com/scientificreports/ CYB5D2 downregulation associates with mutations in the major BC contributing genes. It becomes clear that combination of alterations in gene expression with those in genome more precisely revealed critical events of tumor evolution 42 . For incidence, gene expression plus copy number changes resulted in detailed clarification of BC into 10 sub-groups 6 . We thus analyzed CYB5D2 downregulation-associated genomic alterations. With CYB5D2 reduction at the above three levels (Fig. S4), there were no significant co-alterations in copy number variations defined at q-value (false discovery rate) < 0.05 in the Metabric dataset (n = 2509, cBioPortal). with either CYB5D2-GFP or GFP. Cells were imaged daily at 5 randomly selected fields. The transfection efficiency of CYB5D2-GFP was comparable to that of GFP. Experiments were repeated three times; typical images for one repeat are shown. (B) All GFP-positive cells and possible apoptotic cells (see white arrows for typical apoptotic cells) in 5 randomly selected fields were counted; no less than 300 GFP-positive cells for each transfection per time point were counted. Rounded cells (potential apoptotic cells) were then calculated and graphed; *p < 0.05 (2-tailed Student's t-test) in comparison to the respective GFP transfection. (C) MCF cells were transiently transfected with GFP or CYB5D2-GFP for 24 and 48 hours. TUNEL staining was then performed. Cells positive for both TUNEL (red) and GFP as well as positive for GFP only were counted from 5 randomly selected fields; a total of 300-400 cells for each cell type were counted. TUNEL-positive cells are expressed as % of GFP-positive cells. Experiments were repeated three times. Means ± SD are graphed. *p < 0.05 (2-tailed Student's t-test) in comparison to the respective GFP transfection. (D) Stable knockdown of CYB5D2 in MCF7 cells using either shCTRL (control) or shCYB5D2 lentivirus (inset). MCF7 shCTRL and MCF7 shCYB5D2 cells were seeded in 6-well tissue culture plates (5 × 10 4 cells/well); cell numbers were counted every three days. Experiments were repeated three times. Means ± SD are graphed. *p < 0.05 (2-tailed Student's t-test) in comparison to MCF7 shCTRL cells. (E) 293T cells in 24-well tissue culture plates were transiently transfected in triplicates with a cocktail containing an ER enhancer reconstituted promoter luciferase plasmid plus a β-Gal (galactosidase) vector together with either GFP or CYB5D2 expression plasmid. Luciferase activity was normalized to that of β-Gal. Experiments were repeated three times; means ± SD are graphed; **p < 0.01 in comparison to GFP by 2-tailed Student's t-test (2-tails).
We then examined whether these genomic alterations will strengthen the effects of CYB5D2 downregulation on OS shortening. Based on the association of CYB5D2 downregulations with an OS decrease (Fig. S4), CYB5D2 reduction at the −1.5 SD level was chosen for further analyses. By selection for contributions to CYB5D2 downregulation-associated decreases in OS, we established a signature (Signature #1) consisting of CYB5D2 reduction and the mutations in TP53, CDH1, BRCA1, THSD7A, BIR6, and RB1 (Table S2; Fig. 3A). Signature #1 significantly correlates with a reduction of OS in BC and ER-positive breast cancer (the Curtis dataset, n = 1980) (Fig. 3B,C). Mutations in TP53 occur most frequently in Signature #1 (Fig. 3A) and contribute to Signature #1's correlation with OS reductions. Removal of TP53 decreased the signature's potency, nonetheless, the signature retains association with OS shortening (control cases n = 1145, deaths n = 640, median months survival 169, and 95% CI: 159-181; risk individuals n = 361, deaths n = 241, median months survival 124, and 95% CI: 114-149; p = 4.33e-5). Removal of other individual components also decreased the association (data not shown), supporting their unique contributions to Signature #1. Additionally, Signature #1 independently predicts the risk of BC fatality (HR 1.328, 95% CI 1.131-1.560, p = 5.3e-4) following adjusting for cellularity, age at diagnosis, Neoplasm Histologic Grade, Integrative Cluster, tumor size, Nottingham prognostic index, and tumor stage. The signature remains an independent risk factor for BC deaths after removal of TP53 (HR 1.217, 95% CI 1.041-1.422, p = 0.01379). Additionally, Signature #1 associates with decreases in OS and DFS (disease-free survival) during 80-months follow-up in an independent TCGA-Cell cohort (n = 817) 34 (Fig. S6). Collectively these findings suggest an important connection between CYB5D2 and other oncogenic factors involved in BC pathogenesis.

Identification of CYB5D2-related genes.
To further characterize CYB5D2's involvement in BC pathogenesis, we have determined the DEGs (differentially expressed genes) relevant to CYB5D2 downregulation (−1.5 SD) in both the TCGA-Cell (n = 817) 34 and Metabric (n = 2509) datasets within cBioPortal. DEGs have been selected at q-value < 0.001. In the Metabric and TCGA-Cell datasets, 4981 and 660 DEGs were respectively  (Table S3C). Furthermore, these 471 genes display the same directionality of alterations (downregulation and upregulation) in both datasets (data not shown).
Building a 21-gene signature. To further analyze those DEGs, we selected DEGs that displayed alterations no less than that of CYB5D2. From the Metabric-DEGs, 98 genes were selected with 64 downregulation and 34 upregulation (Table S3D). We then imputed these DEGs in the Cox model using either addition (forward) or elimination (backward) of covariates to select DEGs based on their contributions to hazard ratio (HR). A 21 gene signature (Signature #2) was the result (Table 1; Fig. S7). Signature #2 robustly correlates with decreases of OS in BCs in the Curtis cohort (n = 1980, p = 5.72e-12) and the ER-positive subpopulation (n = 1560, p = 9.32e-12) (Fig. 4A,B). In the ER-negative subpopulation, Signature #2 is on a border line (p = 0.077) for a significant association with OS reductions (Fig. S8A). We noticed a large proportion of ER-negative tumors (432 of 474) being positive for Signature #2 (Fig. S8A), suggesting that this imbalance in distribution was a cause for the observed non-significant association (Fig. S8A).
We subsequently examined the impact of Signature #2 on OS in patients with PAM50-classified intrinsic subtypes using the Curtis dataset (cBioPortal). As the luminal subtypes are essentially ER-positive breast tumors 2-5 , we have focused on other intrinsic subtypes: claudin-low, normal-like, basal-like, and HER2-enriched BCs. Among 209 basal-like BCs, 198 were Signature #2-positive; of that, 110 patients died. Signature #2-positive basal-like BCs exhibited comparable OS compared to those of signature-negative (p = 0.659). In HER2-enriched and particularly claudin-low BCs within the follow up period of 160 months, Signature #2-positive BCs display a significantly shorter OS (Fig. S8b,C, right panel). Impressively, Signature #2 significantly correlates with OS decreases in normal-like BCs (Fig. 4C). Similarly, in the largest TCGA provisional cohort (n = 1101), a population containing the TCGA-Cell cohort (n = 817) 34 , Signature #2 associates with decreases in both OS and DFS within the follow-up period of 100 months (Fig. 4D,E). The TCGA cohorts contain an average of 18.8% lobular BC and has been used to profile the genomic and expression landscape of lobular BCs 34 . Of note, Signature #2 correlates with decreases in OS and DFS in lobular BC in TCGA (Fig. 4F,G) and the Curtis cohorts (Fig. 4H). Furthermore, Signature #2 independently predicts OS once adjusting for clinical features (Table 2). Collectively, these observations reveal associations of CYB5D2-derived 21-gene signature with decreases of OS and DFS in HER2-enriched, normal like molecular subtype, and lobular BC.

CYB5D2-associated DEGs affect pathways regulating cell proliferation and DNA metabolism.
To gain insight on the potential mechanisms contributing to the tumor suppression functions of CYB5D2 in BC, we analyzed pathways affected by the DEGs relative to CYB5D2 downregulation using the GAGE 44 and Reactome 45 packages in R as well as Ingenuity Pathway Analysis. Analysis of the 471 DEGs of both the TCGA-Cell (n = 817) and Metabric cohorts (common-DEGs, Table S3C) with the GAGE package using the KEGG gene sets identified three upregulated gene sets functioning in progesterone-regulated oocyte maturation, cell cycle, and oocyte meiosis (Table 4; Fig. S9; Tables S4A-C). Analysis of the entire 660 DEGs obtained from the TCGA-Cell (TCGA-DEGs) cohort recapitulated the enrichment (Table S4D), indicating all critical DEGs that occurred in TCGA-DEGs were present in Metabric-DEGs. Indeed, the same cell cycle KEGG gene set (hsa04110 cell cycle) was also enriched and upregulated in Metabric-DEGs (Table S4E). Two additional gene sets of DNA replication and ribosome biosynthesis were also enriched and upregulated in the Metabric-DEGs (Table S4E); pathways regulated by both gene sets belong to the core machinery of cell proliferation. In accordance with cell proliferation being critical for tumorigenesis, a cancer-related pathway was upregulated in common-DEGs (Table 4).
By using the GO gene set and the go.sets.hs datasets 44 , gene sets related to multiple aspects of cell cycle, mitotic phase, DNA metabolism, DNA replication, DNA repair, checkpoint activation and others were also significantly enriched in common-DEGs (Table 4; Table S5A). Similar gene sets of GO terms were enriched in TCGA-DEGs (Table S5B,C) and Metabric-DEGs (Table S5D,E).
Consistent with the above gene-set enrichment analyses, pathway enrichment determination using the Reactome package in R 45 yielded pathways regulated by the above enriched gene sets in common-DEGs, TCGA-DEGs, and Metabric-DEGs ( Fig. S10; Table S6A-C). The enriched pathways derived from common-DEGs  Table 3. Co-alteration of mutations with Signature #2 a . a CYB5D2 mRNA reduction at levels < −1.5 SD; mutations in co-alteration with Signature #2 were determined using the Curtis dataset (n = 817). b these mutations were co-altered with CYB5D2 downregulation (see Table S1). c log2-based ratio of percentage in altered group/percentage in unchanged group; positive and negative ratios are for co-occurrence and mutual exclusiveness, respectively.

Discussion
Evidence suggesting a role of CYB5D2 in BC suppression includes its reported tumor suppression in cervical cancer 11 and its chromosomal localization at 17p13.2; loss of which occurs frequently in BC 29 . While homodeletion of CYB5D2 can be detected in both the TCGA (8/1093 = 0.7%) and Metabric (4/2051 = 0.2%) datasets, the frequency is low; nonetheless, in both cohorts homodeletion in TP53 was not detected. However, mutation in TP53 occurs in 33% (659/1980) of the population in the Curtis and 34% (281/816) in the TCGA-Cell population. Although mutation in the CYB5D2 genes was undetectable in both, its expression is significantly reduced in vitro, in vivo and in primary BCs (Fig. 1); importantly, we have shown that CYB5D2 may regulate MCF7 cell viability or proliferation (Fig. 2). Although different genetic alterations are in play in the inactivation of p53 and CYB5D2, both events occur in a strong concordance; no less than 80% of tumors with CYB5D2 reduction < −1.5 SD

Gene sets b
Set size c p-value q-value  www.nature.com/scientificreports www.nature.com/scientificreports/ contain TP53 mutations (Table S1), and 85% of TP53 mutations are in tumors positive for CYB5D2-derived 21-gene signature (Table 3). This concordance suggests collaboration between p53 inactivation and CYB5D2 downregulation in BC tumorigenesis.
Significant co-occurrence of RB1 with CYB5D2 reduction (Table S1) and particularly Signature #2 in which 79.6% (39/49) of RB1 mutations were detected in Signature #2-positive tumors (Table 3) agrees well with the enrichment of CYB5D2 reduction-associated DEGs in cell cycle progression, mitotic phase events, and cell proliferation. On the other hand, the mutual exclusiveness of CYB5D2 downregulation with other major oncogenic events of PIK3CA, MAP3K1, GATA3, and CDH1 (Table 3) implies that CYB5D2 downregulation is a component of the oncogenic processes involving these oncogenic factors. While the association with above genomic alterations suggest mechanisms contributing to CYB5D2-derived tumor suppression in breast cancer, the potential mechanisms underlying the co-occurrence and mutual exclusion of CYB5D2 expression with the above gene mutations are likely complex. DEGs accompanied with CYB5D2 downregulation affect pathways regulating DNA repair (Table 4). This suggests a scenario that downregulation of CYB5D2 leads to genome instability, which facilitates gene mutations. However, this will not explain the selectivity of CYB5D2 downregulation with mutations in tumor suppressor and oncogenes as aforementioned above. With the current knowledge, we suggest that CYB5D2 downregulation selectively associates with genomic alterations via CYB5D2-derived tumor suppression activity. This possibility is consistent with CYB5D2 inducing apoptosis in MCF7 cells (Fig. 2C).
CYB5D2 reduction is associated with a large number of DEGs in primary BC (Table S3). While the underlying mechanisms affecting the gene expression have not been studied, it is possible that CYB5D2 indirectly alters gene expression though its tumor suppression function. Structural wise, CYB5D2 does not have motifs known to directly modulate gene expression 11,16 . Stable expression of CYB5D2 in HeLa cells did not substantially affect gene expression 16 .
Currently, there are approximate 140 driver genes functioning in 12 signaling pathways involving PI3K, MAPK, cell cycle, and DNA damage regulation 48 . These pathways are enriched in CYB5D2 downregulation-associated DEGs, supporting CYB5D2 as a novel tumor suppressor in BC. Nonetheless, it is important to directly examine this notion and to investigate the underlying mechanisms. Currently, we are studying these avenues. However, our observed correlation of CYB5D2 with BC progression and the results generated in silico using more than 3,000 patients with BC from two of the most comprehensive BC datasets provide a strong basis to explore CYB5D2-derived novel tumor supressing activities.
It appears that Signature #1 and #2 exhibit a more robust correlation with OS shortening in the Curtis population (Figs 3, 4A-C) compared to the TCGA-Cell cohort (Fig. 4D-H). A potential factor for these observations might be attributable to the difference in cohort composition (Table S7). While additional research is certainly required to determine the biomarker potential of the 21-gene signature, its predictive value in multiple intrinsic subtypes of BC (Figs 4A-C, S8) is not only appealing but also in accordance with its extensive overlap with TP53 mutations, an event that occurs in all intrinsic subtypes of BC. While evidence supports the unique biomarker value of CYB5D2-associated 21-gene signature, its clinical potential needs to be further tested both retrospectively and prospectively in future.
There are several promising multigene signatures available commercially to assess disease recurrence for newly diagnosed patients with different BC types, including Oncotype DX, MammaPrint, EndoPredict, and Prosigna 49,50 . Our signatures were constructed to predict OS and thus could be used together with these multigene panels to improve decision making and patient management.

Data Availability
We are committed to have all materials and data available to the research community upon publication.