Cytochrome b561 Serves as a Potential Prognostic Biomarker and Target for Breast Cancer.

PURPOSE: Cytochrome b561 (CYB561) is a transmembrane protein and participates in ascorbate recycling and iron homeostasis. However, its role in breast cancer remains unclear. PATIENTS AND METHODS: In this study, we explored the expression pattern and prognostic value of CYB561 in breast cancer through The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), PrognoScan and Kaplan-Meier Plotter and confirmed its mRNA expression in human breast cell lines. LinkedOmics, Metascape and Gene Expression Profiling Interactive Analysis (GEPIA2) databases were applied to investigate the co-expression genes and construct microRNA (miRNA) networks associated with CYB561. The correlations between CYB561 and immune infiltration cells and genes were also illustrated.
RESULTS: The CYB561 expression was upregulated in breast cancer tissues and cell lines and significantly correlated with the clinical features of breast cancer patients. High CYB561 expression was associated with poor survival and was an independent risk factor for overall and disease-specific survival. Functional enrichment analysis showed that CYB561 and its co-expressed genes were mainly enriched in lipid biosynthetic process, Wnt signaling pathway, Hippo signaling pathway, etc. The miRNA network analysis suggested that hsa-miR-497 was negatively correlated with CYB561 expression and was predicted to direct target CYB561. CYB561 expression was positively correlated with infiltrating levels of CD4+ T cells, neutrophils and dendritic cells in breast cancer. Subsequent analysis found that B cells could predict the outcome of breast cancer. Also, CYB561 showed strong correlations with diverse immune marker sets in breast cancer.
CONCLUSION: CYB561 may serve as a potential prognostic biomarker and target for breast cancer. Our findings laid foundation for future research on molecular mechanisms of CYB561 in breast cancer.

Introduction

Breast cancer is the most common carcinoma and has surpassed lung cancer ranking first in the incidence of malignant among women worldwide.1 Although comprehensive treatment including surgery, chemotherapy, radiotherapy, endocrine therapy, target therapy and immunotherapy, have significantly improved the prognosis of breast cancer patients, breast cancer remains the leading cause of cancer death among women due to the lack of specific biomarkers for early diagnosis.1 Moreover, there are still challenges to understanding the biological behavior of breast cancer and formulating personalized treatments. Biomarkers are particularly useful in identification, diagnosis, determining prognosis and guiding treatment in breast cancer.2 Thus, digging for novel biomarkers is important to improve early diagnosis and develop better treatment strategies and ultimately improve breast cancer patients’ outcome.

The development of next-generation sequencing technology has provided us with a large amount of cancer multi-omics data. While the analysis and integration of these raw data can only be carried out with the application of bioinformatics tools. Currently, bioinformatics plays a vital role in both biological sciences and medical research. In cancer research, through bioinformatics analysis, we can obtain the molecular mechanisms and pathways of cancer onset and metastasis, identify new potential therapeutic targets, and study drug resistance mechanisms.

Cytochrome b561 (CYB561/CYB561A1), also known as chromaffin granule cytochrome b561, is a member of the CYB561 protein family, whose gene is located at 17q23.3. CYB561 protein family are transmembrane proteins consisting of six transmembrane helices and two b-type hemes on each side of the membrane. They act as monodehydroascorbate reductases and ferric reductases, participating in ascorbate recycling and iron homeostasis.3 Apart from CYB561, members of the CYB561 protein family in mammals include duodenal CYB561 (DCYTB/CYBRD1/CYB561A2), lysosomal CYB561 (CYB561A3), stromal cell-derived receptor 2 (SDR2/FRRS1) and 101F6 (CYB561D2/TSP10), which is related to suppression of tumor cell growth.4

Ascorbate and iron are important cellular components and are needed as co-factors in numerous metabolic reactions. Accumulating evidence indicates that ascorbate promotes iron uptake and regulates its metabolism.4 Unlike normal cells, the growth of cancer cells strongly depends on the trace element iron. Previous studies have implicated diverse roles of iron in cancer, such as tumorigenesis, tumor development and metastasis.5 These findings led us to conjecture that CYB561 may play a role in the biological process of breast cancer.

Loss of function mutation of CYB561 was identified in patients with lifelong disabling orthostatic hypotension which suggested its pivotal role in sympathetic function and cardiovascular regulation.6 In contrast, CYB561 was less studied in cancer. Few studies suggested that CYB561 was highly expressed in several cancer cell lines and may play an important role in cancer progression, invasion and metastasis.7,8 However, the biological function and mechanism of CYB561 in breast cancer remain unclear. In the present study, we combined the data of The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), PrognoScan, Kaplan–Meier plotter and breast cell lines to comprehensively analyze the CYB561 expression and correlation with prognosis of breast cancer patients. Moreover, the multidimensional analysis was used to construct the co-expression and functional network associated with CYB561 in breast cancer, and the correlation of CYB561 with tumor-infiltrating immune cells and immune signatures via Tumor Immune Estimation Resource (TIMER) and Gene Expression Profiling Interactive Analysis (GEPIA2). Our findings may reveal potential therapeutic targets and strategies for breast cancer diagnosis and treatment.

