Hepcidin expression levels in various types and subtypes of cancers in TCGA project and healthy controls were obtained from the database of Tumor IMmune Estimation Resource version 2 (TIMER2, timer.cistrome.org) (Li et al., 2017; Li et al., 2020). The types of common male genitourinary system tumors that we focused on in this study were as follows: bladder urothelial carcinoma (BLCA), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), prostate adenocarcinoma (PRAD), and testicular germ cell tumors (TGCT). For expression analysis, box plots of Gene Expression Profiling Interactive Analysis version 2 (GEPIA2, gepia2. cancer-pku.cn) were applied for types of cancer with limited information of healthy control (Tang et al., 2019). Also, hepcidin expression in different stages of cancer was summarized as in violin plots using GEPIA2. The transformed transcriptional levels were calculated as log2 of transcripts per million (TPM) +1 and were used in box and violin plots.

The protein-level analysis of the data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) (Edwards et al., 2015) was conducted using the UALCAN portal (ualcan.path.uab.edu/analysis-prot.html) (Chandrashekar et al., 2017). We explored the expression level of hepcidin between primary tumor and normal tissues, respectively. The correlation between protein expression of hepcidin and patients’ age, gender, disease stage, tumor histology, and body weight was further analyzed. The available datasets were the data of KIRC.

Both overall survival (OS) and disease-free survival (DFS) of hepcidin in various TCGA cancers were obtained using the survival map section of GEPIA2. The threshold for different expression levels was set as high expression (50%) and low expression (50%). A hypothesis test was carried out using the log-rank test. Pooled analysis (univariate Cox regression) of OS, disease-free interval (DFI), progression-free interval (PFI), and disease-specific survival (DSS) was conducted across male genitourinary system tumors. The forest plot was drawn without hazard ratio (HR) merging. Then, the receiver operator characteristic curve (ROC) was performed to determine whether the expression of hepcidin could exactly predict the 1-, 3-, and 5-year survival (OS) of patients in selected types of tumors from TCGA database. The best cutoff value for the expression of hepcidin was calculated. Patients of each type of tumor studied were divided into the high-expression group and low-expression group by using the best cutoff value. Survival analysis was used to validate the OS between groups that were stratified by cutoff value.

