A new ferroptosis-related gene model for prognostic prediction of papillary thyroid carcinoma

ABSTRACT Papillary thyroid carcinoma (PTC) is a highly heterogeneous malignancy with diverse prognoses. Ferroptosis is a new type of cell death dependent on iron. Nevertheless, the predictive ability of ferroptosis-related genes for PTC is unclear. Based on the mRNA expression information from The Cancer Genome Atlas, we compared tumor and normal tissues in terms of the gene expression, for identifying differentially expressed genes (DEGs). Then, the risk score of a 5-gene signature was calculated and a prognostic model was established to test the predictive value of this gene signature by virtue of the LASSO Cox regression. The 5 genes were validated in PTC tissues by RT-qPCR.At last, functional analysis was implemented to investigate the underlying mechanisms. We found a total of 45 ferroptosis-related genes expressed differentially between tumor and normal tissues. 6 DEGs exhibited a significant relevance to the overall survival (OS) (P< 0.05). We classified patients into group with high risk and group with low risk based on the median risk score of a 5-gene signature. Patients in the group with low risk presented a remarkably higher OS relative to the group with high risk (P< 0.01). The Cox regression analysis displayed the independent predictive ability of the risk score. The receiver operating characteristic analysis helped to validate the predictive power owned by the gene signature. After validation, the 5 genes were abnormally expressed between PTC and normal tissues. Functional analysis showed two groups had different immune status. A new ferroptosis-related gene signature can predict the outcomes of PTC patients.


Introduction
Papillary thyroid carcinoma (PTC) is considered the most common endocrine malignancy with multiple prognoses. The incidence of PTC has steadily increased over the last few years in many countries [1]. The survival rates of PTC patients vary greatly. While most PTC patients had an optimistic prognosis with a 10-year survival rate of 80-90%, some only had survival rates of 25% at 5-year and 10% at 10-year [2,3]. Therefore, it is necessary to develop novel predictive models for PTC.
Ferroptosis -characterized by lipid peroxidation and iron accumulation -is a form of programmed cell death [4,5]. Regulating the ferroptosis pathway has the potential to slow down cancer progression [6,7]. Studies showed that some genes, such as SLC7A11 [8], GPX4 [9] and p53 [10], regulated ferroptosis in cancer cells and affected the prognosis of hepatocellular carcinoma [11] and breast cancer [12]. However, it remains unclear whether these ferroptosis-related genes have an effect on the prognosis of PTC.
The purpose of this study was to explore the effect on the prognosis of PTC of ferroptosisrelated genes and facilitate the development of prognostic models for PTC patients. We obtained PTC patients' clinical data as well as mRNA expression information from The Cancer Genome Atlas (TCGA). The ferroptosis-related differentially expressed genes (DEGs) were used for constructing a prognostic model. Then, we validated the model. Finally, we explored the underlying mechanisms using functional enrichment analysis.

