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

Background The thyroid cancer subtype that occurs more frequently is papillary thyroid carcinoma (PTC). Despite a good surgical outcome, treatment with traditional antitumor therapy does not offer ideal results for patients with radioiodine resistance, recurrence, and metastasis. The evidence for the connection between iron metabolism imbalance and cancer development and oncogenesis is growing. Nevertheless, the iron metabolism impact on PTC prognosis is still indefinite. Methods Herein, we acquired the medical data and gene expression of individuals with PTC from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database. Typically, three predictive iron metabolism-related genes (IMRGs) were examined and employed to build a risk score (RS) model via the least absolute shrinkage and selection operator (LASSO) regression, univariate Cox, and differential gene expression analyses. Then we analyzed somatic mutation and immune cell infiltration among RS groups. We also validated the prognostic value of two IMRGs (SFXN3 and TFR2) by verifying their biological function through in vitro experiments. Results Based on RS, all patients with PTC were stratified into low- and high-risk groups, where Kaplan-Meier analysis indicated that disease-free survival (DFS) in the high-risk group was much lower than in the low-risk group (P < 0.0001). According to ROC analysis, the RS model successfully predicted the 1-, 3-, and 5-year DFS of individuals with PTC. Additionally, in the TCGA cohort, a nomogram model with RS was developed and exhibited a strong capability to anticipate PTC patients’ DFS. In the high-risk group, the enriched pathological processes and signaling mechanisms were detected utilizing the gene set enrichment analysis (GSEA). Moreover, the high-risk group had a significantly higher level of BRAF mutation frequency, tumor mutation burden, and immune cell infiltration than the low-risk group. In vitro experiments found that silencing SFXN3 or TFR2 significantly reduced cell viability. Conclusion Collectively, our predictive model depended on IMRGs in PTC, which could be potentially utilized to predict the PTC patients’ prognosis, schedule follow-up plans, and provide potential targets against PTC.


INTRODUCTION
The most common endocrine malignancy identified is thyroid carcinoma (TC), and recently its frequency has been steadily rising (Ferrari et al., 2020). Papillary thyroid carcinoma (PTC), the most prevalent pathological TC type, accounts for 85% of cases that are detected (Lim et al., 2017). Although PTC has a good overall prognosis, it has been reported that between 5% and 21% of cases have a recurrent tumor, lymph node metastasis, or distant postoperative metastasis, resulting in a significant reduction in survival rate (Nixon et al., 2014). Therefore, it is crucial to clarify the pathogenesis of PTC and detect effective prognostic markers.
For cells to keep functioning and maintain homeostasis, iron is a crucial component. Tumor incidence, progression, and metastasis are closely correlated with iron metabolism imbalance (Fischer-Fodor et al., 2015;Leftin et al., 2017;Steegmann-Olmedillas, 2011). Markedly, it is informed that iron metabolism has two-fold impacts on cancer cells. On one side, iron addiction refers to the phenomenon in which tumor cells, compared to healthy cells, become more dependent on iron intake for development and are more prone to deplete iron (Dev & Babitt, 2017). Alternatively, iron overload in tumor cells causes a unique type of programmed cell death named ferroptosis, where lipid peroxidation and excessive reactive oxygen species (ROS) cause cell death (Chen et al., 2021;Dixon et al., 2012). Several studies have identified that iron metabolism-related genes (IMRGs) are important predictors of malignancy prognosis. Nevertheless, the connection between IMRGs and PTC patients' prognosis has not been well examined (Li et al., 2020;Xu et al., 2021;Yuan, Liu & Zhang, 2021).
In this investigation, IMRGs were investigated in PTC. Depending on transcriptome and medical data available from the public database, we conducted extensive bioinformatics studies. We first defined the iron metabolism-related differentially expressed genes (IMR-DEGs) and conducted Pearson's correlation and network analyses. Subsequently, we developed three predictive IMRG signatures, assessed them, and verified them using the least absolute shrinkage and selection operator regression (LASSO, Tibshirani) regression and Cox regression analyses. Using this information as a foundation, we established a risk score (RS) system for PTC, performed multivariate and univariate Cox regression analyses, and then built and assessed a nomogram to anticipate the prognosis. We also assessed somatic mutation and immune cell infiltration among RS groups. Finally, we performed in vitro experiments to validate the function of two candidate genes.

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.

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: 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.

