COAD gene expression profiles and corresponding clinical information were downloaded from the TCGA database (n = 512) and GEO database (GSE39582, n = 576). Seventy iron metabolism-related genes are listed in Supplementary Table S1, which were derived from a previous study (Wei et al., 2020). The GSE39582 dataset was used to validate the established signature.

The TCGA dataset was used as the training dataset to build the iron metabolism-related genes risk score. Univariable Cox regression and Least absolute shrinkage and selection operator (LASSO) Cox analysis were performed to select iron metabolism-related genes associated with prognosis (p < 0.05). Then, multivariate Cox regression was used to further determine candidate iron metabolism-related genes associated with prognosis in COAD. Following this, the risk score for each patient was calculated as follows: R i s k   S c o r e = c o e f f i c i e n t   1 E ∗ x p 1 + c o e f f i c i e n t   2 E ∗ x p 2 + ....... + c o e f f i c i e n t   n E ∗ x p   n where Exp represented the expression level of the candidate iron metabolism-related gene. Using the median of the above risk score as the cutoff point, patients with COAD in the training cohort and validation cohort were divided into low-risk and high-risk subgroups.

Kaplan–Meier analysis was carried out to compare the prognostic difference between the two subgroups. Univariable and multivariate Cox regression were conducted to further evaluate the iron metabolism-related gene signature’s prognostic value in the training cohort. Then, the prognostic gene signature was validated externally in the validation dataset, followed by further Kaplan–Meier analysis.

Univariate and multivariate Cox regression analysis were carried out to determine the prognostic factors including the risk score established above, and age, sex, tumor stage, and tumor-node-metastasis (TNM) stage in COAD. The factors with a significant difference (p < 0.05) were selected to construct a nomogram model. A calibration curve was plotted to verify the accuracy of the nomogram.

Enriched pathways in different iron metabolism risk score datasets were assessed using Gene Set Enrichment Analysis (GSEA) software (GSEAv4.0.1, https://www.gsea-msigdb.org/). The Hallmark (v7.1) gene set collections were used as references. p < 0.05 and a false discovery rate of q < 0.25 were considered significant.

The abundance of 22 TICs in each tumor sample in the COAD cohort were estimated using the CIBERSORT module in the R package. CIBERSORT (Chen et al., 2018) was performed to determine the relative mRNA expression levels in the high- and low-risk subgroups to characterize the cell composition of tumor tissues.

Total RNA from colon cancer cells were extracted using the TRIzol reagent (Takara Biotechnology, Dalian, China). Next, cDNA was prepared using a Revert Aid First Strand cDNA Synthesis kit (Thermo Fisher Scientific, Waltham, MA, United States). The cDNA was then used as a template in a quantitative real-time polymerase chain reaction (qPCR) to determine the expression levels of the iron metabolism-related genes using Japan). The qPCR amplification reactions conditions were as follows: 95°C for 15 min; followed by 40 cycles of 95°C for 30 s, 55°C for 1 min, and 72°C for 30 s. The expression levels were normalized to those of ACTB (encoding β-actin). All primers were synthesized by Sangon Biotech (Shanghai, China) and are listed in Supplementary Table S2. All PCR reactions were performed in triplicate, and the relative expression levels of mRNA were quantified using the 2-∆∆Ct method.

This protocol was approved by the Ethics Committee of Sun Yat-Sen University Cancer Center (Guangzhou, China). Written informed consent was obtained from patients at their first visit. The surgically resected colon cancer tissues of eight patients at Sun Yat-sen University Cancer Center were included in this study. Tumors and corresponding nontumorous tissue were fixed in 4% paraformaldehyde, embedded in paraffin blocks, and processed into 4 µm-thick continuous sections. Immunohistochemical staining was performed to determine the distribution of SLC39A8 and SLC48A1. The antibodies used were: anti-SLC39A8 (1:200; 20459-1-AP; Proteintech) and anti-SLC48A1 (1:200; NBP1-91563; Novus Biologicals).