Data collection
The clinical data as well as the mRNA expression information regarding 507 PTC patients were obtained from the TCGA database (https://portal. gdc.cancer.gov/, up to December 2020). We used the 'limma' R package to normalize gene expression data. This study was exempted from ethics reviews since all data we used were publicly available, and we followed the TCGA Ethics & Policies.

Construction and validation of a signature with ferroptosis-related genes
We compared gene expression levels between tumor and adjacent normal tissues to identify the DEGs using the 'limma' R package [16]. The criterion of DEGs was a false discovery rate (FDR) < 0.05. Cox analysis assisted in examining the ability of ferroptosis-related genes to predict overall survival (OS). Benjamini-Hochberg adjusted p-values were used to decrease FDR. The 'glmnet' R package served for a LASSO Cox regression -a powerful technique for variable selection and regularization -to analyze if DEGs could predict the OS and the status of the PTC patients [17]. The optimal value of the penalty parameter (λ), which corresponds to the minimum of the partial likelihood deviance, was identified by ten-fold cross-validation. We calculated the risk score using the following formula: risk score = esum (the normalized expression level regarding each gene× its regression coefficient) [11]. We then used the median risk score as the standard for classifying patients into group with high risk and group with low risk. We depicted the gene distribution in two groups by performing PCA and t-SNE using the 'stats' and the 'Rtsne' R package [18]. The 'surv_cutpoint' function of the 'survminer' R package served for the survival analysis to identify the optimal cutoff values for gene expression [19]. We then used the 'survivalROC' R package to perform time-dependent ROC analysis to estimate if the gene signature has predictive ability [20]. We then used the univariate and multivariate Cox regression analyses for determining if the risk score could independently predict patients' prognosis in terms of the OS.

Functional enrichment analysis
We conducted the Gene Ontology (GO) enrichment as well as the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the DEGs between two groups with the 'clusterProfiler' R package and |log2 (foldchange)| ≥ 1 and FDR <0.05 were considered as the criteria for the DEGs [21]. The 'gsva' R package was applied to the single-sample gene set enrichment analysis (ssGSEA) for calculating the enrichment score regarding immune cells as well as immune-related pathways [22].

Patients and Specimens
We selected 15 papillary thyroid carcinoma tissues and 10 normal tissues from the First Hospital of Jiaxing between 2020 and 2021. Postoperative pathological examination confirmed papillary thyroid carcinoma. This study was approved by the Ethics Committee of the First Hospital of Jiaxing. Each patient signed an informed consent. All specimens were placed in liquid nitrogen and stored at −80°C immediately.

Quantitative Real-Time PCR
Trizol extracts total RNA from tissues and then reverse transcripts it into cDNA. PCR was performed using TB Green®Premix Ex Taq™II kit (Takara, China). The reaction process is as follows: 94°C for 30 s, 58°C for 30s, 72°C for 60 s, 40 cycles. GAPDH as internal control. The relative expression level was determined by 2− ΔΔ CT.

Statistical analyses
Student's t-test served for comparing the gene expressions between tumor issues and adjacent normal tissues. Chi-square test assisted in examining the differences in proportions between group with high risk and group with low risk. We performed the Mann-Whitney test for comparing the enrichment scores regarding immune cells as well as immunerelated pathways between the two groups. The Kaplan-Meier method together with the log-rank test were conducted to depict the survival curves. We then conducted Cox regression analysis to figure out the factors that could predict OS. We performed all statistical analyses in R (Version 4.0.3). A two-tailed P value < 0.05 exhibited statistical significance.

Results
In our research, we set up a prognostic model by using ferroptosis-related DEGs to investigate the effect of ferroptosis-related genes on the prognosis of PTC, and then we verified the model by internal validation and tissues validation. At last, we explored the underlying mechanisms. The results were as follows: Baseline data regarding PTC patients A total of 507 PTC patients were selected from the TCGA cohort whose median age was 46 (range 15-89) years. Most subjects were female and were in stage Ι. Distant metastasis occurred in a few patients. They had the median OS of 714 (range 6-5423) days. Table 1 lists the detailed characteristics of the subjects.

Identification of ferroptosis-related DEGs with predictive value
We found 45 ferroptosis-related genes expressed differentially between tumor and adjacent normal tissues, 6 of which presented a significant association with OS ( Figure 1a, p < 0.05). The forest plots displayed the associations of gene expression with OS ( Figure 1b). The heatmap indicated tumor and adjacent normal tissues could be distinguished by the DEGs (Figure 1c). Figure 1d showed the associations between these genes.

Construction of a prognostic model
Aforementioned 6 genes were employed for establishing a prognostic model by virtue of LASSO Cox regression analysis, which were associated with OS. A 5-gene signature was then determined according to the optimal value of λ. We used the following formula to calculate the risk score: e ( . We used the median risk score as the standard for classifying patients into group with high risk and group with low risk (as shown in Figure 2a). Between the two groups, the stages of cancer were significantly different ( Table 2, P< 0.05). PCA and t-SNE analysis showed the patients with PTC were distributed in two directions between the two groups (Figure 2c-d). Patients in group with high risk exhibited a higher mortality rate relative to group with low risk (Figure 2b). Similarly, the Kaplan-Meier curve displayed that group with low risk had remarkably higher OS relative to group with high risk (Figure 2e, p< 0.05). Figure 2f showed the OS predictive power of the risk score, with the area under the curve (AUC) of 0.621, 0.728, and 0.875 at 1 year, 2 years, 3 years, respectively.

Independent predictive ability of the risk score
Cox regression analysis served for confirming if the risk score could predict OS independently. As found by the univariate Cox regression analysis, the risk score exhibited a significant association with OS (HR = 10.697, 95% CI = 1.328-86.173, P = 0.026) (Figure 3a). Consistently, the multivariate Cox regression analysis discovered a significant association between the risk score and OS (HR = 11.682, 95% CI = 1.454-93.878, P = 0.021) (Figure 3b).

Functional analysis in the TCGA
We conducted the GO enrichment as well as KEGG pathway analysis of the DEGs for identifying the biological functions and pathways related to the risk score. Figure 4a showed the top 10 BP terms, CC terms, and MF terms. The main enriched Go terms were digestion, neuronal cell body, synaptic membrane, passive transmembrane transporter activity, and channel activity. Figure 4b showed the top 6 KEGG pathways, namely the neuroactive ligand-receptor interaction, the ytokine-cytokine receptor interaction, the cAMP signaling pathway, the ECM-receptor interaction, and the PPAR signaling pathway.
The immune status was quantified by ssGSEA using the enrichment score and its association with the risk score was analyzed. Between two groups, the enrichment scores of Treg, TIL, Th1-cells, T-helpercells, pDCs, NK-cells, Neutrophils, Mast-cells, Macrophages, IDCs, and DCs were significantly different (adjusted P< 0.05, Figure 5a). The two groups  presented an obvious difference in terms of scores of MHC-class-I, Parainflammation, HLA, CCR, T-cellco-stimulation, T-cell-co-inhibition, Inflammationpromoting, APC-co-stimulation, APC-co-inhibition, Check-point, Type-I-IFN-Response, as well as Type-II-IFN-Response (adjusted P< 0.05, Figure 5b).

Internal validation of the 5-gene signature
We randomly selected half of the data to perform time-dependent ROC analysis for internal validation, the results showed that the area under the curve (AUC) of 0.613, 0.689, and 0.815 at 1 year, 2 years, 3 years, respectively ( Figure 6), which were close to our results (Figure 2f) .

Validation of the 5 ferroptosis-related genes in PTC tissues
We further validated the mRNA expression levels of the 5 ferroptosis-related genes in PTC tissues and normal tissues by RT-qPCR, the results showed that DPP4 (Figure 7a

Discussion
Ferroptosis was first proposed by Dixon [4] in 2012 and was crucial in the occurrence and development of many tumors. Previous studies [23,24] found several ferroptosis-related genes may be involved in PTC. However, their associations with the OS of PTC patients remained unknown.
In this research, we explored the correlation between 60 ferroptosis-related genes and the prognosis of PTC. We further constructed a new predictive model with 5 genes related to ferroptosis.  In our study, most ferroptosis-related genes (75%, 45/60) expressed differentially between tumor and adjacent normal tissues and six genes were related to the OS, which implied ferroptosis was involved in PTC and showed the possible predictive value of these genes.
We used 5 ferroptosis-related genes (DPP4, GSS, HMGCR, PGD, TFRC) to construct a predictive model. These genes act through different mechanisms [7,13]. DPP4, which is involved in glutathione and amino acid metabolism, promotes lipid peroxidation through NOXs. Ferroptosis could be limited by the p53 by blocking DPP4 activity in colorectal cancer cells [25]. GSS is the target gene of NRF2, which promotes resistance to ferroptosis [26]. The inhibition of HMGCR results in enhancing FIN-56-induced ferroptosis and blocking mevalonic acid synthesis [27]. PGD acts through pentose phosphate pathway and is an important regulator in tumor cells.
In non-small-cell lung cancer (NSCLC) cells Calu-1, knockdown of PGD inhibits erastin-triggered ferroptosis [27]. For iron metabolism, silencing TFRC can inhibit erastin-triggered ferroptosis [28]. Similarly, knockdown of TFRC inhibits ferroptosis caused by the deprivation of amino acid/cystine [29,30].To sum up, three of the genes (DPP4, PGD, TFRC) in the prognostic model can accelerate ferroptosis, while the other two genes (GSS, HMGCR) can protect cells from ferroptosis. Some studies [12,31] reported that DPP4 was upregulated in PTC patients and TFRC was associated with unfavorable prognoses. In our study, all the 5 ferroptosis-related genes were upregulated among PTC patients and were related to unfavorable clinical outcomes. However, whether these genes affecting the prognosis of PTC through regulating ferroptosis deserves further research.
As the mechanisms of ferroptosis in tumors was still elusive, we conducted GO enrichment and KEGG pathway analysis. We found digestion, neuronal cell body, synaptic membrane, passive transmembrane transporter activity, and channel activity were enriched, which provide a new direction for future research. Interestingly, the higherrisk group had remarkably higher contents of macrophages and Treg cells than group with low risk in this study. Considering tumor-associated macrophages [32,33] and Treg cells [34] are significantly related to cancer cell migration and poor prognoses in PTC patients, the unfavorable clinical outcomes in group with high risk might be caused by the impaired antitumor immunity.
Admittedly, our prognostic model was limited by using public databases. It is necessary to conduct more in vitro and in vivo studies for verifying the clinical application value owned by this model.

Conclusion
Our research identified a new predictive model using 5 ferroptosis-related genes. The risk score could predict the OS of PTC independently. The underlying mechanism of the association of the ferroptosis-related genes and PTC is unclear and deserves further study.