Materials and Methods

Patient Datasets

UCSC Xena is a web-based tool that assembles multi-omics data from public databases and associated annotations, such as TCGA and GTEx.9 TCGA is a public database providing researchers open access to the multiple cancer genomic profiles for analyses and publications.10 Breast cancer patients’ gene and miRNA expression profiles as well as their clinical data were downloaded from the UCSC Xena and cBioPortal websites ( and ), respectively.

Two independent validation datasets GSE109169 and GSE1456 were also downloaded at the GEO website ().11–13 The dataset GSE109169 contains mRNA expression data by microarray from 25 paired breast cancer tissues and adjacent normal tissues, which was utilized to validate CYB561 expression patterns in breast and normal tissues. The dataset GSE1456, which provides microarray data of 318 breast cancer patients with follow-up information, was analyzed to confirm the prognostic value of CYB561 expression.

Cancer Cell Line Encyclopedia (CCLE) Database

The CCLE database () was used to investigate and visualize the relative mRNA levels of CYB561 in a variety of breast cancer cell lines.14 The RNA sequencing data and annotation data of cell lines obtained from CCLE website were utilized to plot a heatmap using the OmicShare tools ().

The Human Protein Atlas (HPA)

The HPA database () consists of six separate parts and uses the integration of various omics technologies to provide human proteins in cells, tissues and organs.15–17 The Tissue Atlas and The Pathology Atlas include immunohistochemistry (IHC) data, using tissue microarray-based analysis to analyze 44 different normal tissue types and to perform proteome analysis of 17 major cancer types. The staining intensity, quantity, location and patient information of patients with various cancer types can be obtained online. Protein expression of CYB561 between human breast cancer and normal tissue was compared by IHC images from HPA database.

PrognoScan Database and Kaplan–Meier Plotter Analysis

PrognoScan database () provides meta-analysis of the prognostic value of genes in a large number of publicly available cancer microarray datasets.18 Kaplan–Meier Plotter () is a web application utilized to evaluate the impact of 54k genes on the survival rate of 21 cancer types, including breast, ovarian, lung and stomach cancer.19 The correlation between CYB561 expression and breast cancer survival was analyzed in PrognoScan and Kaplan–Meier Plotter. The analyses from PrognoScan were visualized by GraphPad software.

LinkedOmics Database Analysis

LinkedOmics database () is an open portal for analyzing multi-omics data from all 32 TCGA cancer types.20 CYB561 co-expression genes and associated miRNAs were analyzed statistically using Pearson correlation test and the results were displayed in the form of volcano plot, heat maps, or scatter plots.

GEPIA2 Database Analysis

The GEPIA2 database () is a network server used to analyze the RNA sequencing expression data of tumor and normal samples from the TCGA and GTEx projects.21 GEPIA2 was used to plot correlation analysis and survival curves of co-expression genes, and survival heat maps of top 50 positively and negatively correlated immune marker genes.

Metascape Database Analysis

Metascape () is an open online tool for gene annotation and analysis.22 In this study, Metascape was used for the Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Terms with P-value <0.01, minimum count of 3, and enrichment factor of >1.5 were considered as significant. The following databases were used for protein–protein interaction enrichment analysis: STRING,23 BioGrid,24 OmniPath,25 InWeb_IM.26 If the generated network contains 3 and 500 proteins, the Molecular Complex Detection (MCODE) algorithm27 is used to identify densely connected network components, and the three terms with the highest p-score are retained as functional descriptions of the pathway or process enrichment components.

Immune Infiltration Analysis

The TIMER database () is an integrated resource for systematically analyzing the immune infiltration of different cancer types.28,29 TISIDB () is a comprehensive repository portal for tumor-immune system interactions.30 Gene signatures of 28 tumor-infiltrating lymphocytes, chemokines, receptors, immunomodulators and major histocompatibility complex (MHC)31 were downloaded from the database. TIMER was used to evaluate the correlation between CYB561 expression and the abundance of infiltration of six types of immune cells and immune signatures in breast cancer with partial Spearman correlation corrected for tumor purity. Kaplan–Meier analysis was performed to evaluate the prognostic value of each immune infiltrate in breast cancer. Cox analysis was used to assess how CYB561 and six types of immune cells affect breast cancer outcomes.

Cell Lines

The human breast epithelial cell line MCF-10A and cancer cell lines (MCF-7, MDA-MB-231, MDA-MB-468, SK-BR-3 and HCC70) were purchased from Procell Life Science&Technology Co., Ltd. MCF-10A was grown in MCF-10A cell special medium (CM-0525, Procell, Wuhan, China). Those breast cancer cells were grown in DMEM medium (Gibco, Grand Island, NY, USA) containing 10% fetal bovine serum, 100 ug/mL streptomycin and 100 U/mL penicillin at 37°C with an atmosphere of 5% CO2 and 95% humidity.

RNA Isolation and Quantitative Polymerase Chain Reaction (RT-qPCR)