Colon cancer cells (DLD-1) were cultured in Dulbecco’s modified Eagle Medium (Gibco, Grand Island, NY, United States) supplemented with 10% fetal bovine serum (Invitrogen, Carlsbad, CA, United States) at 37°C with 5% CO₂. Small interfering RNAs (siRNAs) targeting SLC39A8 and SLC48A1 were synthesized by Sangon Biotech. The sequences of the siRNAs are listed in Supplementary Table S4. Transient transfection was performed using the Lipofectamine 2000 reagent (Invitrogen, Shanghai, China) according to the producer’s protocol.

DLD-1 Cells were seeded to a 96-well plate and transfected with negative control (NC)-siRNA, SLC39A8-siRNA, or SLC48A1-siRNA for 24, 48, or 72 h, after which the cells were further incubated with 20 μL of MTS reagent for another 2 h. Cell viability was detected as the optical density (OD) value at 490 nm.

To study the effects of SLC39A8 and SLC48A1 on cell proliferation, DLD-1 Cells transfected with NC-siRNA, SLC39A8-siRNA, or SLC48A1-siRNA were seeded into 6-well plates and incubated for 14 days. Cells were then stained using Crystal Violet Staining Solution, and the number of colonies was detected using light microscopy.

All statistical analyses were conducted using the R software. Kaplan–Meier analysis was used to compare the overall survival differences between the low- and high-risk subgroups. Given the possibility of multicollinearity, we used LASSO Cox regression analysis to identify the most valuable prognostic genes among all iron metabolism-related genes. An L1 penalty was set in the LASSO Cox model to shrink some regression coefficients to exactly zero and a 10-fold cross-validation with minimum criteria was performed to find the optimal λ value. Univariate and multivariate Cox proportional hazard regression analyses were carried out to evaluate the relationship between the risk score and overall survival. In all instances, a value of p < 0.05 was regarded as statistically significant.

As shown in Figure 1, a systematic study was carried out for the pivotal roles and the latent prognostic values of iron metabolism-related genes in COAD. The mRNA expression profiles in tissues from patients with COAD (n = 512) were downloaded from the TCGA database as the training dataset. The GSE39582 dataset (n = 576) was used as validation cohort. Of the 56,753 genes in the TCGA expression data, 70 iron metabolism-related genes were selected.

In the TCGA training dataset, single-factor Cox analysis was performed to analyze comprehensively the prognostic value of iron metabolism-related genes in COAD. We found that six genes [HAMP (encoding hepcidin antimicrobial peptide), SFXN3 (encoding sideroflexin 3), SLC22A17 (encoding solute carrier family 22 member 17), SLC39A14 (encoding solute carrier family 39 member 14), SLC39A8, and SLC48A1] were significantly related to the prognosis of patients with COAD (Supplementary Figure S1). These iron metabolism-related genes were subsequently subjected to LASSO Cox regression analysis to avoid colinear influences, and regression coefficients were calculated (Figure 2A). When the six genes were incorporated, the model achieved the best performance (Figure 2B). These genes and their related coefficients are shown in Supplementary Table S3. The six iron metabolism-related genes were further analyzed by multivariate Cox regression. As shown in Figure 2C; Table 1, SLC39A8 and SLC48A1 were identified as independent predictors of prognosis for patients with COAD.

The two hub iron metabolism-related genes identified in Construction of the Iron Metabolism-Related Gene Signature in the TCGA Cohort were used to construct a prognosis-related genetic risk score. The risk score of each patient with COAD was calculated as follows: R i s k   S c o r e =   − 0.052106 ×   E x p S L C 39 A 8 + 0.173944 ×   E x p S L C 48 A 1

