Skip to main content

A seven-autophagy-related gene signature for predicting the prognosis of differentiated thyroid carcinoma

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.

Introduction

Autophagy involves autophagosome induction, phagophore nucleation and expansion, and lysosomal degradation [1, 2]. During cellular stress, autophagy is activated to maintain cellular homeostasis and energy balance [3]. Autophagy has been implicated in numerous cancer types [4] but its role in cancer is controversial. For instance, inhibition of autophagy may contribute to tumour growth and metastasis but promote chemotherapeutic sensitivity [5,6,7].

Thyroid cancer (THCA) is the most common endocrine malignancy [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 [9,10,11]. In 2017, the incidence of thyroid cancer in Korean women was 78.5/105, second only to breast cancer among female cancers [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 [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 [14, 15]. However, distant metastasis, refractoriness to radioactive iodine, or locoregional recurrence significantly reduce overall DTHCA survival [16,17,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 [19]. Recent studies have shown that baicalein induces autophagy through ERK/PI3K/Akt signalling to suppress the growth of thyroid cancer cells [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 [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 Ca2+ and FOS, thus mediating the upregulation of NIS and restoration of RAI uptake [21, 22]. Wang et al. also found that inhibition of autophagy increased the sensitivity of BRAF-mutated thyroid cancer cells to the chemotherapeutic drug vemurafenib [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 (http://www.autophagy.lu) and GO_AUTOPHAGY gene datasets from the Gene Set Enrichment Analysis (GSEA, http://www.gsea-msigdb.org/gsea/msigdb), respectively, as described previously [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 (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 log2(fold change) | ≥ 1.

Gene function analysis

To elucidate the DEARGs’ biological function, they were uploaded into Metascape (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 (https://www.string-db.org) and visualized with Cytoscape v.3.8.2 (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=\({\sum}_{i=1}^n\left( Coefi\ast Expi\right)\). 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 S1). Gene expression data and clinical information on 58 normal and 510 DTHCA tissues were downloaded from the TCGA database (Table 1). Application of the criteria FDR < 0.05 and | log2(FC) | ≥ 1 resulted in 54 DEARGs. Of these, 23 were downregulated and 31 were upregulated (Fig. 1A–C)

Table 1 Clinical characteristics of the TCGA-THCA dataset
Fig. 1
figure 1

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. 2A, B). KEGG pathway enrichment analysis revealed that the DEARGs were closely related to pathway autophagy-animal and apoptosis (Fig. 2C, D).

Fig. 2
figure 2

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. 3A; red = upregulation, green = downregulation). Next, 30 genes with > 2 edges were selected as the hub DEARGs for further analysis (Fig. 3B). Enrichment analysis of the module with the highest score revealed that the DEARGs in the module were associated with autophagy (Fig. 3C, D).

Fig. 3
figure 3

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. 4A). LASSO-penalized Cox regression analysis was performed to establish a risk score prognostic model that included 7-genes (Fig. 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. 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. 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. 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. 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. 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. 5D). Together, these data show that the prognostic model has good predictive accuracy.

Fig. 4
figure 4

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
figure 5

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 S2), 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. 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. 6C–E). These data indicate that the model can be applied to predict prognosis independently.

Fig. 6
figure 6

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. 7A), 0.871 (Fig. 7B), and 0.912 (Fig. 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
figure 7

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. 8A. The nomogram’s concordance index (C-index) of 0.911 (se = 0.044) showed that the nomogram’s predictive ability was moderate (Fig. 8B, C).

Fig. 8
figure 8

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. 9).

Fig. 9
figure 9

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. 10).

Fig. 10
figure 10

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 [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 [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 [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 [28,29,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 [31,32,33]. However, CDKN2A upregulation may enhance autophagy in non-endometrioid endometrial carcinoma, which protects cells against unfavourable conditions, promoting tumorigenesis [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 [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 [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 [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 [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 [39]. DNAJB1, which belongs to the heat shock 40 protein family is associated with autophagy and participates in tumour progression [40]. In autophagy, DNAJB1 acts as a link between WIPI2 and ATG2A and WIPI2 depletion reduces the autophagosome number [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 [42]. DAPK2 is located on chromosome 12p11.21 and encodes a calcium/calmodulin (Ca2+/CaM) dependent protein kinase [43]. Beclin-1, a core autophagic protein is a direct DAPK2 substrate. Activated DAPK2 induces autophagy via Beclin-1 phosphorylation [44]. 2,3′,4,4′,5-pentachlorobiphenyl (PCB118) enhances autophagy and damages thyroid ultrastructure via the DAPK2/PKD/VPS34 pathway [45]. Jiang et al. revealed that DAPK2 could promote I-κBα degradation through inducing selective autophagy and further contributing to anaplastic thyroid carcinoma development [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 [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 [48, 49]. Diagnostic and therapeutic potential has been reported for HAP1 in breast and pancreatic cancer, but its role in thyroid cancer needs further study [50, 51]. In mammals, ITPR1 participate in encoding IP3-gated calcium channel, regulating calcium homeostasis in endoplasmic reticulum [52]. ITPR1-mediated induction of autophagy is thought to promote renal cell carcinoma by suppressing NK cells [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 [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 [55]. CTSB upregulation in papillary thyroid cancer correlates with positive lymph node metastasis and cancer cell migration via p38-mediated EMT [56]. CTSB regulates autophagy via enhancing the lysosomal membrane permeabilization process and reducing the functional lysosomes required by autophagy flux [57]. High CTSB expression has been associated with suberoylanilide hydroxamic acid (SAHA)-induced autophagy, which suppresses breast cancer growth [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 [59,60,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) [62, 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.

Availability of data and materials

The data that support the findings of this study are available in The Cancer Genome Atlas databases at https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga, reference number: TCGA-THCA.

Abbreviations

OS:

Overall survival

THCA:

Thyroid cancer

DTHCA:

Differentiated thyroid carcinoma

AUC:

Area under ROC curve

ARGs:

Autophagy-related genes

TCGA:

The Cancer Genome Atlas

PPI:

Protein-protein interaction

FDR:

False discovery rate

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

HR:

Hazard ratio

GSEA:

Gene set enrichment analysis

RAI:

Radioiodine

NIS:

Sodium iodide symporter

References

  1. Sun T, Li X, Zhang P, Chen WD, Zhang HL, Li DD, et al. Acetylation of Beclin 1 inhibits autophagosome maturation and promotes tumor growth. Nat. Commun. 2015;6:7215 https://doi.org/10.1038/ncomms8215.

    PubMed  Article  Google Scholar 

  2. Xie Y, Kang R, Sun X, Zhong M, Huang J, Klionsky DJ, et al. Posttranslational modification of autophagy-related proteins in macroautophagy. Autophagy. 2015;11(1):28–45 https://doi.org/10.4161/15548627.2014.984267.

    PubMed  Article  CAS  Google Scholar 

  3. Feng Y, He D, Yao Z, Klionsky DJ. The machinery of macroautophagy. Cell Res. 2014;24(1):24–41 https://doi.org/10.1038/cr.2013.168.

    CAS  PubMed  Article  Google Scholar 

  4. Onorati AV, Dyczynski M, Ojha R, Amaravadi RK. Targeting autophagy in cancer. Cancer. 2018;124(16):3307–18 https://doi.org/10.1002/cncr.31335.

    PubMed  Article  Google Scholar 

  5. Nassour J, Radford R, Correia A, Fusté JM, Schoell B, Jauch A, et al. Autophagic cell death restricts chromosomal instability during replicative crisis. Nature. 2019;565(7741):659–63 https://doi.org/10.1038/s41586-019-0885-0.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. Dower CM, Wills CA, Frisch SM, Wang HG. Mechanisms and context underlying the role of autophagy in cancer metastasis. Autophagy. 2018;14(7):1110–28 https://doi.org/10.1080/15548627.2018.1450020.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. Pecoraro A, Carotenuto P, Franco B, De Cegli R, Russo G, Russo A. Role of uL3 in the crosstalk between nucleolar stress and autophagy in colon cancer cells. Int. J. Mol. Sci. 2020;21(6):2143 https://doi.org/10.3390/ijms21062143.

    CAS  PubMed Central  Article  Google Scholar 

  8. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. Cancer J Clin. 2021;71(1):7–33 https://doi.org/10.3322/caac.21654.

    Article  Google Scholar 

  9. Kitahara CM, Sosa JA. Understanding the ever-changing incidence of thyroid cancer. Nature reviews. Endocrinol. 2020;16(11):617–8 https://doi.org/10.1038/s41574-020-00414-9.

    Google Scholar 

  10. Lim H, Devesa SS, Sosa JA, Check D, Kitahara CM. Trends in thyroid cancer incidence and mortality in the United States, 1974-2013. JAMA. 2017;317(13):1338–48 https://doi.org/10.1001/jama.2017.2719.

    PubMed  PubMed Central  Article  Google Scholar 

  11. Li M, Dal Maso L, Vaccarella S. Global trends in thyroid cancer incidence and the impact of overdiagnosis. Lancet Diabetes Endocrinol. 2020;8(6):468–70 https://doi.org/10.1016/S2213-8587(20)30115-7.

    PubMed  Article  Google Scholar 

  12. Hong S, Won YJ, Park YR, Jung KW, Kong HJ, Lee ES, et al. Cancer statistics in Korea: incidence, mortality, survival, and prevalence in 2017. Cancer Res. Treat. 2020;52(2):335–50 https://doi.org/10.4143/crt.2020.206.

    PubMed  PubMed Central  Article  Google Scholar 

  13. Asa SL. The current histologic classification of thyroid cancer. Endocrinol. Metab. Clin. North Am. 2019;48(1):1–22 https://doi.org/10.1016/j.ecl.2018.10.001.

    PubMed  Article  Google Scholar 

  14. Yan KL, Li S, Tseng CH, Kim J, Nguyen DT, Dawood NB, Livhits MJ, Yeh MW, Leung AM. Rising incidence and incidence-based mortality of thyroid cancer in California, 2000-2017. J Clin Endocrinol Metab. 2020;105(6):dgaa121. https://doi.org/10.1210/clinem/dgaa121

  15. Ibrahim EY, Busaidy NL. Treatment and surveillance of advanced, metastatic iodine-resistant differentiated thyroid cancer. Curr. Opin. Oncol. 2017;29(2):151–8 https://doi.org/10.1097/CCO.0000000000000349.

    PubMed  Article  Google Scholar 

  16. Kreissl MC, Janssen M, Nagarajah J. Current treatment strategies in metastasized differentiated thyroid cancer. J Nuclear Med. 2019;60(1):9–15 https://doi.org/10.2967/jnumed.117.190819.

    CAS  Article  Google Scholar 

  17. Jin Y, Van Nostrand D, Cheng L, Liu M, Chen L. Radioiodine refractory differentiated thyroid cancer. Crit. Rev. Oncol. Hematol. 2018;125:111–20 https://doi.org/10.1016/j.critrevonc.2018.03.012.

    PubMed  Article  Google Scholar 

  18. Chereau N, Oyekunle TO, Zambeli-Ljepović A, Kazaure HS, Roman SA, Menegaux F, et al. Predicting recurrence of papillary thyroid cancer using the eighth edition of the AJCC/UICC staging system. Br. J. Surg. 2019;106(7):889–97 https://doi.org/10.1002/bjs.11145.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  19. Hou X, Shi X, Zhang W, Li D, Hu L, Yang J, et al. LDHA induces EMT gene transcription and regulates autophagy to promote the metastasis and tumorigenesis of papillary thyroid carcinoma. Cell Death Dis. 2021;12(4):347 https://doi.org/10.1038/s41419-021-03641-8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. Wang M, Qiu S, Qin J. Baicalein induced apoptosis and autophagy of undifferentiated thyroid cancer cells by the ERK/PI3K/Akt pathway. Am. J. Transl. Res. 2019;11(6):3341–52.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. Oh JM, Ahn BC. Molecular mechanisms of radioactive iodine refractoriness in differentiated thyroid cancer: impaired sodium iodide symporter (NIS) expression owing to altered signaling pathway activity and intracellular localization of NIS. Theranostics. 2021;11(13):6251–77 https://doi.org/10.7150/thno.57689.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. Tesselaar MH, Crezee T, Swarts HG, Gerrits D, Boerman OC, Koenderink JB, et al. Digitalis-like compounds facilitate non-medullary thyroid cancer redifferentiation through intracellular Ca2+, FOS, and autophagy-dependent pathways. Mol. Cancer Ther. 2017;16(1):169–81 https://doi.org/10.1158/1535-7163.MCT-16-0460.

    CAS  PubMed  Article  Google Scholar 

  23. Wang W, Kang H, Zhao Y, Min I, Wyrwas B, Moore M, et al. Targeting autophagy sensitizes BRAF-mutant thyroid cancer to vemurafenib. J. Clin. Endocrinol. Metab. 2017;102(2):634–43 https://doi.org/10.1210/jc.2016-1999.

    PubMed  Article  Google Scholar 

  24. Wang Y, Zhao W, Xiao Z, Guan G, Liu X, Zhuang M. A risk signature with four autophagy-related genes for predicting survival of glioblastoma multiforme. J. Cell. Mol. Med. 2020;24(7):3807–21 https://doi.org/10.1111/jcmm.14938.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  25. Jayanthi V, Das AB, Saxena U. Recent advances in biosensor development for the detection of cancer biomarkers. Biosens Bioelectron. 2017;91:15–23 https://doi.org/10.1016/j.bios.2016.12.014.

    CAS  PubMed  Article  Google Scholar 

  26. Jang EK, Song DE, Sim SY, Kwon H, Choi YM, Jeon MJ, Han JM, Kim WG, Kim TY, Shong YK, Kim WB. NRAS codon 61 mutation is associated with distant metastasis in patients with follicular thyroid carcinoma. Thyroid. 2014;24(8):1275–81. https://doi.org/10.1089/thy.2014.0053.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. Enumah S, Fingeret A, Parangi S, Dias-Santagata D, Sadow PM, Lubitz CC. BRAFV600E mutation is associated with an increased risk of papillary thyroid cancer recurrence. World J Surg. 2020;44(8):2685–91 https://doi.org/10.1007/s00268-020-05521-2.

    PubMed  PubMed Central  Article  Google Scholar 

  28. Zhou H, Xie X, Chen Y, Lin Y, Cai Z, Ding L, Wu Y, Peng Y, Tang S, Xu H. Chaperone-mediated autophagy governs progression of papillary thyroid carcinoma via PPARγ-SDF1/CXCR4 signaling. J Clin Endocrinol Metabol. 2020;105(10):dgaa366. https://doi.org/10.1210/clinem/dgaa366

  29. Yi H, Ye T, Ge M, Yang M, Zhang L, Jin S, et al. Inhibition of autophagy enhances the targeted therapeutic effect of sorafenib in thyroid cancer. Oncol Rep. 2018;39(2):711–20 https://doi.org/10.3892/or.2017.6118.

    CAS  PubMed  Google Scholar 

  30. Liu K, Yu Q, Li H, Xie C, Wu Y, Ma D, et al. BIRC7 promotes epithelial-mesenchymal transition and metastasis in papillary thyroid carcinoma through restraining autophagy. Am J Cancer Res. 2020;10(1):78–94.

    PubMed  PubMed Central  Google Scholar 

  31. Mao X, Li B, Liang Y, Li S, Zhou J, He Q, et al. Auxiliary diagnostic value of p16 amplification combined with the detection of heterozygous and homozygous loss for urothelial carcinoma. Oncol Lett. 2018;15(5):6533–40 https://doi.org/10.3892/ol.2018.8137.

    PubMed  PubMed Central  Google Scholar 

  32. Serrano M. The tumor suppressor protein p16INK4a. Experiment Cell Rese. 1997;237(1):7–13 https://doi.org/10.1006/excr.1997.3824.

    CAS  Article  Google Scholar 

  33. Appay R, Dehais C, Maurage CA, Alentorn A, Carpentier C, Colin C, et al. CDKN2A homozygous deletion is a strong adverse prognosis factor in diffuse malignant IDH-mutant gliomas. Neuro-Oncol. 2019;21(12):1519–28 https://doi.org/10.1093/neuonc/noz124.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Devis-Jauregui L, Eritja N, Davis ML, Matias-Guiu X, Llobet-Navàs D. Autophagy in the physiological endometrium and cancer. Autophagy. 2021;17(5):1077–95 https://doi.org/10.1080/15548627.2020.1752548.

    CAS  PubMed  Article  Google Scholar 

  35. Wang P, Pei R, Lu Z, Rao X, Liu B. Methylation of p16 CpG islands correlated with metastasis and aggressiveness in papillary thyroid carcinoma. J Chin Med Assoc. 2013;76(3):135–9. https://doi.org/10.1016/j.jcma.2012.11.007.

    CAS  Article  PubMed  Google Scholar 

  36. Lezmi G, Verkarre V, Khen-Dunlop N, Vibhushan S, Hadchouel A, Rambaud C, et al. FGF10 Signaling differences between type I pleuropulmonary blastoma and congenital cystic adenomatoid malformation. Orphanet J Rare Dis. 2013;8:130 https://doi.org/10.1186/1750-1172-8-130.

    PubMed  PubMed Central  Article  Google Scholar 

  37. Belleudi F, Scrofani C, Torrisi MR, Mancini P. Polarized endocytosis of the keratinocyte growth factor receptor in migrating cells: role of SRC-signaling and cortactin. PloS one. 2011;6(12):e29159 https://doi.org/10.1371/journal.pone.0029159.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  38. Kondo T, Zhu X, Asa SL, Ezzat S. The cancer/testis antigen melanoma-associated antigen-A3/A6 is a novel target of fibroblast growth factor receptor 2-IIIb through histone H3 modifications in thyroid cancer. Clin Cancer Res. 2007;13(16):4713–20. https://doi.org/10.1158/1078-0432.CCR-07-0618.

    CAS  Article  PubMed  Google Scholar 

  39. Guo M, Liu W, Serra S, Asa SL, Ezzat S. FGFR2 isoforms support epithelial-stromal interactions in thyroid cancer progression. Cancer Res. 2012;72(8):2017–27 https://doi.org/10.1158/0008-5472.CAN-11-3985.

    CAS  PubMed  Article  Google Scholar 

  40. Ren H, Luo M, Chen J, Zhou Y, Li X, Zhan Y, et al. Identification of TPD52 and DNAJB1 as two novel bile biomarkers for cholangiocarcinoma by iTRAQ-based quantitative proteomics analysis. Oncol Rep. 2019;42(6):2622–34 https://doi.org/10.3892/or.2019.7387.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. Behrends C, Sowa ME, Gygi SP, Harper JW. Network organization of the human autophagy system. Nature. 2010;466(7302):68–76 https://doi.org/10.1038/nature09204.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  42. Cui X, Choi HK, Choi YS, Park SY, Sung GJ, Lee YH, et al. DNAJB1 destabilizes PDCD5 to suppress p53-mediated apoptosis. Cancer Lett. 2015;357(1):307–15 https://doi.org/10.1016/j.canlet.2014.11.041.

    CAS  PubMed  Article  Google Scholar 

  43. Kawai T, Nomura F, Hoshino K, Copeland NG, Gilbert DJ, Jenkins NA, et al. Death-associated protein kinase 2 is a new calcium/calmodulin-dependent protein kinase that signals apoptosis through its catalytic activity. Oncogene. 1999;18(23):3471–80 https://doi.org/10.1038/sj.onc.1202701.

    CAS  PubMed  Article  Google Scholar 

  44. Shiloh R, Gilad Y, Ber Y, Eisenstein M, Aweida D, Bialik S, et al. Non-canonical activation of DAPK2 by AMPK constitutes a new pathway linking metabolic stress to autophagy. Nature Commun. 2018;9(1):1759 https://doi.org/10.1038/s41467-018-03907-4.

    Article  CAS  Google Scholar 

  45. Zhou Q, Wang L, Chen H, Xu B, Xu W, Sheng Y, et al. 2,3′,4,4′,5-Pentachlorobiphenyl induced autophagy of the thyrocytes via DAPK2/PKD/VPS34 pathway. Arch Toxicol. 2019;93(6):1639–48 https://doi.org/10.1007/s00204-019-02458-x.

    CAS  PubMed  Article  Google Scholar 

  46. Jiang Y, Liu J, Xu H, Zhou X, He L, Zhu C. DAPK2 activates NF-κB through autophagy-dependent degradation of I-κBα during thyroid cancer development and progression. Ann Transl Med. 2021;9(13):1083 https://doi.org/10.21037/atm-21-2062.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. Liao M, Shen J, Zhang Y, Li SH, Li XJ, Li H. Immunohistochemical localization of huntingtin-associated protein 1 in endocrine system of the rat. J Histochem Cytochem. 2005;53(12):1517–24. https://doi.org/10.1369/jhc.5A6662.2005.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. Wang T, Wang J, Wang J, Mao L, Tang B, Vanderklish PW, Liao X, Xiong ZQ, Liao L. HAP1 is an in vivo UBE3A target that augments autophagy in a mouse model of Angelman syndrome. Neurobiol Dis. 2019;132:104585. https://doi.org/10.1016/j.nbd.2019.104585.

  49. Wong YC, Holzbaur EL. The regulation of autophagosome dynamics by huntingtin and HAP1 is disrupted by expression of mutant huntingtin, leading to defective cargo degradation. J Neurosci. 2014;34(4):1293–305. https://doi.org/10.1523/JNEUROSCI.1870-13.2014.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  50. Zhu L, Song X, Tang J, Wu J, Ma R, Cao H, et al. Huntingtin-associated protein 1: a potential biomarker of breast cancer. Oncol Rep. 2013;29(5):1881–7 https://doi.org/10.3892/or.2013.2303.

    CAS  PubMed  Article  Google Scholar 

  51. Kakolyris S, Kaklamanis L, Giatromanolaki A, Koukourakis M, Hickson ID, Barzilay G, et al. Expression and subcellular localization of human AP endonuclease 1 (HAP1/Ref-1) protein: a basis for its role in human disease. Histopathol. 1998;33(6):561–9 https://doi.org/10.1046/j.1365-2559.1998.00541.x.

    CAS  Article  Google Scholar 

  52. Roy D, Sin SH, Damania B, Dittmer DP. Tumor suppressor genes FHIT and WWOX are deleted in primary effusion lymphoma (PEL) cell lines. Blood. 2011;118(7):e32–9 https://doi.org/10.1182/blood-2010-12-323659.

    PubMed  PubMed Central  Article  Google Scholar 

  53. Messai Y, Noman MZ, Hasmim M, Janji B, Tittarelli A, Boutet M, et al. ITPR1 protects renal cancer cells against natural killer cells by inducing autophagy. Cancer Res. 2014;74(23):6820–32 https://doi.org/10.1158/0008-5472.CAN-14-0303.

    CAS  PubMed  Article  Google Scholar 

  54. Peng D, Li W, Zhang B, Liu X. Overexpression of lncRNA SLC26A4-AS1 inhibits papillary thyroid carcinoma progression through recruiting ETS1 to promote ITPR1-mediated autophagy. J. Cell. Mol. Med. 2021;25(17):8148–58 https://doi.org/10.1111/jcmm.16545.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. Mijanović O, Branković A, Panin AN, et al. Cathepsin B: A sellsword of cancer progression. Cancer Lett. 2019;449:207–14. https://doi.org/10.1016/j.canlet.2019.02.035.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. Kim EK, Song MJ, Jang HH, Chung YS. Clinicopathologic analysis of cathepsin B as a prognostic marker of thyroid cancer. Int. J. Mol. Sci. 2020;21(24):9537 https://doi.org/10.3390/ijms21249537.

    CAS  PubMed Central  Article  Google Scholar 

  57. You L, Wang Z, Li H, Shou J, Jing Z, Xie J, et al. The role of STAT3 in autophagy. Autophagy. 2015;11(5):729–39 https://doi.org/10.1080/15548627.2015.1017192.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  58. Han H, Zhou H, Li J, Feng X, Zou D, Zhou W. TRAIL DR5-CTSB crosstalk participates in breast cancer autophagy initiated by SAHA. Cell Death Dis. 2017;3:17052 https://doi.org/10.1038/cddiscovery.2017.52.

    CAS  Article  Google Scholar 

  59. Zhang L, Wang Y, Li X, Wang Y, Wu K, Wu J, et al. Identification of a recurrence signature and validation of cell infiltration level of thyroid cancer microenvironment. Front Endocrinol. 2020;11:467 https://doi.org/10.3389/fendo.2020.00467.

    Article  Google Scholar 

  60. Wang W, Yang Z, Ouyang Q. A nomogram to predict skip metastasis in papillary thyroid cancer. World J Surg Oncol. 2020;18(1):167 https://doi.org/10.1186/s12957-020-01948-y.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  61. Feng JW, Hong LZ, Wang F, Wu WX, Hu J, Liu SY, et al. A nomogram based on clinical and ultrasound characteristics to predict central lymph node metastasis of papillary thyroid carcinoma. Front Endocrinol. 2021;12:666315 https://doi.org/10.3389/fendo.2021.666315.

    Article  Google Scholar 

  62. Ruchong P, Haiping T, Xiang W. A five-gene prognostic nomogram predicting disease-free survival of differentiated thyroid cancer. Dis Markers. 2021;2021:5510780 https://doi.org/10.1155/2021/5510780.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  63. Tang J, Jiang S, Gao Q, Xi X, Gao L, Zhao R, et al. Development and validation of a nomogram based on stromal score to predict progression-free survival of patients with papillary thyroid carcinoma. Cancer Med. 2021;10(16):5488–98 https://doi.org/10.1002/cam4.4112.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We greatly appreciate support by the major subject fund of Zhongnan Hospital of Wuhan University (XKJS202015).

Funding

This work was supported by the major subject fund of Zhongnan Hospital of Wuhan University (XKJS202015).

Author information

Affiliations

Authors

Contributions

CX Li and GS Wu designed the study, and CX Li, QQ Yuan, and GS Wu wrote and revised the manuscript. All authors participated in the data analysis. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Gaosong Wu.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors have no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visithttp://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, C., Yuan, Q., Xu, G. et al. A seven-autophagy-related gene signature for predicting the prognosis of differentiated thyroid carcinoma. World J Surg Onc 20, 129 (2022). https://doi.org/10.1186/s12957-022-02590-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/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