Total RNA was collected from MCF-10A and five breast cancer cells using TRIzol® Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. Complementary DNA was generated from each 500ng RNA sample using Fast All-in-one RT kit (ESscience Biotech, Shanghai, China). Real-time qPCR analysis was carried out on the ABI7600 Prism system using the 2x Super SYBR Green qPCR Master Mix kit (ESscience Biotech, Shanghai, China). Melt curve analysis was always performed to calculate the product melting temperature. The 2−∆∆Ct method was used for the quantification of CYB561 mRNA, with GAPDH mRNA tested as an internal control. The following primer sequences were used: CYB561, sense strand 5ʹ- CCTGCTGGTTTACCGTGTCTTC −3ʹ, antisense strand 5ʹ- CTTCCTGTGGTAGTCGAACACC −3ʹ; GAPDH, sense strand 5ʹ- ATGTTCGTCATGGGTGTGAAC-3ʹ, antisense strand 5ʹ- ATGGACTGTGGTCATGAGTCC −3ʹ. We performed the RT-qPCR assays in duplicate at least three independent experiments.

Statistical Analysis

The analyses of differential mRNA expression of CYB561 and miRNAs expressions in tumor and normal tissues were examined by the Wilcoxon test, including unpaired or paired test. The associations of CYB561 expression with clinical characteristics and breast cancer patient prognosis were examined by a non-parametric test and Cox regression analyses using SPSS software. Survival analyses were conducted between CYB561High and CYB561low groups in the TCGA-BRCA cohort through Kaplan–Meier analysis with the Log rank test, and were visualized by GraphPad software. The correlation between CYB561 expression and hsa-miR-497 expression was assessed by GraphPad software. Data were considered statistically significant at p < 0.05 derived from two-tailed tests. The Receiver operating characteristic (ROC) curve analysis was generated to identify the ability of CYB561 to predict prognosis event by R (3.6.3 version).

Results

Overexpression of CYB561 in Breast Cancer Tissues and Cell Lines

To identify the role of CYB561 in breast cancer, we initially explored CYB561 mRNA expression based on the RNA sequencing expression data from the TCGA-BRCA cohort. As shown in Figure 1A and B, the unpaired and paired tests both indicated that the mRNA expression of CYB561 in breast tumor tissue was significantly elevated (p < 0.001). We also performed transcriptional analysis of another breast cancer cohort of patients from the GEO dataset (GSE109169), consisting of 25 paired breast cancer specimens, and demonstrated that CYB561 was significantly enhanced in tumor tissues compared with matched normal tissues (Figure 1C). We further analyzed CYB561 protein expression in HPA database, discovering it was highly expressed in tumor tissues (Figure 1D). Besides, we examined CYB561 expression among 60 breast cancer cell lines detected in the CCLE database, which provides public access to genomic data, analysis and visualization for over 1100 cell lines. Consistent with the results in breast cancer patients, the expressions of CYB561 in breast cancer cell lines were the highest among all tumor cell lines (Figure 1E). Moreover, to validate CYB561 mRNA expression in breast cancer cell lines, we performed RT-qPCR in human breast epithelial cell line and five breast cancer cell lines. We found dramatically increased CYB561 mRNA expression in breast cancer cell lines than normal breast epithelial cell line (Figure 1F). Taken together, these results suggested that CYB561 was highly expressed at transcriptional and protein levels in breast cancer tissues and cell lines compared with normal tissues and cell line.

Figure 1

CYB561 was significantly overexpressed in breast cancer. (A) CYB561 mRNA expression comparison between unpaired normal and tumor tissues in the TCGA-BRCA cohort. (B) CYB561 mRNA expression comparison between paired normal and tumor tissues in the TCGA-BRCA cohort. (C) CYB561 expression was elevated in breast cancer compared to paired normal tissues in GEO dataset (GSE109169). (D) CYB561 protein expression comparison between normal and tumor tissues obtained from the HPA database. (E) The expression of CYB561 in human cancer cell lines, analyzing by the CCLE dataset. (F) The mRNA expression of CYB561 in human breast epithelial cell line MCF-10A and five breast cancer cell lines. ***p value <0.001, ****p value <0.0001.

Different Clinical Characteristics Between CYB561high and CYB561low Groups

To download the raw data from TCGA-BRCA database, we analyzed the correlation between mRNA expression data and related clinical information. In total, 1009 female breast cancer patients with complete follow-up information were analyzed. Patients were divided into two groups (CYB561high and CYB561low) according to the median value of CYB561 expression level for analysis of the relationship between clinical features and CYB561 mRNA expression (Table 1). According to Chi-square tests, high CYB561 mRNA expression was highly associated with age (p < 0.047), histological type (p < 0.001), molecular subtype (p < 0.001), T classification (p < 0.012), vital status (p = 0.023) and disease-specific survival status (p = 0.024). There were no significant differences between CYB561high and CYB561low groups in N classification, M classification, stage, race, radiation therapy, neoadjuvant treatment, progression-free status and disease-free status (p > 0.05).

Table 1

Correlation of CYB561 Expression with Clinicopathologic Characteristics in Breast Cancer

Characteristics	Variable	Numbers	CYB561 mRNA Expression		χ2	p value
			High n (%)	Low n (%)
Age	<60	548	258 (51.19)	290 (57.43)	3.952	0.047
	≥60	461	246 (48.81)	215 (42.57)
