Identification and verification of prognostic genes related to zinc homeostasis and zinc transport in breast cancer

Data extraction

The Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA) dataset, including survival information, age, and tumor node metastasis (TNM) staging, was downloaded from The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov/). Patients without survival information were excluded, resulting in tissue from 1,082 BC and 113 control specimens (Li et al., 2022). The GSE20685 dataset was obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). It comprises 327 BC tumor tissue specimens with complete survival information, sequenced on the GPL570 platform (Luo et al., 2022). A search was conducted in Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/) for the keywords “WP_ZINC_HOMEOSTASIS” (Human Gene Set, WP3529, 37 genes) and “GOBP_ZINC_ION_TRANSPORT” (Human Gene Set, GO:0006829, 28 genes) to retrieve ZHTGs. After removing duplicates, a total of 40 ZHTGs were retained (Table S1).

Differential expression analysis

The TCGA-BRCA dataset underwent analysis for differential gene expression employing the DESeq2 package (v 1.34.0) (Love, Huber & Anders, 2014) to evaluate the differentially expressed genes (DEGs) between BC and control specimens (p.adjust <0.05 and —log2FoldChange (FC)—>1).

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) is a bioinformatics approach based on gene expression correlation to construct modular networks. It clusters genes with similar expression patterns, analyzes the associations between modules and specific traits/phenotypes, and identifies key functional gene groups (Langfelder & Horvath, 2008). In this study, WGCNA was applied to analyze BC transcriptome data to identify co-expression modules associated with ZHTGs, thereby screening their related genes. First, univariate Cox regression analysis based on the expressions of 40 ZHTGs in the TCGA-BC dataset was performed using the survival package (v 3.5-3) (Lei et al., 2023) (Hazard Ratio (HR) ≠ 1 & p value < 0.05) to screen ZHTGs linked to BC prognosis for subsequent analysis. Following this, utilizing the ZHTGs, single-sample Gene Set Enrichment Analysis (ssGSEA) was employed to compute ZHTGs scores for BC specimens. These BC specimens were subsequently segregated into high and low score groups using the ideal threshold value of the scores. The survminer package (v 0.4.9) (Liu et al., 2021) was used to plot Kaplan–Meier (K–M) survival curves among high/low score cohorts (p value < 0.05). Next, the WGCNA package (v 1.70.3) (Langfelder & Horvath, 2008) was employed to construct a co-expression network using ZHTGs scores as traits, aiming to identify module genes most correlated with the traits. Particularly, all BC samples were grouped to detect and eliminate anomalies. A soft threshold (β) with a connectivity close to 0 and an R2 value great than 0.85 was selected. A scale-free network was built using selected soft threshold, and a hybrid dynamic tree cutting algorithm with a cutting tree parameter of 0.4 (the minimum gene number per module is 100, and the module merging parameter is 0.25) was used to identify the co-expression modules. Subsequently, the correlation coefficients between these modules and the traits were calculated, and the top two modules with the strongest positive and negative correlations were selected as key modules. Finally, by setting the gene significance (GS) to 0.3 and module membership (MM) to 0.3, key module genes were filtered and retained to select more important key module genes.

Enrichment analyses

The obtained DEGs were intersected with more important key module genes to obtain a set of candidate genes. To investigate the function in which the candidate genes were involved, Gene Ontology (GO, p.adjust < 0.05) and Kyoto Encyclopedia of Genes and Genomes (p value < 0.05, KEGG) enrichment analyses were performed using the clusterProfiler package (v 4.2.2) (Wu et al., 2021).

Development and verification of risk models

Within the TCGA-BRCA dataset, candidate genes were subjected to univariate Cox regression analysis via the survival package (Lei et al., 2023) to identify genes associated with BC prognosis, with HR ≠ 1 and p value < 0.01 as the screening criteria. Subsequently, a least absolute shrinkage and selection operator (LASSO) regression analysis was executed using the glmnet package (v 4.1-4) (Sasikumar et al., 2022) grounded on the findings of the previous step to further refine prognostic genes. The risk score for each individual was computed using the subsequent equation: ${\displaystyle \text{risk score}=\sum _{\mathrm{i}=1}^{\mathrm{n}}\text{coef}\left(\text{genei}\right)\text{*}\text{expr}\left(\text{genei}\right).}$