We classified risk scores using the optimal cutoff points decided by the maximally selected log-rank statistics, in which patients with risk scores above the cutoff value were recognized as the high-risk group, and patients with risk scores below the cutoff value were recognized as the low-risk group. A survival analysis was performed to evaluate the predictive effect of this model. Figure 3A shows the distribution of the TCGA training cohort. Furthermore, dot plots were made to compare survival of patients in the high-risk and low risk-groups, which showed that survival in the high-risk group was worse than that in the low-risk group (Figure 3B). The heat maps in Figure 3C show the comparison of the expression levels of the two iron metabolism-related genes between the groups. The expression of SLC39A8 was higher in the low-risk patients, while SLC48A1 expression was higher in high-risk patients. Kaplan–Meier survival analyses found that patients with COAD in the high-risk group had worse prognosis compared with those in the low-risk group (Figure 3D; p < 0.0001).

To determine the accuracy of the two-gene prognostic model, we used the GSE39582 dataset as an external validation cohort. The distribution of the validation cohort is displayed in Figure 4A, Supplementary Figures S2A, S3A. The dot plot and heat map results were similar to those of the TCGA cohort (Figures 4B,C, Supplementary Figures S2B,C, S3B,C). Survival analysis indicated that the overall survival of the low-risk group was markedly longer than that of the high-risk group (Figure 4D; p = 0.034, Supplementary Figure S2D; p = 0.053 and Supplementary Figure S3D; p = 0.021).

We next analyzed other co-variates using univariate and multivariate Cox regression to determine the prognostic factors in COAD. The results showed that risk score, N stage, and M stage correlated with survival in the multivariate analysis. The hazard ratios of the risk score and above clinical factors are listed in Figure 5A. We then constructed a nomogram using the risk score, N stage, and M stage as variables (Figure 5B). A lower point was related to a better prognostic result on the nomogram. The C-index was 0.8092 [95% confidence interval (CI):0.7543–0.8642] for the nomogram, indicating that it had a good discriminating ability. Outcome was reported as 1-, 2- and 3-years overall survival. The associated calibration curves from the nomograms at 1, 2 and 3 years are displayed in Figures 5C,D, which showed a good performance of our nomogram in predicting 1-, 2- and 3- year survival. Therefore, this nomogram could be a useful model to predict the survival of patients with COAD.

To explore the molecular functions of the identified iron metabolism-related genes in this study, we conducted GSEA analysis to determine the gene expression profile. The results are shown in Figure 6. Multiple functional gene sets were enriched significantly in the low-risk group containing HALLMARK_ANDROGEN_RESPONSE, HALLMARK_KRAS_SIGNALING_UP, HALLMARK_PI3K_AKT_MTOR_SIGNALING, HALLMARK_PROTEIN_SECRETION, and HALLMARK_UV_RESPONSE _DN. These findings suggested that the two iron metabolism-related genes were potentially closely correlated with the status of the COAD microenvironment.

To further determine the relationship of the two iron metabolism-related genes with the immune microenvironment, we analyze the proportion of TIC subpopulations and constructed immune cell profiles in COAD using CIBERSORT. As shown in Figures 7A,B, a stacked bar plot and heat map were provided to describe the immune microenvironment in the high-risk and low-risk groups. Furthermore, the proportions of 22 immune cell proportions of COAD are shown in Figure 7C. The results showed that two TICs, regulatory T cells and eosinophils, were related to the two-gene signature risk score. The low- and high-risk score groups showed specific immune cell distributions.

To ascertain the expression levels of SLC39A8 and SLC48A1 in COAD tissues, eight colon cancer tissues and corresponding normal tissues were tested. Immunohistochemistry (Figure 8A) showed that the SLC39A8 level was downregulated in colon cancer tissues, and the SLC48A1 level was upregulated significantly in colon cancer tissues. Furthermore, the mRNA levels of SLC39A8 and SLC48A1 were successfully knocked down using siRNA in colon cancer cells in vitro (Figure 8B, Supplementary Figure S4). Next, we analyzed the potential function of SLC39A8 and SLC48A1 in colon cancer. Knocking down SLC39A8 promoted the proliferation of different colon cancer cell lines (DLD1, HCT116, and HT29), and silencing SLC48A1 suppressed the proliferation of colon cancer cells in vitro (Figures 8C,D, Supplementary Figure S4). These results indicated that SLC39A8 might function as a tumor suppressor gene and SLC48A1 might function as an oncogene in COAD.