Cell viability assay
For cell viability assays, cells (2 × 10 3 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.

RESULTS
To visually and systematically describe our work, we presented the research process in Fig. 1.

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:  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 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 (cindex = 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).

DISCUSSION
Although most PTC patients have ideal postoperative effects and good prognoses, some well-differentiated TC tumors are still more invasive (Wen et al., 2021). They are resistant to standard treatments, such as those that are non-operable, relapse after surgery, and do not respond to radioiodine treatment (Qin et al., 2021). Additionally, PTC is a very heterogeneous condition, and the development of tumors includes a complicated network made up of several signaling mechanisms (Balachandran et al., 2015;Zeiger & Schneider, 2013). Therefore, it is particularly important to detect markers related to tumor diagnosis and prognosis to formulate treatment and follow-up plans.
Iron is vital for cell viability as it exists in proteins that perform diverse functions, including respiration, transport, and homeostasis of oxygen, as well as the synthesis of biomolecules (Dlouhy & Outten, 2013). Moreover, iron acts as a crucial component of many proteins, including the repair and metabolism of nucleic acid as well as the progression of the cell cycle (Zhang, 2014). Since iron is an essential part of physiology and  anatomy with low bioavailability, the human body regulates iron stores strictly to maintain conservation and minimize toxic effects (Dev & Babitt, 2017). The iron metabolism imbalance contributes to many diseases, and it is significant to explore this role in the search for therapeutics (Kell, 2009). Many recent studies have demonstrated that iron metabolism has a function in all cancer development stages (Guo et al., 2021;Recalcati, Gammella & Cairo, 2019). The total disturbance of iron absorption, effusion, storage, and modulation in malignancies suggests that reprogramming iron metabolism will result in the dysregulation of proliferation and survival in tumor cells (Andrews, 2008;Jung et al., 2019;Manz et al., 2016;Wang et al., 2018). Recent investigations have shown that vitamin C activates ferritinophagy to cause ferroptosis in anaplastic thyroid tumor cells . Lin et al. (2021) also used network analysis to examine the prognosis of PTC and ferroptosis genes in immune infiltration. The dynamic profile of genes associated with iron metabolism in TC is still unknown.
In this investigation, to screen seven down-regulated and five up-regulated IMR-DEGs, we employed clinicopathological data and gene expression from the open-access database. Network analyses indicated that these genes were associated with multiple transcription factors, miRNAs, small molecular compounds, and drugs, which may reveal their role in tumor progression and potential as therapeutic targets. Subsequently, we used Cox and LASSO regression analyses to create three predictive IMRGs signatures and established an RS system to construct and evaluate a nomogram for predicting prognosis.
Using LASSO regression and univariate Cox analyses, a total of three IMR-DEGs were found to be possible predictive biomarkers, and these genes were utilized to build a predictive model. Two of them (SFXN3, TFR2) had negative correlations between their expression levels and DFS, but the expression level of the LTF gene correlated positively with DFS. During one-carbon metabolism, SFXN3 is a crucial mitochondrial serine transporter that is involved in tumor cell growth . Recent research has shown that the immunosuppressive microenvironment was linked to the elevated SFXN3 expression mediated by non-coding RNA, which was a predictive marker in head and neck squamous cell cancer (HNSC) . Through performing a competitive endogenous RNA analysis, Zheng et al. (2022) also showed that MIR193HG-miR-29c-3p-SFXN3 participated in HNSC, which significantly influenced the treatment efficacy and prognosis. In addition, Murase et al. (2008) indicated that the serum SFXN3 autoantibodies may act as a new cancer biomarker for oral squamous cell cancer. Consistent with these studies, our research also indicated that the elevated SFXN3 expression represented a poor predictive marker and an oncogene in PTC, and its silence significantly inhibited the proliferation of PTC cells.
Iron is transported to cells from peripheral blood through a transmembrane glycoprotein receptor family called the transferrin receptor family. TFR2, a subtype of TFRs, has been found altered expression in tumor cells. An elevated TFR2 expression was identified in the colon, glioblastoma (GBM), and ovarian cancer cell lines (Calzolari et al., 2009;Calzolari et al., 2010). However, in subjects with myelodysplastic syndrome, acute myeloid leukemia, and GBM, the higher transcription levels of TFR2 were connected to a better prognosis than that in the lower group (Calzolari et al., 2010;Di Savino et al., 2017;Nakamaki et al., 2004). The phenomenon might be due to the fact that TFR2 sensitizes tumor cells to cell-cycle-specific chemotherapy medications by regulating potential cellular signaling . For instance, the increased TFR2 expression in GBM cells led to a high increase in growth, which was highly sensitive to temozolomide (Calzolari et al., 2010). In our work, we illustrated that the increased TFR2 expression is related to a worse prognosis and proliferation of PTC, and this phenomenon is different from that of other tumors. Exactly, it will be of potential value to explore the internal mechanism behind these differences in the future.
LTF is an iron-binding transport glycoprotein identified as a key immune-related gene correlated with PTC prognosis Qin et al., 2021). LTF expression level in cancer tissue was markedly lower than that in healthy tissue. Also, ssGSEA outcomes indicated that LTF expression level was strictly connected to the immune process (Qin et al., 2021). This is consistent with current results that LTF was downregulated in PTC and was related to a better prognosis. Furthermore, Hosseinkhan et al. (2020) determined LTF as a promising biomarker down-regulated in a stage I subgroup and in the common stage IV of PTC. Taken together, this evidence suggests that these genes may have a vital function in the incidence and progression of malignancies and might be ideal prognostic markers of PTC prognosis.
The multigene signature for DFS was an independent predictive variable in PTC patients in subsequent studies. Patients with high risk possessed a worse DFS than those with low risk, according to risk classification by RS. This model was successfully and consistently verified in several patient populations, and multivariate Cox regression analysis validated its independence as a predictive marker. To further predict the DFS of individuals with PTC, we built a predictive nomogram model depending on IMRGs, which included the T and N stages in addition to RS. The calibration curves also demonstrated the nomogram's reliable prognostic value for DFS in the TCGA cohort. This nomogram model could be utilized to examine PTC patients' prognoses and schedule follow-up plans. Interestingly, we observed that compared to the low-risk group, the high-risk group had significantly higher BRAF mutation frequency, TMB, and some immune cell infiltration. BRAF mutation is a common event in PTC, and its presence indicates tumor progression (Guerra et al., 2014;Kim et al., 2014). Moreover, the amount of tumor somatic coding mutations is defined as TMB, which is used to estimate neoantigen burden and to predict immune checkpoint inhibition response across a variety of tumor types (Sha et al., 2020). Overall, these findings indicated the clinical significance of RS grouping beyond predicting prognosis.
The clinical samples were subsequently validated by IHC, which revealed elevated SFXN3 and TFR2 levels in cancerous tissue. It is well known that iron is essential for ribonucleotide reductase's catalytic activity, an enzyme responsible for converting ribonucleotides into deoxyribonucleotides, the rate-limiting stage of DNA synthesis as well as a necessary step in cell replication (Torti & Torti, 2020). Massive hepcidin release is a TC biomarker that causes reduced iron exporter, FPN, expression, and elevated intracellular iron retention, thus promoting tumor growth (Zhou et al., 2018). Therefore, we conducted the cell viability assay to determine the silencing IMRGs' effect on the cell proliferation ability of PTC cell line TPC-1. The outcomes exhibited that both SFXN3 and TFR2 knockdown suppressed the TC growth in vitro. Notably, in our study, SFXN3 and TFR2 could affect iron metabolism through the regulation of FTL, FTH, and TFR expression in iron-related proteins. FTL and FTH belong to ferritin, while TFR belongs to the iron importer. After silencing SFXN3 or TFR2, FTL, and FTH were down-regulated, while TFR was up-regulated, all of which could up-regulate intracellular iron concentration. These results suggested that the silence of SFXN3 or TFR2 may lead to the overload of intracellular Fe 2+ , thereby causing ferroptosis Lu et al., 2021). Collectively, these investigations also suggest that SFXN3 and TFR2 may be viable therapeutic targets and markers of clinical prognosis in PTC.
However, there are still restrictions on our work. First, this was a retrospective work with bias in selecting variables such as clinicopathological features. Therefore, there is still a need to improve the accuracy of the data. Second, our prediction model depends on the survival function estimation after analyzing multiple influencing factors in TCGA comprehensively. Hence, there are limited clinical implications from this study. It is essential to obtain more prospective information to validate its clinical significance. Finally, although we revealed the characteristics of three IMRGs, it is still crucial to carry out further trials to clarify the roles and pathways in PTC development.

CONCLUSIONS
In conclusion, this study provided, as far as we know, the first application of IMRGs to evaluate their prognostic prediction ability in PTC. We examined IMRGs in PTC and found three genes linked to the predictive and clinicopathological features of PTC, of which two were validated in 40 PTC and 38 NT samples by IHC, and conducted in vitro experiments to explore the biological function. Additionally, we created and verified an RS system for risk categorization and prognosis. Finally, a nomogram model was built that demonstrated high predictive accuracy for 1-, 3-and 5-year DFS rate predictions. These findings indicate that IMRGs possess the predictive ability of PTC prognosis, and our prognostic model could be included in large-scale prospective research in the future to further verify its clinical value.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
Our research was funded by the National Natural Science Foundation of China under Grant (U20A20382). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: National Natural Science Foundation of China: U20A20382.