Abstract Background Numerous studies have implicated autophagy in the pathogenesis of thyroid carcinoma. This investigation aimed to establish an autophagy-related gene model and nomogram that can help predict the overall survival (OS) of patients with differentiated thyroid carcinoma (DTHCA). Methods Clinical characteristics and RNA-seq expression data from TCGA (The Cancer Genome Atlas) were used in the study. We also downloaded autophagy-related genes (ARGs) from the Gene Set Enrichment Analysis website and the Human Autophagy Database. First, we assigned patients into training and testing groups. R software was applied to identify differentially expressed ARGs for further construction of a protein-protein interaction (PPI) network for gene functional analyses. A risk score-based prognostic risk model was subsequently developed using univariate Cox regression and LASSO-penalized Cox regression analyses. The model’s performance was verified using Kaplan-Meier (KM) survival analysis and ROC curve. Finally, a nomogram was constructed for clinical application in evaluating the patients with DTHCA. Finally, a 7-gene prognostic risk model was developed based on gene set enrichment analysis. Results Overall, we identified 54 differentially expressed ARGs in patients with DTHCA. A new gene risk model based on 7-ARGs (CDKN2A, FGF7, CTSB, HAP1, DAPK2, DNAJB1, and ITPR1) was developed in the training group and validated in the testing group. The predictive accuracy of the model was reflected by the area under the ROC curve (AUC) values. Univariate and multivariate Cox regression analysis indicated that the model could independently predict the prognosis of patients with THCA. The constrained nomogram derived from the risk score and age also showed high prediction accuracy. Conclusions Here, we developed a 7-ARG prognostic risk model and nomogram for differentiated thyroid carcinoma patients that can guide clinical decisions and individualized therapy. Supplementary Information The online version contains supplementary material available at 10.1186/s12957-022-02590-6. Keywords: Differentiated thyroid carcinoma, Autophagy-related genes, Prognostic risk model, Nomogram, Gene set enrichment analysis, The Cancer Genome Atlas Introduction Autophagy involves autophagosome induction, phagophore nucleation and expansion, and lysosomal degradation [[37]1, [38]2]. During cellular stress, autophagy is activated to maintain cellular homeostasis and energy balance [[39]3]. Autophagy has been implicated in numerous cancer types [[40]4] but its role in cancer is controversial. For instance, inhibition of autophagy may contribute to tumour growth and metastasis but promote chemotherapeutic sensitivity [[41]5–[42]7]. Thyroid cancer (THCA) is the most common endocrine malignancy [[43]8]. Since the 1980s, advances in medical imaging technology and fine-needle aspiration technology for thyroid nodules have led to the detection of an increasing thyroid cancer incidence [[44]9–[45]11]. In 2017, the incidence of thyroid cancer in Korean women was 78.5/10^5, second only to breast cancer among female cancers [[46]12]. Based on pathological, clinical, and genetic features, malignant thyroid cancer can be classified into 5 subtypes: papillary, follicular, Hürthle cell, poorly differentiated, and anaplastic thyroid carcinoma [[47]13]. Follicular thyroid carcinoma and papillary thyroid carcinoma are known as differentiated thyroid carcinomas and account for > 96% of thyroid carcinomas; these subtypes have good overall survival after standardized treatment with surgery, radioiodine ablation, and hormone therapy [[48]14, [49]15]. However, distant metastasis, refractoriness to radioactive iodine, or locoregional recurrence significantly reduce overall DTHCA survival [[50]16–[51]18]. Numerous reports have implicated autophagy in thyroid carcinoma growth and treatment. For instance, inhibiting lactate dehydrogenase A (LDHA) activates autophagy via AMPK signalling, which has anti-cancer effects in papillary thyroid carcinoma [[52]19]. Recent studies have shown that baicalein induces autophagy through ERK/PI3K/Akt signalling to suppress the growth of thyroid cancer cells [[53]20]. Furthermore, autophagy can function as a therapeutic target for thyroid cancer treatment. For radioiodine (RAI)-refractory differentiated thyroid cancer, increasing the iodide uptake and concentration ability of thyroid follicular cells can significantly improve therapeutic efficacy [[54]21]. Some researchers have focused on mediating the active iodine uptake function of the sodium iodide symporter (NIS), finding that autophagy-activating compounds activate the accumulation of intracellular Ca^2+ and FOS, thus mediating the upregulation of NIS and restoration of RAI uptake [[55]21, [56]22]. Wang et al. also found that inhibition of autophagy increased the sensitivity of BRAF-mutated thyroid cancer cells to the chemotherapeutic drug vemurafenib [[57]23]. The above research results reflect the great potential of autophagic mechanism in the field of thyroid cancer treatment. The development of new autophagy-related diagnostic biomarkers may improve the early diagnosis and individualized treatment of DTHCA. Here, autophagy-related genes (ARGs) were retrieved from the HADb database and GSEA website while THCA patients’ clinical features and gene expression data were downloaded from The Cancer Genome Atlas (TCGA). After differential expression analysis (between tumour and normal tissue) and functional enrichment analysis, prognosis-related genes were identified through Cox regression analyses. Next, a prognostic model was developed using LASSO-penalized COX regression analysis and its accuracy and independence validated. Finally, a prognostic nomogram for estimating patients’ overall survival based on independent clinical features and risk scores was constructed. Materials and methods Acquisition of ARGs In total, 232 and 394 ARGs used in this study were acquired from the HADb database ([58]http://www.autophagy.lu) and GO_AUTOPHAGY gene datasets from the Gene Set Enrichment Analysis (GSEA, [59]http://www.gsea-msigdb.org/gsea/msigdb), respectively, as described previously [[60]24]. Selection of ARGs that were common to the 2 datasets resulted in the identification of 531 ARGs for further analysis. Clinical features and gene expression data of THCA patients in TCGA Clinical characteristics and transcriptome data were obtained from TCGA ([61]https://portal.gdc.cancer.gov). In particular, the RNA-seq dataset contained data from 510 DTHCA tumour and 58 non-tumour tissues. The DTHCA patients were randomly assigned to the training group (n = 252) and testing group (n = 249) by random selection in SPSS (v.22.0). Patient data included age, gender, survival status, and pathological stage. Selection of differentially expressed autophagy-related genes The Wilcoxon signed-rank test on R software (v.4.0.3) was employed to select differentially expressed autophagy-related genes (DEARGs). The cut-off criteria were | false discovery rate (FDR) < 0.05 and log[2](fold change) | ≥ 1. Gene function analysis To elucidate the DEARGs’ biological function, they were uploaded into Metascape ([62]https://metascape.org). Subsequently, enrichment analysis based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses was performed. The cut-off criteria were set as follows: p-value ≤ 0.05, Min Enrichment > 1.5, and Min Overlap = 3. Analysis of the PPI network A PPI network was designed on the STRING v.11.0 database ([63]https://www.string-db.org) and visualized with Cytoscape v.3.8.2 ([64]https://cytoscape.org). MCODE v2.0, a Cytoscape plugin was used to identify densely connected modules (MCODE score ≥ 4.0, degree cut-off = 2). Establishment of an ARG prognostic model The DEARG expression data were subjected to univariate Cox regression analysis to assess the associations of ARGs with OS (p ≤ 0.05 was applied as the inclusion criterion). RNA expression was normalized by log2 transformation. The hazard ratio (HR) was used to classify ARGs as protective or risk factors, with HR > 1 indicating protective factors and HR < 1 indicating risk factors. Based on LASSO-penalized Cox regression analysis, the prognosis-associated genes were used to develop a prognostic model. The risk score was the linear combination of gene expression and was calculated for each patient as follows: risk score= [MATH: ∑i=1nCoefi*Expi :MATH] . Coefi stands for the regression coefficient of an ARG. Expi indicates the relative expression values of that ARG. The patients were then classified into the high- and low-risk groups using the median risk score of the training group as cutoff. Thereafter, Kaplan-Meier (KM) analysis and two-sided log-rank test were applied to compare OS between the two groups. The R package “timeROC” was used for ROC curve analysis to evaluate the prognostic model’s specificity and sensitivity. The ability of the prognostic risk model and clinical characteristics to independently predict OS was elucidated based on univariate and multivariate Cox regression analysis. Subgroup analysis was used to evaluate the model’s applicability. Multivariate ROC analysis To compare the prognostic accuracy of our ARG signature to the TNM staging system, multivariate ROC analysis was performed by R package “timeROC”, “survival”, and “survminer”. Construction and validation of the nomogram To enable clinicians to conveniently use the prognostic model in assessing the 3- and 5-year OS of patients with DTHCA, we used the independent risk factors (risk score and age) to create a nomogram and then validated its predictive ability by plotting calibration curves for comparing differences between predicted and actual survival. The Harrell concordance index (C-index) was used to assess the nomogram’s predictive accuracy. Gene set enrichment analysis Gene set enrichment analysis (GSEA) for both the high- and low-risk groups was performed based on the C2.CP.KEGG.v7.4 gene set using GSEA software (version 4.1.0). The statistical significance criteria were set as NOM p < 0.05 and FDR < 0.25. Statistical analysis Data analysis was performed using R version 4.0.3. KM analysis was performed using the R packages “survival” and “survminer”, and differences were assessed by the log-rank test. ROC curves were plotted using the R package “timeROC”. Univariate and multivariate Cox regression analysis was used to test the model’s prognostic independence. p ≤ 0.05 indicates statistical significance. Cell culture The human normal thyroid epithelial cell line N-thy-ori-3-1 and the thyroid cancer cell lines FTC-133 and TPC-1 were obtained from Procell (Wuhan, China). All cell lines were cultured in DMEM (Irvine Scientific, Carlsbad, CA) supplemented with 10% foetal bovine serum (FBS). RNA extraction and qPCR analysis RNA was extracted from the cell using HiPure Total RNA Mini Kit (R4111-03, Magen, China). HiScript II QRT SuperMix Kit (Vazyme, China) was used for reverse transcription. qRT-PCR was conducted using the SYBR GREEN MIX Kit (Vazyme, China) by the CFX96 Real-time PCR Detection System (Bio-Rad, USA). GAPDH was selected as the internal housekeeping gene, and relative gene expression was calculated by the 2^−ΔΔCt method. Each qRT-PCR was repeated three times. Results Identification of DEARGs in THCA patients A total of 531 identified ARGs were used in the study (Table S[65]1). Gene expression data and clinical information on 58 normal and 510 DTHCA tissues were downloaded from the TCGA database (Table [66]1). Application of the criteria FDR < 0.05 and | log[2](FC) | ≥ 1 resulted in 54 DEARGs. Of these, 23 were downregulated and 31 were upregulated (Fig. [67]1A–C) Table 1. Clinical characteristics of the TCGA-THCA dataset Clinical characteristic TCGA-THCA (n = 501) Percentage Survival status  Alive 485 96.8%  Dead 16 3.2% Age  ≥ 65 76 15.2%  < 65 425 84.8% Gender  Female 366 73.1%  Male 135 26.9% T  T1 142 28.3%  T2 164 32.7%  T3 170 33.9%  T4 23 4.6%  TX 2 0.5% N  N0 229 45.7%  N1 222 44.3%  NX 50 10.0% M  M0 282 55.3%  M1 9 1.8%  MX 210 42.9% AJCC stage  Stage I 281 56.1%  Stage II 52 10.4%  Stage III 111 22.2%  Stage IV 55 10.9%  Unknown 2 0.4% [68]Open in a new tab Fig. 1. [69]Fig. 1 [70]Open in a new tab Differentially expressed autophagy related genes (DEARGs) between thyroid carcinoma tissues and non-tumour normal tissues. (A)Volcano plot for the DEARGs blue: downregulate red: upregulate. (B) Heatmap of 54 DEARGs. The depth of blue indicates the level of low expression, and the depth of red indicates the level of high expression. (C) Boxplot of the DEARGs Functional annotation of DEARGs To elucidate the biological function and mechanisms of the 54 DEARGs, GO biological process and KEGG pathway enrichment analyses were done. The DEARGs showed strong enrichment in the GO terms autophagy, intrinsic apoptotic pathway in response to endoplasmic reticulum stress and positive regulation of cell death (Fig. [71]2A, B). KEGG pathway enrichment analysis revealed that the DEARGs were closely related to pathway autophagy-animal and apoptosis (Fig. [72]2C, D). Fig. 2. [73]Fig. 2 [74]Open in a new tab GO and KEGG functional enrichment analysis. (A) The circus plot of GO biological processes enrichment analysis. (B) Heatmap for the GO biological processes enrichment analysis. (C) The bubble plot of KEGG pathway enrichment analysis. (D) Heatmap for the KEGG pathway enrichment analysis Analysis of the PPI network Examination of the PPI network of interaction between the DEARGs revealed 47 nodes and 208 edges (Fig. [75]3A; red = upregulation, green = downregulation). Next, 30 genes with > 2 edges were selected as the hub DEARGs for further analysis (Fig. [76]3B). Enrichment analysis of the module with the highest score revealed that the DEARGs in the module were associated with autophagy (Fig. [77]3C, D). Fig. 3. [78]Fig. 3 [79]Open in a new tab Protein-Protein Interaction (PPI) network analysis of DEARGs in THCA: (A) PPI network plot for DEAGRs Red: upregulate genes Green: downregulate genes. (B) Hub genes with more than 2 interactions. (C) Subnet modules of PPI network analysis. (D) GO biological processes enrichment analysis for the subnet modules Design of the ARG-based prognostic model Univariate Cox regression analysis of the correlations of ARGs with OS in THCA patients identified 10 prognosis-related genes (p < 0.05, Fig. [80]4A). LASSO-penalized Cox regression analysis was performed to establish a risk score prognostic model that included 7-genes (Fig. [81]4B, C). The risk score for each patient was determined using the equation: Risk score = (0.46026029*CDKN2A expression) − (0.29526375*CTSB expression) + (0.36900263*FGF7 expression) + (0.36720186*HAP1 expression) − (0.03896433*DAPK2 expression) + (0.20590309*DNAJB1 expression) + (0.46793624*ITPR1 expression) and median risk score was set as the cut-off value point to divide the training group into a high-risk group (n = 126) and low-risk group (n = 126). Results of the Kaplan-Meier analysis confirmed that patients with high-risk scores had poorer prognosis (p = 0.002, Fig. [82]5A). ROC analysis of the training group revealed AUC values for 1-, 3-, and 5-year OS of 0.765, 0.773, and 0.897, respectively (Fig. [83]5C). Additionally, THCA patients in the training group were ranked by risk score to reveal the relationship between risk score and prognosis. This analysis showed that most of the dead patients were from the high-risk group and that higher risk scores correlated with increasing mortality rate (Fig. [84]5E, G, I). To validate the prognostic model, similar analyses were done on the testing group and testing group on the basis of the risk score and the same cut-off criteria used for analysis in the training group (Fig. [85]5F, H, J). Consistent with the training group, patients with high-risk scores showed poorer prognosis than those with low scores (p = 0.013, Fig. [86]5B). The ROC curve AUC values of the testing group were 0.967, 0.97, and 0.751 for 1-, 3-, and 5-year survival (Fig. [87]5D). Together, these data show that the prognostic model has good predictive accuracy. Fig. 4. [88]Fig. 4 [89]Open in a new tab Key prognostic-related genes and construction of prognostic risk model: (A) Forest plot of 10 prognostic-related genes (B) Lasso Cox regression of 7 genes used in the prognostic risk model. (C) Lasso filters variables Fig. 5. [90]Fig. 5 [91]Open in a new tab Performance and validation of the prognostic risk model: (A) Kaplan-Meier analysis of the training group. (B) Kaplan-Meier analysis of the testing group. (C) The ROC curve of overall survival for training group. (D) The ROC curve of overall survival for testing group. (E) The number of different risk group patients in training group. (F) The number of different risk group patients in testing group. (G) Survival status in training group. (H) Survival status in testing group. (I) Expression level of risk prognostic model gene in the training group. (J) Expression level of risk prognostic model gene in the testing group Independent prognostic testing of the prognostic model Using univariate and multivariate Cox regression analysis, we assessed whether clinical features such as age, gender, tumour stage, N (lymph node), T (tumour), M (metastasis) (Table S[92]2), as well as the risk score prognostic model independently correlated with OS. Univariate regression analysis found that age, tumour stage, T, and risk score influenced the OS. Multivariate regression analysis revealed that age (p = 0.01, HR = 1.134 95%, CI=1.031–1.248) and risk score (p = 0.001, HR = 2.143 95% CI = 1.346–3.413) were independent risk factors for THCA prognosis (Fig. [93]6A, B). Clinical subgroup analysis with the prognostic model showed that prognosis was poorer in the high-risk group for both female (p = 0.006) and male (p = 0.003), disease stage III–IV (p < 0.001), age < 65 group (p = 0.042), and age ≥ 65 group (p = 0.035) (Fig. [94]6C–E). These data indicate that the model can be applied to predict prognosis independently. Fig. 6. [95]Fig. 6 [96]Open in a new tab Independent prognostic testing and clinical subgroup analysis of the prognostic model: (A, B) Univariate and multivariate Cox regression analysis of clinical characteristic and risk score. (C) Gender. (D) Tumour stage. (E) Age Comparison between ARG signature and TNM staging system We plotted the ROC curves to compare the ARG signature and TNM staging system in the prediction of thyroid cancer prognosis. The ROC curves’ AUC values of ARG signature for 1-, 3-, and 5-year OS were 0.823 (Fig. [97]7A), 0.871 (Fig. [98]7B), and 0.912 (Fig. [99]7C) higher than the AUC values of the TNM staging system. Multivariate ROC analysis revealed that our ARG signature displayed high prognostic value compared with the TNM staging system. Fig. 7. [100]Fig. 7 [101]Open in a new tab Multivariate ROC analysis compare the predictive ability of ARG signature to TNM staging system: (A) Multivariate ROC curves for 1-year OS. (B) Multivariate ROC curves for 3-year OS. (C) Multivariate ROC curves for 5-year OS Construction and validation of the nomogram To make the prognostic model convenient for clinicians, we quantified it and included other independent prognostic factors to develop a nomogram. The nomogram for 3- and 5-year OS is shown in Fig. [102]8A. The nomogram’s concordance index (C-index) of 0.911 (se = 0.044) showed that the nomogram’s predictive ability was moderate (Fig. [103]8B, C). Fig. 8. [104]Fig. 8 [105]Open in a new tab Nomogram and calibration curve of prognostic risk model: (A) Nomogram to predict 3- and 5-year overall survival of thyroid cancer (THCA) patients. (B, C) The calibration curves for 3-year (B) and 5-year (C) survival Pathways associated with the prognostic genes In the high-risk group, the genes were linked to autophagy and cancer pathways, including MAPK, MTOR, PPAR, regulation of autophagy, and WNT signalling (Fig. [106]9). Fig. 9. [107]Fig. 9 [108]Open in a new tab Gene set enrichment analysis Expression level of signature mRNAs in the thyroid cancer cell We performed qRT-PCR to validate the signature mRNA expression level. The results showed that ITPR1 was highly expressed in both FTC-133 and TPC-1 thyroid cancer cell lines, and the hazard ratio of ITPR1 indicated that ITPR1 was an oncogene, and showed research potential of its biological mechanism. The expression of some mRNAs (FGF7, DAPK2) remains unknown due to missing primer sequence (Fig. [109]10). Fig. 10. Fig. 10 [110]Open in a new tab RT-qPCR analysis for the signature mRNAs. (ns: not significant, *p < 0.05, **:p < 0.01, ***p < 0.001, ****p < 0.0001) Discussion Rapid advances in sequencing technology have led to the development of molecular markers to facilitate the detection and monitoring of the clinical course of many cancers [[111]25]. However, the use of individual genetic markers in determining thyroid cancer prognosis has significant limitations. For example, RAS mutations, especially NRAS codon 61 mutations are associated with high incidence of distant metastasis but do not affect overall survival in THCA patients [[112]26]. A study involving 599 PTC patients found that although the BRAF V600E mutation significantly increases PTC recurrence, it is not suitable as an independent prognostic factor [[113]27]. Thus, an effective, gene-based independent prognostic model is needed to guide personalized therapy. Numerous studies have implicated that autophagy is crucial for the tumour progression and efficacy of targeted therapy in thyroid carcinoma. Thus, we assume autophagy-related genes have potential prognostic value in DTHCA [[114]28–[115]30]. To confirm our hypothesis, a multiple-step bioinformatic analysis was performed. First, we acquired autophagy-related genes from the HADb database and GO_AUTOPHAGY gene datasets. Differential expression analysis and PPI network analysis were performed to identify the hub autophagy-related genes in DTHCA. Next, we used hub autophagy-related genes to screen out the autophagy-related genes associated with prognosis by univariate Cox regression analysis and constructed a seven-autophagy-related gene signature to predict DTHCA patients’ prognosis through LASSO-penalized Cox regression analysis. Here, we constructed a prognostic model for DTHCA comprising CDKN2A, FGF7, CTSB, HAP1, DAPK2, DNAJB1, and ITPR1 genes. CDKN2A, also known as multiple tumour suppressor 1, is an inhibitor of CDK4 that suppresses cancer cell proliferation by arresting the cell cycle at phase G1 [[116]31–[117]33]. However, CDKN2A upregulation may enhance autophagy in non-endometrioid endometrial carcinoma, which protects cells against unfavourable conditions, promoting tumorigenesis [[118]34]. A previous study revealed that CDKN2A CpG islands’ methylation results in the inactivation of the CDKN2A gene, the phenomenon of hyper-methylation was more present in the PTC patients with metastases and high AMES risk [[119]35]. This research result indicated CDKN2A might function as a tumour suppressor in PTC, but the autophagy mediated by CDKN2A effect on the thyroid has not been reported until now. FGF7, a member of the FGF family, exerts biological functions by interacting with its high-affinity binding receptor, FGFR2-IIIb [[120]36]. The binding of FGF7 and FGFR2-IIIb could signal to SRC tyrosine kinase inducing phosphorylation of F-actin binding protein and cortactin, which promote the migration of tumour cells [[121]37]. In our bioinformatics analysis, FGF7 is lowly expressed in thyroid cancer; Tetsuo et al. demonstrated that this expression is associated with its DNA promoter methylation [[122]38]. It has been reported that in thyroid epithelial cells, FGFR2-IIIb reduces the expression of FGF7 and increases the ratio of FGF4/FGF7 and couples epithelial signalling with expansion of stomal to favour tumour growth [[123]39]. DNAJB1, which belongs to the heat shock 40 protein family is associated with autophagy and participates in tumour progression [[124]40]. In autophagy, DNAJB1 acts as a link between WIPI2 and ATG2A and WIPI2 depletion reduces the autophagosome number [[125]41]. In addition, DNAJB1 is correlated to the essential tumour suppressor gene p53. Cui et al. found that DNAJB1 interacts with PDCD5 and can facilitate tumour progression via inhibiting the apoptotic function of wtp53 [[126]42]. DAPK2 is located on chromosome 12p11.21 and encodes a calcium/calmodulin (Ca^2+/CaM) dependent protein kinase [[127]43]. Beclin-1, a core autophagic protein is a direct DAPK2 substrate. Activated DAPK2 induces autophagy via Beclin-1 phosphorylation [[128]44]. 2,3′,4,4′,5-pentachlorobiphenyl (PCB118) enhances autophagy and damages thyroid ultrastructure via the DAPK2/PKD/VPS34 pathway [[129]45]. Jiang et al. revealed that DAPK2 could promote I-κBα degradation through inducing selective autophagy and further contributing to anaplastic thyroid carcinoma development [[130]46]. The effect of DAPK2 on differentiated thyroid tumours has remained to be reported. HAP1, encoding the Huntington (HTT) gene, is expressed in the nervous system and endocrine organs, such as the thyroid gland [[131]47]. Several studies have shown that HAP1 function as an autophagy-regulating gene and plays a role in the development of neurodegenerative diseases and neurodevelopmental disorders [[132]48, [133]49]. Diagnostic and therapeutic potential has been reported for HAP1 in breast and pancreatic cancer, but its role in thyroid cancer needs further study [[134]50, [135]51]. In mammals, ITPR1 participate in encoding IP3-gated calcium channel, regulating calcium homeostasis in endoplasmic reticulum [[136]52]. ITPR1-mediated induction of autophagy is thought to promote renal cell carcinoma by suppressing NK cells [[137]53]. In papillary thyroid carcinoma, expression of ITPR1 was promoted via effect of lncRNA SLC26A4-AS1 mediated ETS1 recruitment, suppressing tumour growth by enhancing autophagy [[138]54]. However, regression analysis revealed high ITPR1 expression is associated with poor THCA prognosis. CTSB belongs to the cysteine protease family that affects tumour growth and invasion through interaction with cellular proteins [[139]55]. CTSB upregulation in papillary thyroid cancer correlates with positive lymph node metastasis and cancer cell migration via p38-mediated EMT [[140]56]. CTSB regulates autophagy via enhancing the lysosomal membrane permeabilization process and reducing the functional lysosomes required by autophagy flux [[141]57]. High CTSB expression has been associated with suberoylanilide hydroxamic acid (SAHA)-induced autophagy, which suppresses breast cancer growth [[142]58]. This finding is consistent with our past findings and shows that the role of CTSB in thyroid cancer merits further investigation. Several analyses suggested that our prognostic model has excellent accuracy and independence. Our prognostic model’s predictive ability for DTHCA patient’s OS is better than the widely used TNM staging system by conducting multivariate ROC analysis. However, this result still needs more patients to validate in the clinical practice. After that, a nomogram was constructed to screen out the relatively high-risk patients in clinical practice; and this population might need extended resection and lymphadenectomy. Nomograms are widely used to predict prognosis, skip metastasis, and central lymph node metastasis in thyroid carcinoma [[143]59–[144]61]. To give clinicians a practical tool for individualized prognosis prediction, we constructed a nomogram based on the independent factors in the multivariate Cox regression analysis. It is worth mentioning that our nomogram has better predictive ability than other thyroid prognosis-related nomograms, as quantified by the concordance index (C-index:0.911 95% CI:0.867–0.955), the C-index of Wang’s nomogram was 0.797 (95% CI: 0.730–0.864), and that of Zhang’s nomogram was 0.717 (95% CI, 0.603–0.831) [[145]62, [146]63]. However, this study had two main limitations. First, all transcriptome data and information on clinical characteristics were obtained from public datasets and lacked external validation. Second, further clinical cases and experiments are needed to clarify the mechanism by which ARGs affect thyroid cancer development. Despite these limitations, our prognostic model showed excellent predictive accuracy and clinical applicability and may guide individualized therapy. Conclusion In conclusion, we developed a novel seven-mRNA (CDKN2A, FGF7, CTSB, HAP1, DAPK2, DNAJB1, ITPR1) risk model and nomogram that could help assess the prognosis of patients and guide clinical decision-making and individualized therapy for differentiated thyroid cancer. Supplementary Information [147]Additional file 1. ^(3.9KB, txt) [148]Additional file 2. ^(456B, txt) [149]Additional file 3. ^(622B, txt) [150]Additional file 4. ^(473B, txt) Acknowledgements