Individuals were subsequently classified into high/low risk cohort according to the median risk score. Subsequent Kaplan–Meier (K-M) survival analysis using the survminer package (Liu et al., 2021) to plot K-M curves (p-value < 0.05). The survival receiver operating characteristic (ROC) package (v 1.0.3) (Heagerty, Lumley & Pepe, 2000) was utilized to generate ROC curves for 1-, 3-, and 5-year periods to assess the prognostic performance of BC patients. The external dataset GSE20685 for BC patients was employed to validate the constructed risk model in this study.

Recognition of independent prognostic factors

To further investigate the prognostic implications of clinical pathological factors such as Age, Race, T/N/M stage were included along with the risk score in the prognostic model for univariate Cox regression analysis (HR ≠ 1, p value < 0.05) in BC specimens from the TCGA-BRCA dataset. Subsequently, a proportional hazards (PH) assumption test was conducted to select genes (p value > 0.05). Factors fulfilling the PH presumption were subsequently incorporated in multivariate Cox regression analysis, where factors with a p value < 0.05 were defined as independent prognostic factors.

Construction of nomogram

In order to further investigate the prognostic implications of independent prognostic factors, a nomogram forecasting the longevity likelihood of BC specimens was constructed using the rms package (v 6.5-1) (Sachs, 2017). The anticipatory effectiveness of this model was assessed via calibration curves and ROC.

Correlation analysis of clinical characteristics

In the TCGA-BRCA dataset, rank sum examinations were used to scrutinize the variations in risk scores among subgroups with diverse clinical characteristics and the survival disparities among two risk groups in the subgroups with diverse clinical characteristics.

Gene set enrichment analysis

In the TCGA-BRCA dataset, based on the two risk cohorts, differential analysis was carried out via DESeq2 to calculate log2 fold change (log2FC). The log2FC values were then sorted from largest to smallest, followed by gene set enrichment analysis (GSEA) were carried out via clusterProfiler package (Yu et al., 2012), using the h.all.v2023.1.Hs.symbols. gmt from MisgDB database as the background gene set (FDR < 0.05).

Gene set variation analysis

To delve deeper into the variations in KEGG pathways among two risk cohorts in TCGA-BRCA, gene set variation analysis (GSVA) was executed utilizing the GSVA package (v 1.42.0) (Hänzelmann, Castelo & Guinney, 2013) based on the background gene set “c2.cp.kegg.v2023.1.Hs.symbols.gmt” to obtain enrichment scores for different pathways. The limma package (v 3.54.0) (Love, Huber & Anders, 2014) was then utilized to compare the functional enrichment pathways among two risk cohorts, with a threshold of p value < 0.05 to select key pathways.

Immune-related analysis