Histological type	Infiltrating Ductal Carcinoma	746	398 (78.97)	348 (68.91)	24.903	0.000
	Infiltrating Lobular Carcinoma	198	73 (14.48)	125 (24.75)
	Medullary Carcinoma	5	0 (0)	5 (0.99)
	Mixed Histology	29	14 (2.78)	15 (2.97)
	Mucinous Carcinoma	2	2 (0.40)	0 (0)
	Other	29	17 (3.37)	12 (2.38)
Molecular subtype	Basal	158	35 (6.94)	123 (24.36)	102.027	0.000
	Her2	72	54 (10.71)	18 (3.56)
	Luminal A	482	236 (46.83)	246 (48.71)
	Luminal B	187	132 (26.19)	55 (10.89)
	Normal	32	11 (2.18)	21 (4.16)
	Not available	78	36	42 (8.32)
T classification	T1	264	129 (25.60)	135 (26.73)	12.782	0.012
	T2	586	293 (58.13)	293 (58.02)
	T3	123	55 (10.91)	68 (13.47)
	T4	33	26 (5.16)	7 (1.39)
	TX	3	1 (0.20)	2 (0.40)
N classification	N0	464	230 (45.63)	234 (46.34)	4.253	0.373
	N1	336	161 (31.94)	175 (34.65)
	N2	117	61 (12.10)	56 (11.09)
	N3	72	38 (7.54)	34 (6.73)
	NX	20	14 (2.78)	6 (1.19)
M classification	M0	837	419 (83.13)	418 (82.77)	3.858	0.145
	M1	20	14 (2.78)	6 (1.19)
	MX	152	71 (14.09)	81 (16.04)
Stage	I	169	86 (17.06)	83 (16.44)	7.328	0.120
	II	573	273 (54.17)	300 (59.41)
	III	231	121 (24.01)	110 (21.78)
	IV	18	13 (2.58)	5 (0.99)
	X	13	9 (1.79)	4 (0.79)
	Not available	5	2 (0.40)	3 (0.59)
Race	American Indian or Alaska Native	1	1 (0.20)	0 (0)	4.930	0.177
	Black or African American	170	85 (16.87)	85 (16.83)
	White	706	331 (65.67)	375 (74.26)
	Asian	53	32 (6.35)	21 (4.16)
	Not available	79	55 (10.91)	24 (4.75)
Vital status	Deceased	143	84 (16.67)	59 (11.68)	5.150	0.023
	Living	866	420 (83.33)	446 (88.32)
Radiation therapy	NO	405	202 (40.08)	203 (40.20)	0.078	0.780
	YES	525	257 (51.00)	268 (53.07)
	Not available	79	45 (8.93)	34 (6.73)
Neoadjuvant treatment	NO	1002	499 (99.01)	503 (99.60)	2.683	0.217
	YES	6	5 (0.99)	1 (0.20)
	Not available	1	0 (0)	1 (0.20)
Progression free status	Censored	872	433 (85.91)	439 (86.93)	0.223	0.637
	Progression	137	71 (14.09)	66 (13.07)
Disease free status	Recurred/progressed	81	38 (7.54)	43 (8.51)	0.187	0.665
	Disease free	797	394 (78.17)	403 (79.80)
	Not available	131	72 (14.29)	59 (11.68)
Disease-specific Survival status	Alive or dead tumor free	914	447 (88.69)	467 (92.48)	5.125	0.024
	Dead with tumor	77	48 (9.52)	29 (5.74)
	Not available	18	9 (1.79)	9 (1.78)

Note: Bold values indicate p value <0.05.

CYB561high is Associated with Adverse Outcomes in Breast Cancer

In order to confirm the correlation between CYB561 expression and patient prognosis, Kaplan–Meier survival curve was used to evaluate and compare the survival differences between patients with high and low CYB561 expression (grouped by median). In the TCGA-BRCA cohort, the overall survival (OS) of the CYB561high group was significantly shortened, and the median OS of CYB561high vs that of CYB561low group was 122.8 months vs 216.8 months (Log rank test, p = 0.0004, Figure 2A). The unfavorable disease-specific survival (DSS) rate in CYB561high group was also significantly lower than that in CYB561low group (Log rank test, p = 0.0019, Figure 2B). In addition, we examined the PrognoScan and Kaplan–Meier Plotter and found that high CYB561 expression was associated with poor OS and DSS (Figure 2C–E).

Figure 2

CYB561 is associated with survival outcome. (A and B) Survival analyses of CYB561 in the TCGA-BRCA cohort by Kaplan–Meier estimator with a Log rank test. (C and D) Survival analyses of CYB561 by Kaplan–Meier estimator with Log rank test obtained from PrognoScan web tool. (E) Survival analysis of CYB561 by Kaplan–Meier estimator with Log rank test obtained from the Kaplan–Meier plotter web tool. (F) ROC curve was generated to validate the ability of CYB561 for predicting prognosis. The AUC index was 0.873. Survival differences are compared between patients with high and low expression of CYB561 (grouped by median). p<0.05 was used to assess differences.

High Expression of CYB561 is a Potential Independent Risk Factor

