Identification of iron metabolism-related genes as prognostic indicators for papillary thyroid carcinoma: a retrospective study

Data download and preprocessing

The Cancer Genome Atlas (TCGA) database was employed to acquire the transcriptome data (RNA-seq) of 494 PTC and 59 normal tissue (NT) samples in conjunction with medical features and follow-up outcome data. Consequently, the tumor group was randomly distributed into model training and validation sets at a 7:3 ratio. Finally, 350 cases of tumor tissue were in the training set, and 144 cases of tumor tissue were in the validation set. Subsequent analyses and prognostic model construction used a training set and validated it in a validation set. At the same time, we downloaded the original microarray outcomes (CEL format) of transcriptome dataset GSE33630 of PTC and NT samples with reliable sample sources from the Gene Expression Omnibus (GEO) database (Dom et al., 2012), as well as the annotation file of GPL570[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform. Transcriptome data for 45 NT and 49 PTC were included in the GSE33630 collection.

Expression and correlation analysis of IMRGs

IMRGs were from two databases: Iron uptake and transport (R-HAS-917937) pathway genes from the Reactome database (https://reactome.org/) and cellular iron ion homeostasis genes (GO:0006879) from the AmiGo2 database (http://amigo.geneontology.org/amigo). Finally, 70 IMRGs were screened (Table 1). The IMR-DEGs were identified utilizing the limma package of R language (Ritchie et al., 2015) using a linear model. Adj. P < 0.05 and |log2FC| > 1 were the criteria for DEG screening. We used the ggplot2 package (Villanueva & Chen, 2019) and pheatmap package (Kolde, 2019) to the heatmaps and volcano plots of differentially expressed genes (DEGs) to illustrate the differential expression analysis outcomes. Pearson’s correlation was then analyzed between different IMR-DEGs.

Protein-protein interaction (PPI) network analysis and transcription factor, miRNA, drug-protein interaction network analysis

For IMR-DEGs, we built a PPI network utilizing the STRING (http://www.string-db.org/, Version: 11.0) online method (Szklarczyk et al., 2019) and Cytoscape program (Otasek et al., 2019). We further analyzed transcription factors, miRNA, small molecular compounds, and their drug interaction networks. The TF-gene, Gene-miRNA, Protein-chemical, and Protein-drug interactions modules of NetworkAnalyst were used to examine the uploaded DEGs (https://www.networkanalyst.ca/) (Xia, Gill & Hancock, 2015), respectively. TF databases, ENCODE (http://cistrome.org/BETA/), miRTarBase v8.0 (https://mirtarbase.cuhk.edu.cn), comparative toxicogenomics database (CTD, http://ctdbase.org/) and DrugBank database v5.0 (https://go.drugbank.com/) were the reference database. The analyzed network was visualized by Cytoscape software.

Screening of prognostic markers

For the screened IMR-DEGs, we used Cox regression analysis to assess the connection between disease-free survival (DFS) of individuals with PTC and gene expression. A compression estimation is the LASSO approach (Tibshirani, 1996). By designing a penalty function, which causes it to compress certain coefficients and set others to zero, the model is refined. As a result, it still has the benefit of subset shrinking and is biased when processing data that have complicated collinearity. Considering the findings of survival analysis, we employed LASSO regression to further identify predictive biomarkers. The variables were screened using the glmnet function of the glmnet package (Engebretsen & Bohlin, 2019) and cross-validated using the cv. glmnet function, and then the combination of predictive biomarkers with the lowest CV coefficient was selected.

Construction of RS and evaluation of clinical prognosis predictive ability

The following formula was employed to detect RS for each case:

$RS={\sum }_{i=1}^n\mathrm{C}\mathrm{o}\mathrm{e}{\mathrm{f}}_i\times Ex{p}_i,$where Coef was the LASSO regression coefficient and Exp was the gene expression level (log2 conversion).

The best RS cutoff value for predicting the survival time of individuals with PTC was estimated utilizing the maxstat package (Hothorn, 2017). Subsequently, according to the cutoff value, participants were divided into low- and high-risk groups, and the survival curve was plotted utilizing the Kaplan-Meier technique. Simultaneously, the survival period of 1-, 3- and 5-year for individuals with RS was predicted utilizing the survivalROC package (Heagerty, 2013). Then the anticipated ROC was mapped, and the AUC value was estimated.

We used a Cox proportional hazards model to evaluate the other medical features’ impact on patients’ prognostic survival, including age, sex, tumor stage, and TNM stage, and developed forest maps with the forestmodel package (Kennedy, 2020). Subsequently, to evaluate the independent RS predictive value on patient prognosis, the clinical parameters that significantly affected prognosis were included in multivariate Cox regression as covariates. A forest map was then created. The fitting effect of different models was evaluated by AIC values.

Finally, to display the model outcomes and make the prediction model results more readable, the nomogram and calibration curve of the best multifactor model was created utilizing the rms package (Harrell, 2021). The nomograms’ predictive ability for survival was estimated through the calculation of the concordance index (C-index).

Differential expression analysis and functional enrichment analysis of risk groups

To investigate the potential RS biological importance, we conducted a differential study of gene expression in various risk groups. Adj. P < 0.05 and |log2FC| > 1 were the DEG examining criteria. The heatmap for the DEGs and the volcano plot were created utilizing ggplot2 (Villanueva & Chen, 2019) and pheatmap packages (Kolde, 2019) to show the differential analysis outcomes.

A famous database that includes information on biological mechanisms, genomes, diseases, and medicines is called the Kyoto Encyclopedia of Genes and Genomes (KEGG). Go function annotation analysis, which includes molecular function (MF), cell component (CC), and biological process (BP), is a typical approach for large-scale gene function enrichment study. KEGG enrichment analyses and gene ontology (GO) (Harris et al., 2004; Kanehisa & Goto, 2000) were conducted on DEGs of risk groups utilizing the clusterProfiler package (Yu et al., 2012); P < 0.05 was deemed statistically significant.

Analysis of somatic mutation and immune cell infiltration

Fisher’s exact test was performed on all genes in the high-/low-risk group of the TCGA dataset to detect differentially mutated genes using the “mafCompare” function in the R package Maftools (version: 2.14.0) (Mayakonda et al., 2018). Subsequently, the Wilcoxon test was used to compare tumor mutation burden (TMB) between high- and low-risk groups, and a violin plot was used to present the results.

We used the gene expression matrix from the TCGA database in the CIBERSORTx (https://cibersortx.stanford.edu/) (Newman et al., 2015) online analysis tool to calculate the immune cell infiltration of the samples and filtered out samples with P < 0.05.

The patients and the tissue samples

We collected PTC samples retrospectively from patients, who had undergone one-stage surgery to analyze the expression of SFXN3 or TFR2. A written informed consent form was completed by each patient, and all studies were approved by the ethics committee of Zhejiang provincial people’s hospital (Ethical Approval No: 2021QT251).

Immunohistochemistry (IHC)

Xylene was employed to deparaffinized paraffin-embedded thyroid segments, and ethanol was utilized to rehydrate it. The antigen was then recovered from the segments utilizing 1 mM EDTA, and non-specific binding was reduced by preincubating the segments in TBS with 5% goat serum. The sections were then treated with either the TFR2 (A9845; ABclonal, Wuhan, China) or SFXN3 (15156-1-AP; Proteintech, San Diego, CA, USA) antibodies. Thyroid tissue samples were examined using DAB and counterstained with hematoxylin after treatment with horseradish peroxidase secondary antibodies.

Every specimen’s immune-positive rate and staining levels were incorporated into the IHC score. Three pathologists examined all the slices’ protein-expressing scores in the current investigation. SFXN3 and TFR2’s immunoreactivity score (IRS), which ranged from 0 to 12, was determined as the intensity and positive rate product (Song et al., 2021).

Cell culture

The PTC cell line TPC-1 was given by Dr. Xin Zhu (Zhejiang Cancer Hospital). It was cultured in RPMI-1640 + 10% FBS. Short tandem repeat (STR) profiling has been employed to regularly detect mycoplasma in all cell lines and to verify their authenticity.

Transfection

SiRNAs (small interfering RNA) targeting SFXN3 and TFR2 were obtained from RiboBio (RiboBio Biotechnology, Guangzhou, China), the sequence of SFXN3 siRNA#1 was CTGGAAGCTTCTCGGAACA, the sequence of SFXN3 siRNA#2 was GCGTTGAAGTGGTCTACTA, the sequence of SFXN3 siRNA#3 was GGTGAATTGCCTTTAGACA, the sequence of TFR2 siRNA#1 was GGCTAGTGGTCAACAATCA, the sequence of TFR2 siRNA#2 was CTCAAGGAGTGCTCATATA, and the sequence of TFR2 siRNA#3 was TGGGCAGACTCTCTATGAA. We transfected SFXN3-siRNAs or TFR2-siRNAs using jetPRIME (Polyplus, Berkeley, CA, USA) at final doses of 50 nM after cells had been grown in a six-well plate at 40% confluence. After 2 days, cells were obtained for further trials.

qRT-PCR and western blots (WB)

We performed qRT-PCR and WB as previously described (Pan et al., 2022). Table 2 displayed the primers used for each gene.

For WB, the primary antibodies included: ferritin light chain (FTL, 10727-1-AP; Proteintech, San Diego, CA, USA), ferritin heavy chain (FTH, YT1692; ImmunoWay, Jiangsu, China), transferrin receptor (TFR, YT5374; ImmunoWay, Jiangsu, China), GAPDH (10494-1-AP, Proteintech, San Diego, CA, USA).

Cell viability assay

For cell viability assays, cells (2 × 10³ cells/well) transfected with siRNA were incubated for two days in 96-well plates. Cultures were examined utilizing the Cell Counting Kit-8 (CCK-8; Beyotime, Shanghai, China) following the manufacturer’s directions.

Statistical analyses

The R program (R Core Team, 2022) was employed to perform all statistical analysis and data calculations. The Benjamini-Hochberg (BH) method was utilized for multiple test corrections, and FDR correction was performed in multiple tests to decrease the false positive rate. For the comparison of two groups, the Mann-Whitney U test (i.e., Wilcoxon rank sum test) was employed to compare the differences of continuous variables with non-normal distributions, while the independent student T-test was utilized to determine the statistical significance of continuous variables with normal distributions. To evaluate the prognostic biomarkers’ predictive value, a Cox regression model was used. The pROC package of R was employed to plot the receiver operating characteristic (ROC) curve, and the accuracy of RS in determining prognosis was evaluated by computing the area under the curve (AUC). P-values for all tests were two-sided, with P < 0.05 indicating statistical significance.

Differential expression analysis and correlation analysis of IMRGs

The 70 IMRGs expression levels in the tumor and control groups were compared. Finally, we discovered that the 12 IMRGs expression had significant differences, among which SFXN3, TFR2, TMPRSS6, HMOX1, and LCN2 expressions were elevated in the cancer group, and SCARA5, LTF, STEAP2, SLC39A14, CP, ALAS2, and TF were low expressed in the tumor group (Fig. 2A). The IMRGs expression levels in different samples were displayed by heatmap (Fig. 2B). Figure 2C illustrates the correlation matrix of IMRGs expression level. Among them, HMOX1 exhibited significant positive correlations with TMPRSS6, SFXN3, and LCN2 (P < 0.001); SFXN3 exhibited significant positive correlations with TMPRSS6, LCN2, and TFR2 (P < 0.001); LCN2 exhibited a significant positive correlation with TMPRSS6 (P < 0.001); SCARA5 exhibited a significant positive correlation with CP (P < 0.001). Each correlation coefficient was greater than or equal to 0.45. In contrast, SLC39A14 exhibited a significant negative correlation with TMPRSS6 (P < 0.001), with a correlation coefficient of less than −0.45 (Fig. 3).

Protein, transcription factor, miRNA, small molecule compound, and drug interaction networks of IMR-DEGs

PPI network analysis of IMR-DEGs was conducted in cancer and healthy tissues, and Fig. 4A presents the PPI network. In the PPI network, TFR2, CP, and SLC39A14 had a large weight and a strong connection. The transcription factor analysis showed that SFXN3 was related to 72 transcription factors, LCN2 was related to 45 transcription factors, HMOX1 was related to 42 transcription factors, SLC39A14 was related to 35 transcription factors, CP was related to 24 transcription factors, TF was related to 24 transcription factors, and STEAP2 was related to 11 transcription factors. The miRNA interaction network analysis showed that HMOX1 was related to 47 miRNAs, STEAP2 was related to 43 miRNAs, SLC39A14 was related to 22 miRNAs, SFXN3 was related to 20 miRNAs, LTF was related to 18 miRNAs, and TMPRSS6 was related to 14 miRNA. The results of small molecule compound analysis showed that HMOX1 was related to 545 compounds, TF was related to 89 compounds, CP was related to 80 compounds, LCN2 was related to 41 compounds, SLC39A14 was related to 39 compounds, LTF was related to 33 compounds, STEAP2 was related to 25 compounds, and SFXN3 was related to 22 compounds. The results of the drug interaction analysis showed that LCN2 was related to six drugs, HMOX1 was related to nine drugs, and ALAS2 was related to two drugs (Fig. 4).

Survival analysis of IMRGs and screening of prognostic markers

The connection between IMR-DEGs and the patient’s prognosis was examined utilizing univariate Cox regression, and it found that LTF, SFXN3, and TFR2 were significantly connected with prognosis (Fig. 5A). Subsequently, LASSO regression analysis was performed and three IMRGs were still retained, which could be used as joint prognostic markers (Figs. 5B–5D). The prognostic ROC curve for the LASSO regression revealed an AUC value of 0.718, which indicates good prediction ability (Fig. 5B). Figures 5E–5G show the survival curves for the three prognostic markers grouped by high-/low-expression.

RS and prognostic predictive model construction

The coefficients of candidate predictive biomarkers were detected based on the LASSO regression model outcomes, and the RS was detected with the following formula:

$\mathrm{R}\mathrm{S}=(-0.1201)\ast \mathrm{L}\mathrm{T}\mathrm{F}+(0.0053)\ast \mathrm{S}\mathrm{F}\mathrm{X}\mathrm{N}3+(0.4250)\ast \mathrm{T}\mathrm{F}\mathrm{R}2.$

Figure 6A illustrates the ROC curve indicated by the score for the 1-, 3- and 5-year survival. Among them, the predictive ability for the 3-year survival was the best (AUC = 0.688). We then used the maxstat package to determine that the best RS threshold value for anticipating the survival time of individuals with PTC was 1.2779 (Fig. 6B). According to the cutoff value, we distributed the cases into high- and low-risk groups and patients with increased RS had a significantly lower prognostic survival time than those with low RS (Fig. 6C). The univariate COX regression showed that besides the RS/grouping, cancer, N, and T stages had an influence on the patient’s survival (Fig. 6D).

Therefore, in the validation set, we validated the predictive biomarkers examined by LASSO regression. The TCGA validation set indicated that differential expression was observed between tumor and control groups for the three prognostic markers, and the GEO validation set also verified this. Moreover, the TCGA validation set and the training set’s survival analysis for risk grouping showed agreement (Figs. 7A–7C).

We compared the differences in candidate prognostic markers in various clinical subgroups. The outcomes indicated that there were LTF expression variations in subgroups of tumor stage and N and T stages, SFXN3 in subgroups of tumor stage and M, N, T stage, and TFR2 in subgroups of M, T stage and tumor stage (Figs. 8A–8C).

Subsequently, a multivariate prognostic model was established. Due to the correlation between the tumor stage and the TNM stage, we used a Cox regression model to build a multivariate prognostic prediction model according to the tumor, T, and N stages and RS. Figures 9A and 9B illustrate the forest map. The AIC value of the T and N stages, in addition to the RS model, was 487, while the AIC value of the tumor stage and RS model was 488. Therefore, the model of T and N stages, in addition to RS, was finally selected as the best model. Simultaneously, the nomogram (Fig. 9C) and calibration curve (Figs. 9D–9F) were drawn with an rms package to anticipate the 1-, 3- and 5-year survival possibility of individuals with PTC. The nomogram created with T and N stages (c-index = 0.67) and RS (c-index = 0.68) separately displayed that these variables could be good markers. However, the combined survival prediction model’s c-index was 0.74, indicating a higher capacity for prognosis prediction.

DEGs analysis and functional enrichment analysis in RS grouping

According to the high-/low- RS grouping, the Limma package of R language was employed to screen 465 DEGs from the TCGA mRNA expression matrix, including 289 mRNA down-regulated and 176 mRNA up-regulated. Figures 10A and 10B present the heatmaps of the DEGs and the volcano plot in the TCGA dataset.

KEGG analysis revealed that the mechanisms enriched by DEGs mainly included mineral absorption, tyrosine metabolism, and neuroactive ligand-receptor interaction (Figs. 10C–10F). The statistical outcome was shown in Table 3. GO analysis suggested that DEGs were mainly connected to BP, such as signal release, hormone transport, and regulation of membrane potential; CC, such as cation channel complex, synaptic membrane, and ion channel complex; and MF, such as ion channel, channel, and passive transmembrane transporter activities (Figs. 10G and 10H). The statistical outcome was shown in Table 4.

Analysis of somatic mutation and immune cell infiltration

We subsequently analyzed the difference in somatic mutation among risk groups and the distribution of TMB. Compared to the low-risk group, the high-risk group had a significantly higher mutation proportion of BRAF (85% vs 56%) (Fig. 11A). In addition, compared to the low-risk group, the high-risk group also had a slightly higher TMB, and the significance level is at the critical value (P = 0.054) (Fig. 11B). According to the analysis of immune cell infiltration, the high-risk group exhibited a significantly higher degree of infiltration of immune cells, such as T cells regulatory, macrophages M0, and dendritic cells activated than the low-risk group (P < 0.05) (Fig. 11C).

The correlation analysis between the immune cell infiltration degree and the candidate genes’ expression levels displayed that LTF was positively correlated with macrophages M1, T cells follicular helper, B cells memory, and B cells naïve, respectively (P < 0.001); SFXN3 was positively correlated with dendritic cells resting, macrophages M0, T cells regulatory, and T cells CD4 memory resting, respectively (P < 0.001), while it was negatively correlated with NK cells activated, T cells follicular helper, T cells CD8, and plasma cells, respectively (P < 0.001); TFR2 was positively correlated with T cells regulatory (P < 0.001), and negatively correlated with plasma cells (P < 0.001) (Fig. 11D).

Determination of IMRGs expression level and functional analysis in PTC

To determine the SFXN3 and TFR2 expression levels in PTC tissues, 40 PTC tissues and 38 non-tumorous thyroid tissues were detected, respectively. IHC (Fig. 12A) exhibited that the SFXN3 and TFR2 levels in PTC tissues were up-regulated. Next, to determine whether the elevated expression of SFXN3 and TFR2 was responsible for preserving the TC aggressive phenotype. SFXN3 and TFR2 were knocked down with high efficacy, confirmed at the mRNA levels, respectively (Fig. 12B). Proliferation was decreased in siSFXN3 or siTFR2 transfected TC cells, as demonstrated by a cell viability assay (Fig. 12C). These outcomes suggested that SFXN3 and TFR2 might function as oncogenes in PTC.

To explore the pathways of SFXN3 and TFR2 functioning in intracellular iron metabolism, a profile of iron-related protein mRNA expression was conducted, including ferritin (storage), ferroportin (FPN; iron exporter), divalent metal transporter 1(DMT1; iron importer) and TFR (iron importer). The results indicated that the silence of SFXN3 or TFR2 significantly increased TFR expression while decreasing FTL and FTH expression (Fig. 12D). Further WB experiments also confirmed this (Fig. 12E).