In the TCGA-BRCA dataset, the enrichment scores of 28 immune-infiltrating cells for all samples among risk cohorts were computed using the ssGSEA technique. Subsequently, differential immune infiltrating cells were compared among different risk cohorts. Using the psych package (v 2.1.6) (Revelle, 2021), the Spearman correlation between differential immune infiltrating cells and prognosis genes was calculated (—cor— > 0.3 & p-value < 0.05). In the TCGA-BRCA dataset, we contrasted the levels of eight immune checkpoint molecules (CD274, LAG3, CTLA4, TIMP3, PDCD1, PDCD1LG2, TJAP1, LGALS9) among two risk cohorts. Additionally, the Tumor Immune Dysfunction and Exclusion (TIDE) online database (http://tide.dfci.harvard.edu/) was employed to acquire TIDE score, dysfunction score, and exclusion score for the specimens in the TCGA-BRCA dataset, and inter-cohort differences were compared. TIDE score is a computational method for evaluating tumor immune microenvironment functionality, which reflects the degree of tumor immune evasion. The dysfunction score is an indicator derived from the expression patterns of T cell exhaustion-related genes, assessing whether tumor-infiltrating T cells lose effector functions due to chronic antigen stimulation. The exclusion score evaluates whether a tumor forms an “immune-excluded” microenvironment by analyzing genes related to tumor stromal fibrosis, abnormal angiogenesis, and other factors that inhibit immune cell infiltration. To evaluate the ratio of immune therapy response in two cohorts, chi-square test was executed to examine the variations in immune reaction.

Cell culture and tissue samples collection

MCF-10A, T47D, BT474, MDA-MB-231, SUM-149, SUM-159 and MCF-7 were purchased from American Type Culture Collection. Due to the different molecular types of breast cancer and the existence of heterogeneity, we selected immortalized normal epithelial cells MCF-10A as the control cell line; According to the different expression statuses of ER, PR and HER-2, six breast cancer cell lines, namely T47D, BT474, MDA-MB-231, SUM-149, SUM-159 and MCF-7, were selected as the experimental groups to analyze the expression of related genes in the breast cancer cell lines. MCF-10 A cells were grown in DMEM/F12 medium supplemented with 5% horse serum and growth supplements (Zhong Qiao Xin Zhou Biotechnology, Shanghai, China). T47D, BT474, MDA-MB-231, SUM-149, SUM-159 and MCF-7 cells were cultured in DMEM (Gibco, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, Waltham, MA, USA). All cell lines were cultured in 5% CO2 at 37 °C in a humidified atmosphere. We obtained 10 pairs of BC tissue and adjacent normal tissue from patients without preoperative chemotherapy, endocrine therapy, or radiotherapy who had undergone tumor resection at The Xijing Hospital of the Fourth Military Medical University. The study was approved by the hospital’s ethics committee (KY20232266-C-1), and the content of the study received written informed consent from patients. All methods were performed in accordance with the relevant guidelines and regulations.

RNA isolation and reverse transcription quantitative PCR (RT-qPCR)

Total RNA was extracted from cells or tissues using SPARKeasy Improved Tissue/Cell RNA Kit (Sparkjade, Shandong Sparkjade Biotechnology Co., Ltd., Shandong, China) following the manufacturer’s instructions. cDNA was synthesized using SPARKscript II All-in-one RT SuperMix for qPCR (Sparkjade). The reverse transcriptional reaction program lasted 15 min at 50 °C and 5 s at 80 °C. We performed qPCR using the Quant Studio 7 Pro (Applied Biosystems). RT-qPCR was performed using 2 ×SYBR Green qPCR Mix (Sparkjade). RT-qPCR reaction program: preincubation at 94 °C for 3 min, 40 cycles of amplification with 10 s at 94 °C, 20 s at 55−60 °C, followed by an extension at 72 °C for 30 s. The internal controls were β-actin. Gene expression levels were quantitatively calculated by the 2-ΔΔCt method. Table S2 provides sequences of primers used in this research.

Immunohistochemistry

The Human Protein Atlas (HPA) (https://www.proteinatlas.org/) (Uhlen et al., 2017) contains IHC profiles of normal and tumor tissues.

Statistical analysis

In R software (v 4.1.0; R Core Team, 2021), the data was processed and analyzed. Differences among cohorts were assessed using the Wilcoxon rank-sum test or chi-square test, with a significance threshold of p value < 0.05 indicating statistical significance. Continuous variable data were analyzed using Student’s t-test or Wilcoxon test, and categorical data were analyzed using chi-square test. Survival differences were compared using Kaplan–Meier analysis and log-rank test. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

ZHTGs were associated with survival of BC

In the TCGA-BRCA dataset, there were 4,948 DEGs, consisting of 2,986 genes that were up-regulated and 1,962 genes that were down-regulated (Fig. 1A). Through univariate Cox regression analysis, three ZHTGs associated with BC prognosis were selected, with SLC30A5 identified as a high-risk gene (HR > 1), while SLC1A1 and TMEM163 as low-risk genes (HR < 1) (Fig. 1B). Subsequently, based on the optimal cutoff value of three genes scores, the BC specimens in the TCGA-BRCA dataset were split into two score cohorts. It was found that patients in the high ZHTGs score cohort had significantly higher survival rates than those in the low ZHTGs score cohort (p value = 0.01), suggesting that ZHTGs may impact the outcome of BC patients (Fig. 1C). Then, the ZHTGs score was used as a trait to further screen its related genes through WGCNA. Initially, no obvious outlier specimens were observed in the dataset, hence no specimens needed to be excluded (Fig. S1A). By setting the β value to 6, we achieved an R² approaching 0.85 and a connectivity close to 0 (Fig. 1D). Then, 16 co-expression modules were identified using the mixed dynamic tree-cutting algorithm, excluding the grey module (Fig. S1B). Correlation analysis revealed that the brown module (R = 0.41) and the black module (R =  − 0.36) had the strongest correlations with trait, with the brown module containing 1,496 genes and the black module containing 766 genes (Fig. 1E). Finally, by setting the GS and MM thresholds to 0.3, a total of 443 more important key module genes were obtained, with 313 genes retained in the brown module and 130 genes retained in the black module (Figs. S1C–S1D).

There were six prognosis genes were screened in BC

There were 134 candidate genes were acquired by intersecting 4,948 DEGs with 443 key module genes (Fig. 2A). Subsequently, these candidate genes were further explored for potential enrichment of functions through GO and KEGG analyses. In the GO analysis, these candidate genes were found to be primarily associated with cell cycle and chromosome-related functions, such as mitotic cell cycle phase transition and chromosomal region (Fig. 2B). In the KEGG term, the candidate genes were also involved in pathways related to cell cycle, prostate cancer, and BC (Fig. 2C). Then, a total of six prognosis-related genes (CLIC6, EIF4E3, TFF1, TPRG1, RSPH1, PCSK6) were identified through univariate Cox regression analysis (p-value < 0.01 & HR < 1) (Fig. 2D). Furthermore, these six genes were further confirmed as prognostic genes using the LASSO algorithm (Fig. 2E).

The risk model in TCGA-BRCA and GSE20685 datasets has strong predictive performance

Based on the median of the risk score (Risk score = CLIC6 * (−0.0591) + EIF4E3 * (−0.1226) + TFF1 * (−0.0145) + TPRG1 * (−0.017) + RSPH1 * (−0.0272) + PCSK6 * (−0.0796)), BC specimen of TCGA-BRCA dataset were stratified into two risk cohorts. It was noticed that as the risk scores rose in the BC specimens, the number of deaths notably increased (Fig. 3A). Subsequently, from the K-M curves, it was evident that the survival rate of individuals in the high-risk cohort was lesser in comparison to those in the low-risk cohort (Fig. 3B). Furthermore, ROC curves were plotted at 1, 3, and 5 years as survival time nodes, demonstrating that the AUC values at all three time points were higher than 0.6, suggesting good forecasting ability of the model (Fig. 3C). Finally, the universality of the risk model was further validated in the GSE20685 dataset. The results indicated that similar to TCGA-BRCA, this risk model of this dataset could effectively predict patient survival (Figs. 3D–3F). Additionally, the GSE20685 dataset showed more distinct separation of survival groups than the TCGA-BRCA dataset, which further supported the validity of the model.

There were five independent prognostic factors

Within the univariate Cox regression analysis, the p-values of risk score, age, T/N/M stage were all less than 0.05 and passed the PH assumption (Fig. 4A). Subsequently, riskScore, age, T/N/M stage were screened as independent prognostic factors via multivariate Cox regression analysis (Fig. 4B). To extend the analysis of survival prediction in BC patients by these factors, a nomogram was developed (Fig. 4C). Calibration curve results indicated that the slope of the curve closely approximated the diagonal line, suggesting high prediction accuracy (Fig. 4D). Moreover, with AUC values consistently exceeding 0.7, this suggested that the model possessed robust predictive capability (Fig. 4E).

Clinical feature was associated with risk score and survival rate

By comparing the differences in risk scores between different clinical feature sub-cohorts to explore their relationships, significant differences were observed in the risk scores between the race (white and others) and T stage (T1 and T2, T1 and T3) sub-cohorts (Fig. 5A). Subsequently, the K-M survival differences between the high- and low-risk groups under different clinical characteristics were compared. The results revealed significant survival differences between the high-risk and low-risk groups across various clinicopathological characteristics, including age ≤60, T stage 1–2, N stage 1–3, and M stage 0. Notably, the survival rate of the high-risk group was consistently lower across all these characteristics (Fig. 5B).

The pathways with differences between the high- and low-risk groups

To explore the significantly enriched pathways between the high- and low-risk groups, GSEA and GSVA analyses were performed. The results revealed that DEGs among two risk cohorts were mainly involved in E2F TARGETS, G2M CHECKPOINT, and KRAS SIGNALING DN pathways (Fig. 6A). Subsequently, to further explore the differences in KEGG pathways between the high- and low-risk groups, GSVA was conducted. The results showed that there were 29 differential KEGG pathways between the two cohorts, including JAK-STAT signaling pathway, β-alanine metabolism, ribosome, drug metabolism-cytochrome P450, etc. (Figs. 6B–6C).

Patients in different risk cohorts had different effects on the immune response

For additional investigation into the variances in the immune microenvironment among two risk cohorts, a series of immune-linked analyses were performed. A heatmap was employed to exhibit the spread of 28 immune infiltrating cell enrichment scores among two cohorts in the TCGA-BRCA dataset (Fig. 7A). Following this, differential analysis uncovered substantial differences in the enrichment scores of 24 immune cells among the two risk cohorts. Specifically, central memory CD8 T cell, eosinophil, mature dendritic cell, mast cell, memory B cell, natural killer cells (NK cells), neutrophil, plasmacytoid dendritic cell had higher scores in the low-risk cohort, while other cells were opposite (Fig. 7B). Further correlation analysis revealed significant correlations between six prognostic genes and various differentially expressed immune infiltrating cells. Among them, the strongest negative correlation was observed among RSPH1 and activated CD4 T cell (cor = −0.46) (Fig. 7C, Tables S3–S4). Additionally, we found obvious variances in the expression of seven immune checkpoints (LAG3, CTLA4, TIMP3, PDCD1, PDCD1LG2, TJAP1, LGALS9) between high and low-risk cohorts, with all six immune checkpoints, except TIMP3, displaying elevated expression in the high-risk cohort (Fig. 7D). Furthermore, significant differences were observed in TIDE score and Dysfunction score among two cohorts, with higher scores in the low-risk cohort (Fig. 7E). Furthermore, notable variations were observed in immune treatment responses among two cohorts (Fig. 7F). The results indicated that there might be implications for the effectiveness of immunotherapy in different risk cohorts.

External and experimental validation of the six model related genes

Compared with the paired adjacent tissues, the expression levels of CLIC6, TFF1, TPRG1, RSPH1 and PCSK6 were up-regulated. The expression levels of EIF4E3 were down-regulated in BC specimens (Figs. 8A–8F). In addition, we also investigated the expression of prognostic genes in BC cell lines (including T47D, BT474, MDA-MB-231, SUM-149, SUM-159, and MCF-7), using MCF-10A as the control cell line. Compared with the MCF10A cell line, the gene expression level in most cell lines was overall consistent with the tissue verification results (Figs. 9A–9F). Finally, in order to verify the expression of the six genes involved in building the prognosis model, we downloaded immunohistochemical staining images from the HPA database. CLIC6, TFF1, TPRG1, RSPH1 and PCSK6 were expressed at notably different levels between BC and normal breast tissues. There was no notable difference in the expression of EIF4E3 (Figs. 10A–10F). This may be due to the low consistency between antibody staining and RNA expression data. These results supported our hypothesis and provided evidence for the rationality of selecting these six genes to build the prognosis model.