The genetic alteration features of hepcidin were analyzed using TCGA pan-cancer studies in cBioPortal (cbioportal.org) (Gao et al., 2013). Alteration frequencies, mutation types, and copy number alterations in the HAMP gene in selected types of TCGA cancers were analyzed by the cancer types summary. Mutation sites and 3D structures of hepcidin were also analyzed and shown in schematic illustration. The differences in hepcidin genetic alteration-associated OS, DFS, and progression-free survival (PFS) were analyzed. Kaplan–Meier (KM) plots with p-values obtained from log-rank tests were produced. Association between different DNA methylation probes and expression of hepcidin was carried out by MEXPRESS (mexpress.be) (Koch et al., 2015; Koch et al., 2019). The beta value of each sample was obtained. The p-value (adjusted using the Benjamini–Hochberg procedure) and Pearson correlation coefficient (R) value were obtained. MethSurv (Modhukur et al., 2018) (https://biit.cs.ut.ee/methsurv/) and SurvivalMeth (Zhang et al., 2021) (http://bio-bigdata.hrbmu.edu.cn/survivalmeth/) were used to measure the prognostic significance of DNA methylation of hepcidin in patients with male genitourinary system tumors.

The associations between hepcidin levels and infiltrations of CD8⁺ T cells and cancer-associated fibroblasts in selected types of TCGA cancers were analyzed by the Immune-Gene section of TIMER2. For estimating the immune infiltrations, algorithms including TIMER, CIBERSORT, CIBERSORT-abs, quanTIseq, XCELL, MCPCOUNTER, and EPIC were used. Both the p value and partial correlation (r) value were generated using purity-adjusted Spearman’s rank correlation test. All results were displayed in heatmaps and scatter plots.

The potential associations between the hepcidin level and tumor mutational burden (TMB), microsatellite instability (MSI), checkpoint expressions, immune cell score (MCPcounter), and immune response pathways in malignancy of the male genitourinary system were derived from Sangerbox (sangerbox.com/Tool) based on TCGA database. Spearman’s rank correlation tests were carried out, and the p-value and partial correlation (r) value were obtained. Meanwhile, the assessments of OS and relapse-free survival (RFS) were carried out using the interactive operation interface of the KM plotter (kmplot.com) using different data extracted from TCGA database to examine the correlation between patients’ prognoses and immune cell infiltration in tumor tissue. Multiple cellular contents were included in prognostic analysis such as enrichment of B cells, CD4⁺ memory T cells, CD8⁺ T cells, macrophages, mesenchymal stem cells, natural killer T cells, regulatory T cells, type 1 T-helper cells, type 2 T-helper cells. Median was used in the analysis of associations in BLCA, KIRC, KIRP, and TGCT. The HR, 95% confidence intervals, and log-rank p-value were computed, and the Kaplan–Meier survival plots were generated.

The STRING website (string-db.org) was used for the hepcidin-related gene enrichment analysis. The following parameters were used for the analysis: organism, “Homo sapiens”; minimum required interaction score, “low confidence (0.150)”; meaning of network edges, “evidence”; max number of interactors to show, “no more than 50 interactors in first shell”; and active interaction sources, “experiments”. Proteins that were able to bind or interact with hepcidin were generated according to published data of experimentally confirmed data. Meanwhile, the top 100 genes associated with hepcidin were obtained by the similar gene detection section of GEPIA2 according to the data of selected cancer and healthy controls. The correlation analysis was also performed for the Pearson correlation analysis of hepcidin and potential hepcidin-related genes. The log2 of TPM was used for the analysis of dot plots. The p value and r value were calculated. Also, the heatmap of related genes was generated by the Gene-Corr section of TIMER2, with r and p values calculated by purity-adjusted Spearman’s rank correlation tests. Intersection analyses were conducted for the comparison of hepcidin-related genes using Jvenn (Bardou et al., 2014). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was conducted using combined results obtained from previous analyses. The list of potential hepcidin-related genes was analyzed by the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Dennis et al., 2003). The enriched signaling pathways were listed in tidyr (cran.rproject org/web/packages/tidyr) and ggplot2 (cran.r-project.org/web/packages/ggplot2) of the R package. The cluster profiler R package (http://www.bioconductor. org/packages/release/bioc/html/clusterProfiler) was conducted for gene ontology (GO) enrichment analysis. Results of biological processes, cellular components, and molecular functions were plotted as cnetplots. R language software [R-4.0.4, 64-bit] (www.r-project.org) was used. Two-tailed analysis (p < 0.05) was set as statistical significance.

The studies involving human participants were reviewed and approved by the Ethics Committee of The Third Affiliated Hospital of Soochow University. The formalin-fixed paraffin-embedded (FFPE) KIRC tissues (N = 10) were obtained from the tissue bank of The Third Affiliated Hospital of Soochow University. The patients had provided their written informed consent when their tissues were stored in the tissue bank. Thus, it is not necessary to re-obtain written informed consent for this study. The FFPE KIRC tissues were used to verify the correlation between the expression of hepcidin and checkpoints in KIRC tissue. Immunofluorescent (IF) staining was performed in FFPE KIRC tissues. Common associated immune checkpoint molecules in male genitourinary system tumors such as CD28, CD48, LGALS9, CTLA4, and PDCD1 were selected for verification (for each aforementioned checkpoint molecular, one FFPE KIRC tissue for testing and another FFPE KIRC tissue for validating). Briefly, the sections were incubated with hepcidin (ab190775, Abcam, MA, United States, 1/100 dilution), CD28 (ab185759, Abcam, MA, United States, 2 μg/ml), CD48 (ab185759, Abcam, MA, United States, 2 μg/ml), LGALS9 (ab153673, Abcam, MA, United States, 5 μg/ml), CTLA4 (ab19792, Abcam, MA, United States, 1/200 dilution), and PDCD1 (ab52587, Abcam, MA, United States, 1/50 dilution) at 4°C overnight and then incubated with Alexa^®488-conjugated goat anti-rabbit secondary antibody or Alexa^®576-conjugated goat anti-mouse secondary antibody (Thermo Fisher, CA, United States). The nuclear stain Hoechst 34580 (5 μg/ml; Molecular Probes, Thermo Fisher, CA, United States) was added prior to final washes after the incubation of the secondary antibody. Finally, the sections were dehydrated, cleared, and mounted using a Zeiss confocal microscope. The resulting area and cell measurements were quantified using ImageJ software analysis.

Hepcidin mRNA expression in common male genitourinary system tumors from TCGA—including BLCA, KICH, KIRC, KIRP, PRAD, and TGCT—was analyzed using TIMER2 and GEPIA2 (Figure 1A, B). Compared to the corresponding healthy controls, elevated levels of hepcidin were found in kidney and testicular tumors: KICH, KIRC, KIRP, and TGCT (p < 0.001). For KIRC samples, this trend appears to reverse the trend seen in the protein levels.

Next, based on searches of the CPTAC database, we then compared the differences in total hepcidin protein levels between normal tissues and primary tumor tissues. Only the proteome data from patients with KIRC were analyzed. Interestingly, normal tissue expressed more of the protein than primary tumor tissues. This difference was also confirmed in subgroup analyses, controlling for patients’ age, gender, disease stage, tumor histology, and body weight. For patients with KIRC, there was no significant difference in hepcidin protein expression among different age-groups (21–40; 41–60; 61–80; 81–100 years). Similar results were confirmed in gender groups (male vs. female) and tumor grade groups (grade 1, 2, 3, and 4). It was found that patients with stage 1 cancer had higher hepcidin levels than patients with stage 3 cancer. Patients with normal weight expressed more hepcidin protein than those who were overweight (Figure 1C). The “Pathological Stage Plot” module of GEPIA2 was used to observe the correlation between hepcidin expression and the pathological stages of cancers. Hepcidin was found to be significantly correlated with the stages of BLCA and KIRC among the selected tumors. The mRNA expression of hepcidin is higher in patients with stages III and IV BLCA than those with stage II BLCA, while in KIRC, the mRNA expression of hepcidin increased with the advance of tumor stage. These differences were statistically significant (Figure 1D, p < 0.01).

Pooled analyses between hepcidin expression and patient OS, DFI, DSS, and PFI were conducted across the male genitourinary system tumors for which data were available in TCGA. Expression of hepcidin was an independent prognostic factor for OS, DSS, and PFI in KIRC (HR = 1.52, 95% confidence interval (95% CI): 1.31–1.76, p < 0.001; HR = 1.69, 95% CI: 1.42–2.01, p < 0.001; and HR = 1.44, 95% CI: 1.23–1.68, p < 0.001) but not for any clinical outcomes in any other selected tumors (Figure 2A; Table 1 ).

ROC analysis was then performed to determine the prediction efficiency of hepcidin expression for the 1-, 3-, and 5-year prognosis for OS, DFI, PFI, and DSS of patients from TCGA database. The best cutoff value for hepcidin expression was then calculated for each metric. Patients were divided into the high-expression and low-expression groups according to the best cutoff values. Survival analysis was then conducted to validate the prognosis between groups. We defined good prediction efficiency of hepcidin as the area under the curve (AUC) > 0.6, and a p-value of survival analysis <0.05. Hepcidin was found to have good prediction efficiency for OS, PFI, and DSS in KIRC (cutoff value 0.57, 1-year AUC 0.68, 3-year AUC 0.63, 5-year AUC 0.62; cutoff value 0.57, 1-year AUC 0.62, 3-year AUC 0.61, 5-year AUC 0.60; and cutoff value 0.52, 1-year AUC 0.67, 3-year AUC 0.66, 5-year AUC 0.67; respectively, all p < 0.001) and of DFI in BLCA (cutoff value 0.69, p = 0.0053, 3-year AUC 0.64, 5-year AUC 0.61) (Figure 2B). Hepcidin thus shows potential as a prognostic marker of KIRC.

The genetic alterations of the hepcidin gene in common male genitourinary system tumors (BLCA, KICH, KIRC, KIRP, PRCA, and TGCT) were also explored. As a result of mutation frequency summarizing, we found alterations in these tumors are rare, with a general occurrence rate of 1%. Amplification was the most common genetic alteration found in the studied cancer patients, while patients with BLCA had the highest proportion of altered hepcidin genes at about 2.5%. Detailed information on all the genetic alterations found is displayed in Figure 3A. S4_17del, W8G, and E36* mutations were the most frequent site mutations, which were found in one case each of BLCA, KIRC, and KIRP. These mutations were found to induce in-frame, missense, and nonsense mutations, respectively (Figure 3B).

The possible correlation between genetic alterations of the hepcidin gene and the prognosis of patients was then carried out. As shown in Figure 3C, the results suggest that the adverse prognosis of DFS in selected tumors was linked to the patients with hepcidin alterations. However, for each type of tumor, an adverse prognosis of DFS was only observed in patients with BLCA (Figure 3D) as the hepcidin alteration proportion of patients with other selected tumors is low. The association of the DNA methylation levels in the hepcidin gene and gene expression in selected tumors from TCGA datasets was analyzed using the visualizer MEXPRESS. The DNA methylation level was found significantly correlated with the hepcidin transcription level at different probes. Data from this analysis are summarized in Table 2 and Supplementary Figure S2. The heatmaps (MethSurv) of DNA methylation of hepcidin in male genitourinary system tumors (BLCA, KICH, KIRC, and KIRP) are shown in Figure 4A. Sites cg02131995 and cg27273033 were found highly methylated in all selected male genitourinary system tumors, and methylations of cg02131995 were associated with prognosis in BLCA and KIRC. Results are summarized in Table 3. Figures 4B–E show different expression patterns of hepcidin between low- and high-risk groups in BLCA (Figure 4B), KIRC (Figure 4C), KIRP (Figure 4D), and PRAD (Figure 4E) that were calculated by SurvivalMeth. The heatmap revealed the CpG methylation level of hepcidin in BLCA (Figure 4B), KIRC (Figure 4C), KIRP (Figure 4D), and PRAD (Figure 4E) which are based on TCGA database. Moreover, statistically prognostic association was only found between the high- and low-risk groups in KIRC patients. The results in the same type of cancer were in keeping with those from MethSurv in general.

As one of the main features of the tumor microenvironment, immune cell infiltration is tightly related to oncogenesis. Tumor-associated fibroblasts and CD8⁺ T cells located in the stroma of the tumor microenvironment were found to contribute to the regulation of the functions of different tumor-infiltrated immune cells. In this study, the possible correlation between immune cell infiltration and the hepcidin transcription level was analyzed using data from TCGA. In BLCA, KIRC, KIRP, and PRAD, hepcidin was found to be positively correlated with infiltrations of CD8⁺ T cells in most of the tested algorithms. However, in TGCT, most algorithms did not find any correlation between CD8⁺ T-cell infiltration and the level of hepcidin (Figure 5A). Moreover, we found that the hepcidin level was positively correlated with the estimated immune infiltrating value of tumor-associated fibroblasts in BLCA, KIRP, PRAD, and TGCT, while an inverse correlation was found in KICH and KIRC. Every representative correlation plot of the aforementioned cancers as estimated by a single algorithm is shown in Figures 5A, B.

However, there were no correlations found between hepcidin expression and TMB or MSI (Figure 6A). Furthermore, correlation analyses between immune checkpoints and hepcidin revealed positive correlations with most of the checkpoints studied in BLCA, KIRP, KIRC, and KICH (Figure 6B). About half of all checkpoints were found to be positively correlated with hepcidin in BLCA and TGCT. Among all the tested checkpoints, CD28, CD48, CTLA4, LAIR1, LGALS9, PDCD1, TNFRSF18, TNFRSF9, and VSIR were identified to be positively correlated with hepcidin expression in most of the tumors studied.

In the correlation analyses of the immune response pathways, we found positive correlations between hepcidin expression and almost all of the major immune response pathways, with no negative correlations identified. Common malignancies of the male genitourinary system in TCGA were all associated with the aforementioned findings (Figure 6C).

An immune score model was proposed as a prognostic tool for the evaluation of cancers. In this study, we found that the monocytic cell lineage was positively correlated with hepcidin across most of the tumors studied. Hepcidin seemed to be most associated, however, in PRAD and KIRC. In these tumors, lymphocytes and neutrophils were more likely to be correlated with hepcidin expression (Figure 6D). To further study the correlation between patient prognosis and immune cell infiltration in tumor tissues, assessments of patient prognosis were performed using an immune score model. Data extracted from TCGA database provided the cellular content of B cells, CD4⁺ memory T cells, CD8⁺ T cells, macrophages, mesenchymal stem cells, natural killer T cells, regulatory T cells, type-1 T-helper cells, and type-2 T-helper cells in BLCA, KIRC, KIRP, and TGCT tissues. The results revealed that in KIRC, the association between hepcidin and infiltration of immune cells had a significant impact on OS, with high expression indicating worse survival overall. In this model, KIRC patients who had enriched counts of the following were found: macrophages, CD4⁺ memory T cells, CD8⁺ T cells, and regulatory T cells. Meanwhile, there were also patients who had the number of macrophages, natural killer T cells, B cells, CD4⁺ memory T cells, CD8⁺ T cells, mesenchymal stem cells, type 1 T-helper cells, type 2 T-helper cells, and regulatory T cells decreased. All these were found in at least 90 cases, with all p-values <0.05. Specifically, macrophage enriched [cases included in analysis (n) = 434, HR: 1.79, 95% CI: 1.28–2.52, p = 0.00056], macrophage decreased (n = 93, HR: 3.26, 95% CI: 1.38–7.67, p = 0.0042), natural killer T cells decreased (n = 401, HR: 2.45, 95% CI: 1.62–3.69, P = 1e-05), B cells decreased (n = 285, HR: 2.43, 95% CI: 1.49–3.96, p = 0.00025), CD4⁺ memory T cells enriched (n = 348, HR: 1.64, 95% CI: 1.14–2.38, p = 0.0076), CD4⁺ memory T cells decreased (n = 179, HR: 2.91, 95% CI: 1.63–5.19, p = 0.00015), CD8⁺ T cells decreased (n = 143, HR: 3.15, 95% CI: 1.57–6.31, p = 0.00063), CD8⁺ T cells enriched (n = 384, HR: 1.67, 95% CI: 1.18–2.37, p = 0.0037), mesenchymal stem cells decreased (n = 491, HR: 2.25, 95% CI: 1.59–3.81, p = 2.6e-06), type 1 T-helper cells decreased (n = 483, HR: 2.1, 95% CI: 1.5–2.94, p = 8.6e-06), type 2 T-helper cells decreased (n = 476, HR: 1.9, 95% CI: 1.35–2.69, p = 0.00019), regulatory T cells decreased (n = 329, HR: 1.5, 95% CI: 1.04–2.17, p = 0.031), and regulatory T cells enriched (n = 198, HR: 2.09, 95% CI: 1.18–3.69, p = 0.0097). KIRC patients all showed increased hepcidin expression and were correlated with poor OS. A high expression of hepcidin was also linked to poor OS in patients with mesenchymal stem cell-decreased BLCA (n = 117, HR: 1.93, 95% CI: 1.07–3.51, p = 0.027) and to poor RFS in patients with regulatory T-cell-enriched KIRP (n = 72, HR: 4.45, 95% CI: 1.18–16.84, p = 0.017). The impact of hepcidin expression was thus more pronounced in patients with KIRC, and a higher expression of hepcidin may predict a worse prognosis in patients with KIRC.

Screening for potential hepcidin-binding proteins and genes associated with hepcidin expression was conducted to investigate the potential mechanisms involved in the effects of hepcidin in oncogenesis. As analyzed using the STRING database, a total of seven hepcidin-binding proteins were extracted. The interaction network between hepcidin and these seven proteins is shown in Figure 7A. GEPIA2 was then used to combine the expression data of the selected tumors from TCGA with the top 100 extracted genes associated with hepcidin levels. Positive correlations were found between the hepcidin gene and ALOX5AP (R = 0.48), LY86 (R = 0.47), PRAM1 (R = 0.55), FCGR2B (R = 0.5), and TREM2 (R = 0.61) (Figure 7B). As demonstrated in the heatmap of these genes (Figure 7C), they were all found to be positively correlated with hepcidin expression in common male genitourinary system tumors.

The KEGG and GO enrichment analyses were then conducted using the combined results of the aforementioned studies. In the GO analysis, the hepcidin-associated genes were found to be involved in the innate immune response, neutrophil degranulation, inflammatory response, IgG binding, integral component of the plasma membrane, regulation of immune response, phagocytosis, engulfment, G-protein-coupled adenosine receptor activity, and plasma membrane (Figure 7D). In the KEGG analysis, the top 10 signaling pathways were as follows: Staphylococcus aureus infection, osteoclast differentiation, lysosome, complement and coagulation cascades, pertussis, hematopoietic cell lineage, phagosome, acute myeloid leukemia, and leishmaniasis (Figure 7E). Detailed results are summarized in Supplementary Table S1. Based on the GO and KEGG analyses, the function of hepcidin-associated genes is mainly focused on the innate immunity response. However, we also noticed that several signaling pathways related to oncogenesis have been identified, including transcriptional misregulation in cancer (p = 0.012), cell adhesion molecules (CAMs) (p = 0.021), and the Toll-like receptor signaling pathway (p = 0.012), all of which may contribute to the oncogenesis of cancer (Figure 7E). These might explain why hepcidin-associated immune cells play a critical role in male genitourinary system tumors.

To verify the correlation between hepcidin and immune checkpoint molecule in Figures 6B-F,IF staining was carried out in FFPE KIRC tissues. Checkpoint molecules such as CD28, CD48, LGALS9, CTLA4, and PDCD1 were selected for validation. These checkpoint molecules were previously found to be positively correlated with the hepcidin expression in not only KIRC but also other common male genitourinary system tumors through bioinformatics analyses. As shown in Figure 8, the results of IF in FFPE KIRC tissues showed that the expression of hepcidin was weak; meanwhile, all the selected checkpoint molecules expressed less either. This result indicated that the expression of hepcidin and these checkpoint molecules were positively associated in KIRC tissue, in line with the conclusion of bioinformatics analyses.