In addition, Cox regression was performed to assess the prognostic value of CYB561 expression as an independent predictor of survival in breast cancer. Univariate Cox analyses in OS illustrated that age, T classification, N classification, tumor stage and CYB561 expression were the potential risk roles in breast cancer. More importantly, the multivariate Cox analyses validated the critical value of age, N classification, tumor stage and CYB561, proving that they can predict tumor prognosis independently of other factors in OS (Table 2). Additionally, the Cox analysis based on DSS showed that age, tumor stage and CYB561 could serve as independent predictors of shorter DSS in breast cancer patients (Table 2). Receiver operating characteristic (ROC) curve analysis was generated to identify the ability of CYB561expression to predict prognosis event. The area under ROC curve (AUC) of CYB561 was 0.873, demonstrating that CYB561 has relatively high sensitivity and specificity for breast cancer prognosis prediction (Figure 2F).

Table 2

Cox Regression Analyses of the Correlation of CYB561 Expression and Important Clinical Characteristics with OS and DSS Among TCGA-BRCA Cohort Patients

Parameters	OS				DSS
	Univariate Analysis		Multivariate Analysis		Univariate Analysis		Multivariate Analysis
	HR (95% CI)	P value	HR (95% CI)	p value	HR (95% CI)	p value	HR (95% CI)	p value
Age	2.000 (1.435–2.788)	0.000	2.257 (1.599–3.187)	0.000	1.471 (0.940–2.303)	0.091	1.631 (1.036–2.567)	0.035
Histological type	1.033 (0.900–1.186)	0.642			0.991 (0.810–1.212)	0.930
Molecular subtype	1.115 (0.942–1.319)	0.642			0.967 (0.768–1.217)	0.772
Race	0.213 (0.540–1.148)	0.207	0.825 (0.555–1.228)	0.343	0.793 (0.479–1.314)	0.368
T classification	1.252 (1.044–1.501)	0.124	0.852 (0.677–1.072)	0.172	1.370 (1.071–1.752)	0.012	0.908 (0.698–1.180)	0.471
N classification	1.552 (1.360–1.771)	0.015	1.313 (1.081–1.595)	0.006	1.731 (1.452–2.062)	0.000	1.304 (0.981–1.734)	0.067
M classification	1.228 (0.965–1.563)	0.000	0.759 (0.553–1.043)	0.090	1.209 (0.874–1.672)	0.252
Stage	1.653 (1.411–1.936)	0.096	1.575 (1.184–2.095)	0.002	1.949 (1.590–2.388)	0.000	1.584 (1.098–2.286)	0.014
CYB561	1.817 (1.300–2.539)	0.000	1.722 (1.217–2.437)	0.002	2.052 (1.292–3.261)	0.002	1.904 (1.195–3.032)	0.007

Note: Bold values indicate p value <0.05.

Abbreviations: OS, overall survival; DSS, disease-specific survival; BRCA, breast invasive adenocarcinoma; TCGA, The Cancer Genome Atlas; HR, hazard ratio; CI, confidence interval.

CYB561 Co-Expression Networks in Breast Cancer

To gain insights into the biological function of CYB561 in breast cancer, the function module of LinkedOmics was used to examine CYB561 co-expression mode in breast cancer. As shown in Figure 3A, 4101 genes (dark red dots) showed a significant positive correlation with CYB561, while 7293 genes (dark green dots) showed a significant negative correlation (false discovery rate, FDR < 0.01). The top 50 significant genes that were positively and negatively correlated with CYB561 were shown in the heat map (Figure 3B and C). See for a complete description of co-expressed genes.

Figure 3

CYB561 co-expression genes in breast cancer. (A) Highly correlated genes with CYB561 identified by the Pearson test in BRCA. Red and green dots represent positively and negatively significantly correlated genes with CYB561, respectively. (B and C) Expression heatmap of the top 50 positively and negatively correlated genes with CYB561 in BRCA. (D–F) Kaplan–Meier analyses of the expression of SRXN1, DLG3, and MRPS23 on overall survival in TCGA-BRCA cohort.

Several positively correlated genes including FTSJ3,32 LLGL2,33 LRRC59,34 PSMD12,35 C17orf37,36 TRIM37,37 MPRS23,38 SRXN1,39 and DLG340 were reported with pro-breast cancer effects. Among them, Kaplan–Meier survival analysis showed that patients with relatively higher SRXN1, DLG3, and MRPS23 expression had shorter OS times than patients with lower gene expression levels in the TCGA-BRCA cohort (Figure 3D, p = 0.014; Figure 3E, p = 0.0018; Figure 3F, p = 0.024; respectively).

The functions of CYB561 and its top 500 co-expression genes were predicted by analyzing GO and KEGG in Metascape. The top 20 GO enrichment items and top 20 KEGG pathways for CYB561 and its co-expression genes are shown in Figure 4A–D respectively. Among these pathways, the Hippo signaling pathway, peroxisome proliferator-activated receptor (PPAR) signaling pathway, MAPK signaling pathway and TGF-β signaling pathway were found to be related to multiple tumor development and were involved in breast cancer tumorigenesis and pathogenesis.41,42 Moreover, a Metascape protein–protein interaction (PPI) enrichment analysis was performed to better understand the relationship between CYB561 and breast cancer (Figure 4E). The results showed that biological function was mainly related to PPAR signaling pathway, Wnt signaling pathway, Epstein-Barr virus infection, human papillomavirus infection, pathways in cancer, regulation of actin cytoskeleton and proteasome (Figure 4F).

