Human ccRCC cell line Caki-1 (catalog #HTB-46) was obtained from ATCC (Manassas, VA). Cells were kept in McCoy’s 5A Modified Medium (ATCC Catalog #30-2007) supplied with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin at 37°C in a 5% CO₂ setting. Tyrosine kinase inhibitor Sorafenib was purchased from Cayman Chemicals (Ann Arbor, MI) and dissolved in DMSO at 10 mM stock concentration. Bone morphogenetic proteins (BMP) and interleukin cytokines were synthesized by PeproTech (Cranbury, NJ) and dissolved in a 1% bovine serum albumin (BSA) carrier solution at a concentration of 10 mg/ml. The enzyme-linked immunoassay (ELISA) kit for the hepcidin assay (Catalog #EKF58061) was purchased from BioMatik (Wilmington, DE).

Caki-1 cells were seeded in a 6-well plate at a density of 4 × 10⁴ cells overnight and then treated with Sorafenib and cytokines for 24 h, as indicated in the figure. Cell culture media (conditioned media) were collected at the end of experiments and spanned for 3 min to remove cell debris. The supernatants were subjected to the ELISA assay for hepcidin analysis. The ELISA assay was performed using the assay protocol provided by the supplier.

Gene expression profiles were assessed at the mRNA levels in malignant and normal kidney tissues using the RNA sequencing (RNA-seq) dataset from the TCGA project (Cancer Genome Atlas Research, 2008) and the cDNA microarray dataset (Yusenko et al., 2009). For the RNA sequencing datasets, fragments per kilobase per million (FPKM) data were downloaded from the TCGA online portal and converted to the log₂ value of Transcript Per Million reads (TPM) before statistical analysis. The RNA-seq data from 72 cases-matched normal and malignant kidney tissues were included for the pairwise comparisons. The cDNA microarray dataset was used to generate graphic images of grouped comparisons on the Oncomine™ platform (Rhodes et al., 2004). A ROC curve analysis compared HAMP expression levels in malignant and normal kidney tissues derived from ccRCC patients.

Patient survival outcomes were assessed using the TCGA-KIRC dataset with the Kaplan-Meier curve approach (Nagy et al., 2021). A minimum p-value approach was utilized to determine the cutoff value of gene expression levels for splitting patients into high or low subgroups (Menyhart et al., 2018). Patient survival analysis was also conducted based on gender, ethnic race, tumor mutation burden, and enriched or decreased immune infiltration types in the tumor microenvironment.

The correlation of gene expression was analyzed using the Spearman coefficient test. The correlation between HAMP expression and DNA methylation levels was analyzed using the TCGA-KIRC RNA-seq dataset (Cerami et al., 2012). The highly correlated genes (Spearman rho ≥ +/− 0.3) were subjected to gene enrichment analysis and visualization on the KOBAS platform (Wu et al., 2006). The significantly correlated genes were visualized using the R-package (version 3.6.3).

All quantitative data were presented as the MEAN with the SEM (standard error of the mean). Differences among multiple groups were analyzed using the Wilcoxon rank-sum test. Case-matched pair comparison was conducted using the Wilcoxon signed-rank test. The odds ratio was calculated using the binomial logistic regression analysis for HAMP expression on tumor progression. The differences in patient survival outcomes between the two subgroups were analyzed using the Kaplan-Meier survival plot. The hazard ratio (H.R.) was calculated using the Log-rank test (Menyhart et al., 2018). Each ROC curve was tested for effectiveness with the Youden index and compared relative potential with DeLong’s test.

We reasoned that a robust biomarker of prognostic factor and therapeutic target should be biologically relevant and tightly correlated with disease aggressiveness and patient survival outcome. With this notion in mind, we examined seventy-one iron-modulating genes as reported (Miller et al., 2011; Brown et al., 2020; Schnetz et al., 2020; Torti and Torti, 2020). Gene expression levels were compared between different groups of ccRCC patients stratified by tumor stage (early vs. late), lymph node invasion, distal metastasis, cancer grade (low vs. high), progression-free interval (PFI), overall survival (OS), and disease-specific survival (DSS). We looked for the genes that exhibited a strong and concordant alteration among all the clinicopathological parameters (Supplementary Figure S1). Among these 71 genes, five genes (HAMP, HMBS, ISCA2, STEAP2, and STEAP3) presented concordant alterations in all comparisons (Supplementary Table S1). Of these five genes, the expression levels of HAMP, HMBS, and STEAP3 genes were significantly increased along with disease progression and correlated with worse survival outcomes. In contrast, STEAP2 and ISCA2 gene expression was significantly decreased along with disease progression and correlated with favorite survival outcomes.

We examined gene expression alterations in a pairwise comparison of case-matched malignant and normal tissues derived from 72 patients in the TCGA-KIRC dataset. As shown in Figure 1A, the expression levels of four genes (HAMP, ISCA2, STEAP2, and STEAP3 but not HMBS) were significantly altered in malignant tissues compared to normal kidney tissues, of which HAMP gene expression displayed a consistent upregulation in all but one case. We then conducted a ROC analysis to determine their relative potentials in distinguishing malignant and normal kidney tissues. As shown in Figure 1B, HAMP expression was highly superior (AUC = 0.911, CI95: 0.879-0.944) over other four genes (AUC 0.594–0.791). These data indicated that HAMP upregulation is a potent diagnostic factor in ccRCC patients.

We then analyzed HAMP expression in pRCC and ChRCC tissues with TCGA RNA-seq datasets. Our analysis revealed a significant increase of HAMP expression in pRCC and ChRCC tissues (Figures 1C,D). We then confirmed HAMP upregulation in ccRCC, pRCC, and ChRCC tissues (Figures 1E–G) with a secondary dataset derived from a cDNA microarray assay (Yusenko et al., 2009). In addition, HAMP upregulation was also evidenced in other human cancer types derived from the breast, colon, lung, head-neck, stomach, and esophagus (Supplementary Figure S2). However, HAMP was remarkably downregulated in human liver cancers, as described in our recent publication (Wang et al., 2021).

We also used a qualitative approach to analyze the associations of HAMP upregulation with clinicopathological parameters. Patients were divided into HAMP^high or HAMP^low groups based on the median level of HAMP expression (Table 1). In ccRCC patients, there were significantly more HAMP^high cases than HAMP^low cases in all clinicopathological parameters. Consistent with the functional role of HAMP in lowing circulating iron levels (Wojciechowska et al., 2021), HAMP^high group showed significantly more cases with lower hemoglobin levels than HAMP^low group. However, in pRCC and ChRCC patients, HAMP expression was not associated with these clinicopathological parameters (Supplementary Tables S2, S3). These results suggested that HAMP upregulation was tightly associated with clinicopathological parameters, in line with a recent report showing HAMP upregulation in chemical-induced rat kidney cancers (Matsumoto et al., 2017).

Epigenetic modifications, such as promoter DNA methylation, are common mechanisms for transcriptional gene regulation (Angulo et al., 2021). We explored if HAMP upregulation was due to promoter demethylation. Our exploration found a significant correlation between HAMP mRNA expression and gene hypomethylation in two datasets derived from two different DNA methylation sequencing platforms, HM27 (Figure 2A) and HM450 (Figure 2B). Further analysis revealed a strong correlation between HAMP gene expression and DNA hypomethylation around the Exon-1/Intron-1 region (Figure 2C). Quantitatively, HAMP promoter methylation levels were significantly lower in ccRCC tissues than in normal kidney tissues (Figure 2D), which further decreased in tumors with lymph node invasion (Figure 2E). On the other hand, HAMP gene amplification was only found in one out of 751 (0.13%) ccRCC cases (Supplementary Figure S3). These data indicate that epigenetic modification (DNA hypomethylation) at the Exon-1/Intron-1 region, but not a structural/genetic abnormality, might be a potential mechanism of HAMP upregulation in ccRCC tumors.

It has been reported that HAMP expression was stimulated by BMP proteins (BMP2, 4, 6, 7), inflammatory cytokine interleukin-6 (IL6), and fibroblast growth factor-6 (FGF6) (Vela and Vela-Gaxha, 2018). We analyzed their correlation with HAMP expression in ccRCC tumors. As shown in Figure 3A and Table 2, only IL6, but not BMP proteins and FGF6, were significantly correlated with HAMP expression in ccRCC tissues. To confirm the IL6 effect on HAMP expression, we treated ccRCC Caki-1 cells with IL6 and measured hepcidin production in the cell culture media. The kinase inhibitor Sorafenib was used as a positive control since it was reported to enhance HAMP expression (Mleczko-Sanecka et al., 2014). As shown in Figure 3B, Sorafenib treatment significantly increased hepcidin levels, as expected, whereas IL6 only had a weak but not significant enhancing effect on hepcidin levels, which was not in line with the correlation analysis data.

We then conducted a whole-transcriptome correlation analysis with HAMP expression to explore other interleukin molecules that were associated with HAMP upregulation. Our analysis found that IL27, IL23A, IL10, and IL34 were strongly correlated with HAMP expression (Spearman coefficient r > 0.35, Table 2). We then examined their causative effects on hepcidin expression in Caki-1 cells. Surprisingly, only IL34 significantly stimulated hepcidin expression (Figure 3B). Further analysis revealed that IL34 but not IL6 gene expression was remarkably increased in ccRCC tissues (Figure 3C). These data suggested that IL34 is a potential driving factor for HAMP upregulation in ccRCC tissues.

We then examined the association of HAMP upregulation with patient survival outcomes. As shown in Figure 4A, ccRCC patients with higher levels of HAMP expression had significantly worse survival outcomes (overall survival, disease-specific survival, and progression-free interval). This negative impact was statistically significant in the White race group with a statistical significance (Figure 4B) but not in the Black/Asian group of patients (Figure 4C). In patients with a higher tumor mutation burden (TMB), HAMP upregulation showed a substantial impact on overall survival than those with a lower TMB (HR 5.29 vs. 3.75, Figure 4D). However, the survival impacts were similar between male and female ccRCC patients (Figure 4E). These results demonstrated HAMP upregulation was a potent prognostic factor of unfavorable survival outcomes in the White race group of ccRCC patients with higher TMB levels.

We compared patient survival outcomes based on immune infiltration enrichment to evaluate the clinical significance of immune infiltration expression concerning HAMP expression. Our analysis revealed that HAMP upregulation-associated worse survival outcomes were diminished in ccRCC patients who showed an intratumoral enrichment of natural killer (NK) cells, T Helper cells (T_H1/2), and mesenchymal stem cells (Figure 5). Conversely, other immune infiltration enrichment, including eosinophils, macrophages, Treg cells, CD4⁺ memory cells, cytotoxic CD8⁺ T-cells, and B-cells, did not affect the HAMP upregulation-associated worse survival outcomes (Supplementary Figure S4). These results suggest that HAMP upregulation-related negative impact on patient survival only occurred when the tumor immune microenvironment was weakened due to reduced infiltration of anti-cancer immune cells.

We next conducted a gene-set enrichment analysis (GSEA) to determine the most relevant pathways correlated with HAMP expression. Based on the whole-transcriptome correlation analysis on the TCGA-KIRC dataset, 1725 genes were positively correlated with HAMP expression, while 1188 genes were negatively correlated with HAMP expression based on a Spearman r cutoff value ± 0.3. These genes were then subjected to the GSEA analysis. Our results showed that the most significantly enriched clusters/pathways were immune regulatory and inflammatory cytokine/chemokines, as shown in Table 3 and Supplementary Figure S5.

We then focused our analysis on immune-related genes (Huang et al., 2019; Simonaggio et al., 2021). A strong correlation (Spearman r > 0.35) was identified between HAMP expression and four major inhibitory checkpoint genes, PDCD1, LAG3, TIGIT, and CTLA4 (Figure 6A). In parallel with HAMP upregulation, these four checkpoint genes were also significantly upregulated in ccRCC tissues compared to normal kidney tissues, either in case-match pairwise (Figure 6B) or group cohort comparison (Figure 6C). Their upregulations were also associated with a significantly worse disease-specific survival and progression-free interval (Figures 6D–G), which were in line with a recent report (Zhang et al., 2019). These results indicate that HAMP and immune checkpoint genes (PDCD1, LAG3, TIGIT, and CTLA4) were coordinately upregulated in ccRCC tissues with a strong clinical impact on patient survival outcomes.