Figure 4

The enrichment analysis of CYB561 and its co-expression genes in BRCA (Metascape). (A) GO analysis of top 500 co-expression genes conducted using online website of Metascape. Heatmap of GO enriched terms colored by p values. (B) Network of GO enriched terms colored by p value, where terms containing more genes tend to have a more significant p value. (C) Heatmap of KEGG enriched terms colored by p values. (D) Network of KEGG enriched terms colored by p value, where terms containing more genes tend to have a more significant p value. (E) Protein–protein interaction network. (F) Independent functional enrichment analysis of three MCODE components.

Associations Between Genome-Wide miRNA Profiles and CYB561 Expression

A genome-wide analysis of miRNA-sequencing data was carried out to identify miRNA profiles showing significant correlation to CYB561 expression using online website of LinkedOmics. As shown in Figure 5A and , 823 miRNAs were found, after filtering profiles with FDR <0.05, 93 miRNAs (dark red dots) were shown significant positive correlations with CYB561, whereas 192 miRNAs (dark green dots) were shown significant negative correlations.

Figure 5

Genome-wide microRNA profiles associated with CYB561 expression in BRCA. (A) Volcano plot of microRNAs positively and negatively correlated with CYB561 expression. (B) Venn results of microRNAs which could target CYB561 predicted by TargetScan (), mirDIP (), miRWalk (), and miRDB (). (C) The has-mir-497 expression comparison between normal and tumor tissues in the TCGA-BRCA cohort (unpaired Wilcoxon test). (D) Correlation analysis of the expression of CYB561 and has-mir-497 in the TCGA-BRCA cohort (Pearson correlation test).

Among the positively correlated miRNAs, hsa-miR-425-5p,43 hsa-miR-454,44 and hsa-miR-126645 had high expression levels in breast cancer and were associated with tumor progression and poor outcome. Increased expression of hsa-miR-375 resulted in loss of cellular organization and acquisition of a hyperplastic phenotype, thus contributes to lobular neoplastic progression.46 Moreover, has-miR-191 is a direct estrogen receptor target and acts as an oncogenic regulator in human breast cancer.47

Among the negatively correlated miRNAs, high expression of hsa-miR-101, hsa-miR-217 and hsa-miR-497 could inhibit breast cancer proliferation.48–50 High expression of hsa-miR-493, hsa-miR-497 and hsa-miR-588 also conferred a better prognosis.50–52 More interestingly, hsa-miR-497 was predicted to direct target CYB561 by online web tools and showed negative correlations with CYB561 (Figure 5B, ). We validated that hsa-miR-497 was low expressed in breast tumor tissues, and was negatively correlated between CYB561 expression in breast cancer samples from TCGA (Figure 5C and D).

Correlation Analysis Between CYB561 Expression and Six Main Infiltrating Immune Cells

We investigated whether the expression of CYB561 was related to the immune infiltration levels in breast cancer from the TIMER database. The results show that the expression of CYB561 was significantly correlated with tumor purity, as well as the infiltrating levels of CD4+ T cells, neutrophils and dendritic cells (Figure 6A).

Figure 6

Correlation analyses between CYB561 expression and six types of infiltrating immune cells and top 50 positively and negatively correlated marker genes in BRCA. (A) Correlation between CYB561 expression and six immune infiltrating cells obtained from TIMER (purity-corrected Spearman test). (B) Overall survival curve of each of the six immune infiltrating cells generated by TIMER’s Kaplan–Meier estimator. Survival differences are compared between patients with high and low infiltrating of each kind of immune cells (grouped by median). (C) Survival heatmaps of the top 50 genes positively and negatively correlated with CYB561 in BRCA. The red and blue blocks indicate higher and lower risks, respectively. Boxed rectangles indicate significant unfavorable and favorable results in prognostic analyses (p < 0.05).

Moreover, we evaluated the prognostic value of each of the six types of immune cells via Kaplan–Meier analysis, finding B cells (p = 0.046 in Log rank test) can predict the outcome of breast cancer (Figure 6B). At last, Cox proportional hazard models were applied to assess the impacts of CYB561 expression and the six types of immune cells on the overall survival of breast cancer. CYB561 (HR = 1.21, 95% CI = 1.016–1.441, p = 0.033) and B cells (HR = 0.067, 95% CI = 0.006–0.74, p = 0.027) showed significant risk in univariate analyses, and multivariate analyses indicated that CYB561 (HR = 1.225, 95% CI = 1.017–1.475, p = 0.033) can predict tumor outcomes in the presence of varying immune cells (Table 3).

Table 3

Cox Analyses of the Correlation Between CYB561 Expression and Six Types of Immune Cells and Prognosis in Patients with Breast Cancer

Parameters	Univariate Analysis			Multivariate Analysis
	Coef	HR (95% CI)	p value	Coef	HR (95% CI)	p value
B cell	−2.705	0.067 (0.006–0.74)	0.027	−1.834	0.160 (0.003–9.864)	0.383
CD8+ T cell	−0.497	0.608 (0.184–2.008)	0.414	−0.672	0.511 (0.050–5.173)	0.570
CD4+ T cell	−0.731	0.482 (0.092–2.509)	0.386	−0.065	0.937 (0.029–30.787)	0.971
Macrophage	1.257	3.515 (0.572–21.623)	0.175	2.142	8.518 (0.687–105.673)	0.095
Neutrophil	−0.769	0.463 (0.069–3.101)	0.428	1.846	6.335 (0.046–873.831)	0.463
Dendritic	−0.285	0.752 (0.417–1.358)	0.345	−0.582	0.559 (0.080–3.922)	0.558
CYB561	0.191	1.21 (1.016–1.441)	0.033	0.203	1.225 (1.017–1.475)	0.033

Note: Bold values indicate p value <0.05.

Abbreviations: coef, regression coefficient; HR, hazard ratio; CI, confidence interval.

Correlation Between CYB561 Expression and Immune Signatures

At last, to broaden the understanding of CYB561 and immune gene crosstalk, we analyzed the correlations between CYB561 expression and various immune signatures, including immune marker genes of 28 tumor-infiltrating lymphocytes (TILs), immune inhibitory or stimulatory, cytokine-related, cancer-testis antigen and MHC.

After the correlation adjustment by tumor purity, the analysis revealed that the expression level of CYB561 in breast cancer was significantly correlated with 66.70% (615/922) immune marker genes, which included 34.47% (212/615) positively related and 65.53% (403/615) negatively related immune marker genes (). In activated CD8 T cell, CYB561 is highly correlated with marker gene PTRH2, which has been identified as a potential metastasis suppressor via its regulation of anoikis and epithelial-to-mesenchymal transition in breast cancer.53 For activated CD4 T cell, CYB561 is significantly correlated with BRIP1. Actually, BRIP1is a moderate risk gene for breast cancer, and has been proved to promote breast cancer cell invasion recently.54 Activated dendritic cell markers such as BCL2L1, SLC25A37, ATP5B, SPCS3, RAB1A, CCNA1 and CD207, etc., were also shown significant correlations with CYB561 expression.

Besides, we found significant correlation between CYB561 and immune marker genes of regulatory T cell and myeloid-derived suppressor cell, such as PELO, LIPA, MARCO, STAB1, CD14, CXCR4 and PARVG. Interestingly, LIPA has a critical role in regulating MDSCs’ ability to directly stimulate cancer cell proliferation and overcome immune rejection of cancer metastasis in allogeneic mice through modulation of the mTOR pathway.55 Hypoxia-induced CXCR4 up-regulation correlates with Treg recruitment and suppresses the anti-tumor immune response in basal-like breast cancer.56

As for immunoinhibitory genes, results showed the expression levels of PVRL2 and KDR have positive correlations with CYB561 expression, while VTCN1, PDCD1, CSF1R, IDO1, LAG3, CTLA4, CD244, CD160, PDCD1LG2, LGALS9, BTLA and CD96 have negatively correlations with CYB561 expression. In immunostimulatory genes, the top 5 positively correlated genes with CYB561 expression were TNFSF13, MICB, TNFSF4, CD276 and ENTPD1, while the top 5 negative marker genes were C10orf54, IL6, TNFRSF25, TNFRSF13C and TNFRSF8. Specifically, we showed CXCL3, CXCL2, CX3CL1, CXCL16, CXCL5, CXCL1, CCL20 and CXCL6 were significantly correlated with CYB561 expression (p < 0.00e-10).

Overall, the top 5 positively correlated marker genes with CYB561 were PTRH2, TMBIM6, GPRC5C, NOL11 and NUCB2. And the top 5 negatively correlated marker genes with CYB561 were RPS7, STAC, HAPLN3, IL34 and UBA52. Survival map analysis suggested the high risk of CYB561 positively correlated marker genes and the low risk of CYB561 negatively correlated marker genes (Figure 6C). Together, the results confirm that CYB561 is specifically correlated with immune infiltrating in breast cancer.

Discussion

Iron is a vital trace element and plays a major role in many crucial biological processes, in particular DNA synthesis, cell growth and proliferation.57 The biological versatility of iron owes to its ability to pass electrons by converting between ferric (Fe3+) and ferrous (Fe2+) oxidation states. Epidemiological studies have shown that increased body iron levels are associated with elevated risk of several types of cancer, including breast cancer.58 Iron acts as a tumor initiator and a growth factor in breast cancer.5 CYB561, which functions as a transmembrane electron transfer carrier and ferric reductase, is probably involved in iron metabolism. However, the biological function of CYB561 in breast cancer had remained to be elucidated.

In the present study, we demonstrated that CYB561 expression was remarkably increased in breast cancer tissues and cell lines, and high expression of CYB561 was significantly associated with age, histological type, molecular subtype, T classification and poor clinical outcome. Univariate and multivariate Cox analyses indicated that CYB561 might be a potential independent biomarker for breast cancer prognosis. Subsequently, we examined the co-expression genes and miRNA networks of CYB561. At last, we conducted a correlation analysis between CYB561 and immune infiltration or immune signatures, indicating that CYB561 was associated with the most of the immune marker genes. These results provide guidance for future research on CYB561.

Previous study has analyzed the expression patterns of the CYB561 transcripts in a variety of cancer lines and detected significant levels of CYB561 mRNA expression in cancer cell lines compared to peripheral blood leukocyte.7 Overexpression of CYB561 was also observed in highly metastatic and androgen refractory prostate cancer cells in comparison to androgen responsive and non-metastatic prostate cancer cells, which indicated CYB561 may be involved in prostate cancer progression and metastasis.8 On the contrary, a recent meta-analysis based on bioinformatics analysis identified that low mRNA expression of CYB561 was prognostic of poor outcome in ovarian serous carcinomas.59 Our results indicate that CYB561 was upregulated in breast cancer and was a potential prognostic marker, which is worthy of future clinical verification.

We further explored the molecular signatures associated with CYB561 in breast cancer to get better understanding of breast cancer biological behavior. For CYB561 co-expression network analysis, we found that CYB561 dysregulation was associated with lipid biosynthetic process, Wnt signaling pathway, proteoglycans in cancer and Hippo signaling pathway. Enhanced lipogenesis is an important metabolic hallmark of cancer cells. Elevated levels of enzymes responsible for fatty acid biosynthesis are correlated with poor prognosis in breast cancer patients.60 The Wnt signaling pathway is involved in numerous pathological processes, including carcinogenesis.61 Proteoglycans, hybrid molecules consisting of a protein core to which sulfated glycosaminoglycan chains are bound, are significant components of the extracellular matrix that are implicated in all phases of tumorigenesis.62 Moreover, for miRNAs, we found CYB561 expression was negatively correlated with several miRNAs such as hsa-miR-101,48 hsa-miR-217,49 hsa-miR-497,50 hsa-miR-49351 and hsa-miR-588,52 which were reported to be associated with inhibition of breast cancer or patients’ better prognosis by previous studies. Of these miRNAs, hsa-miR-497 was identified as predicted miRNA that could direct target CYB561. Obviously, further studies are needed to confirm the direct connection of CYB561 with hsa-miR-497.

The results of TIMER analyses showed that CYB561 expression levels had significant correlation with infiltrating levels of CD4+ T cells, neutrophils and dendritic cells. The subsequent Kaplan–Meier analysis found that B cells could predict the outcome of breast cancer. Cox analyses showed that CYB561 was a significant independent risk factor among all parameters. The association between CYB561 and immune signatures was further investigated, and results showed that CYB561 is specifically correlated with immune infiltrating in breast cancer. Our study suggests that CYB561 may be involved in the immune response to breast cancer development, resulting in a poor prognosis in patients with breast cancer.

Ferroptosis is a new defined form of regulated cell death coined in 2012, which is characterized by iron-dependent lipid peroxidation and subsequent cell membrane damage.63 Increased iron accumulation, free radical production, fatty acid supply and lipid peroxidation are essential for the occurrence of ferroptosis.64 Ferroptosis plays a dual role in tumorigenesis, including promoting and inhibiting tumors. The different effects rely on the release of damage-associated molecular patterns and the activation of immune response in the tumor microenvironment.64 Ferroptosis is becoming as a promising therapeutic method to treat cancer. In breast cancer, ferroptosis has been applied to overcome therapy resistance. CYB561, functions as a ferric reductase, is associated with iron homeostasis, lipid biosynthetic process and breast cancer immunity, and is very likely to be involved in ferroptosis pathway, making it a potential therapeutic target.

At present, how CYB561 affects the prognosis of breast cancer and the biological function of CYB561 in breast cancer is still in its infancy. In our study, we conducted a comprehensive and multi-dimension analysis of CYB561 by combining multiple databases to explore the role of CYB561 in breast cancer. Immunity plays a crucial role in the development of tumors. There is no research declaring how CYB561 affects breast cancer immunity. We not only expounded the relationship between CYB561 and six immune cells, but also discussed in detail the immune signature genes related to CYB561 that affect the prognosis. Such a research design for CYB561 and breast cancer has never been reported before, which is aiming to provide much more of a possibility and direction for future breast cancer research.

There were also some limitations in our study. First, although high mRNA expression of CYB561 was an independent prognostic factor for shorter OS and DSS of breast cancer patients, the data analyzed in this study was obtained from online databases. Prospective studies were needed to validate our findings and explore the clinical application of CYB561 in treatment of breast cancer. Second, there are no wet-lab experimental data in this study investigating the potential mechanisms of CYB561 in breast cancer. Further exploration is required to reveal the detailed mechanism of CYB561 in breast cancer as well as in other carcinomas.

Conclusion

In conclusion, the current study suggests that elevated CYB561 expression was significantly correlated with poor survival and immune infiltrations in breast cancer patients from multiple cohorts. CYB561 is a clinically promising biomarker that can be used to predict outcome and guide treatment. Our study provides new and promising insights for subsequent research to elucidate the molecular mechanisms of CYB561 in breast cancer.