Abstract
Background/Aim: Thyroid carcinoma (THCA) is a cancer of the endocrine system that most commonly affects women. Aging-associated genes play a critical role in various cancers. Therefore, we aimed to gain insight into the molecular subtypes of thyroid cancer and whether senescence-related genes can predict the overall prognosis of THCA patients. Materials and Methods: Thyroid carcinoma (THCA) transcriptome-related expression profiles were obtained from The Cancer Genome Atlas (TCGA) database. These profiles were randomly divided into training and validation subsets at a ratio of 1:1. Unsupervised clustering algorithms were used to compare differences between the two subtypes; prognosis-related senescence genes were used to further construct our prognostic models by univariate and multivariate Cox analyses and construct a nomogram to predict the 1-, 3-, and 5-year overall survival probability of THCA patients. In addition, we performed gene set enrichment analysis (GSEA) to predict the immune microenvironment and somatic mutations between the different risk groups. Finally, real-time PCR was used to verify the expression levels of key model genes. Results: The ‘ConsensusClusterPlus’ R package was used to cluster thyroid cancer into two categories (Cluster1 and Cluster2) on the basis of 46 differentially expressed aging-related genes (DE-ARGs); patients in Cluster1 demonstrated a better prognosis than those in Cluster2. Cox analysis was used to screen six prognosis-related DE-ARGs. Finally, our real-time PCR results confirmed our hypothesis. Conclusion: Differences exist between the two subtypes of thyroid cancer that help guide treatment decisions. The six DE-ARG genes have a high predictive value for risk stratifying THCA patients.
Thyroid cancer, also known as THCA, ranks seventh among the most prevalent cancers in women globally. It comprises approximately 3% of all new cancer cases in women worldwide. It was estimated that around 31,180 female individuals would have this disease in 2022 (1). As quality of life has improved over the past few decades, the incidence of thyroid carcinoma has rapidly increased (2). This cancer has multiple pathological subtypes, most of which are papillary. Standardized pharmacological, surgical, and radioiodine treatments have resulted in longer survival times; however, several patients still experience early-onset distant metastasis with high lethality (3, 4). The current treatment for this small group of papillary thyroid carcinoma (PTC) patients remains inadequate. Therefore, we aimed to find reliable prognostic markers for thyroid cancer patients through a complete bioinformatics approach to enhance surveillance and early prevention and maximize survival times.
In recent years, with the diversification of research tools, substantial progress has been made in the study of aging-related genes, and more is known about their structure, function, and complex network regulatory mechanisms. Malignant tumor cells have three distinct characteristics that differentiate them from normal cells: unlimited proliferation, loss of contact inhibition, and uncontrolled invasion (5). This unlimited proliferation is often referred to as immortality. All normal regulated cells in the human body have a certain life cycle, and once they reach a certain threshold, cellular senescence and death occurs (6). However, malignant cells gain unlimited proliferation capacities due to impaired aging regulation and lack of programmed cell death (7, 8). Therefore, we can easily associate thyroid cancer with genes related to the regulation of senescence. Aging has long been a research focus worldwide, and researchers have broadly summarized 12 hallmarks of aging, including genomic instability and cellular senescence (9, 10), that are intricately related to tumor development and show extremely strong two-sidedness (11, 12). On the one hand, chronic inflammatory mediators produced by senescent cells induce epithelial-mesenchymal transition and stimulate the proliferation of malignant epithelial cells (13). On the other hand, senescent cells enhance immune surveillance and provide a solid foundation for the clearance of premalignant tumor cells by immune cells (14). Surprisingly, senescent cancer cells enhance tumor aggressiveness and lymph node metastasis in thyroid carcinoma mice (15). However, the exact invasive mechanism is unknown.
Many prediction models based on senescence-related genes have been developed to predict patient prognosis and drug treatment effects, indicating that senescence-related genes require further attention. However, few studies have explored the molecular subtypes of thyroid cancer and aging-related genes, and no reliable model has been developed to precisely predict the prognosis of thyroid cancer patients based on aging-related genes. Therefore, we aimed to comprehensively analyze the potential association between aging-related genes and thyroid carcinoma through bioinformatic approaches to identify poor prognostic factors and maximize patient life expectancy.
Materials and Methods
Data availability. Our research workflow is depicted in Figure 1. Thyroid carcinoma (THCA) RNA-seq samples and clinical information were downloaded from The Cancer Genome Atlas (TCGA) [GDC (cancer.gov.)], comprising 512 tumor samples. Eight samples with incomplete clinical information were excluded, leaving 504 THCA samples and 59 normal samples for the follow-up analysis. Single nucleotide mutation information was downloaded from the TCGA, which recorded the mutation information of the 498 samples. 307 aging-related genes were obtained from the Human Aging Genomic Resources (HAGR) database (16). The whole TCGA-THCA cohort was randomly divided into a training set and a validation set at a ratio of 1:1, with 252 samples in each set.
The overall pipeline of the project. After analyzing the HAGR database, which contains genes related to aging, and downloading the TCGA thyroid cancer dataset, we performed differential analysis and molecular typing to identify 46 differential genes (DEGs). Using this information, we constructed a prognostic model. Furthermore, we evaluated the high-risk and low-risk groups for differences in immune function and somatic mutations. Ultimately, we validated our hypotheses by qPCR. HAGR: Human Ageing Genomic Resources; TIDE: Tumor Immune Dysfunction and Exclusion; qPCR: quantitative real-time PCR; GSEA: Gene Set Enrichment Analysis.
Analysis of differentially expressed aging-related genes (DE-ARGs). First, the limma package was used for the differential analysis of gene expression in our normal and tumor groups to screen DE-ARGs with a|log2-fold change (FC)|>2 and false discovery rate (FDR) <0.05. Second, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed for the above DE-ARGs to explore their potential regulatory pathways. Protein-protein interactions (PPIs) were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (17), and hub genes were selected by the “cytohubba” plugin in Cytoscape (18). Using the ‘corrplot’ R package (version 0.92), we focused on the correlation of the top 10 up- and down-regulated genes.
Molecular subtype analysis. To further reveal the potential expression patterns and molecular subtypes of aging-related genes in thyroid carcinoma, we applied an unsupervised clustering method with 100 X resampling using the above-screened aging-related differential genes and the K-Means algorithm and Euclidean distance parameters in the “ConsensusClusterPlus” R package (19). The best clustering pattern (k=2) was selected by the consensus clustering matrix and cumulative distribution function (CDF). The corresponding samples were clustered into two clusters. Survival analyses were performed on the clustered samples, and clinically relevant information was further explored.
Analysis of the immune landscape between the two subtypes. Since immunity is irreplaceable in tumorigenesis and progression, we focused on immune differences between the two clusters. Human leukocyte antigen (HLA) expression on antigen-presenting activated CD8+ T cells is critical in predicting patient response to immunotherapy (20). Therefore, we compared HLA expression differences between the two clusters. The Estimation of Stromal and Immune cells in Malignant Tumor tissues using the Expression data (ESTIMATE) algorithm (21) was used to determine the differences between tumor and stromal cells in both clusters. Although immunotherapy is not routine for thyroid carcinoma, immunotherapy shows substantial promise in oncology overall. Immunotherapy for iodine-refractory differentiated thyroid carcinoma is a potential future approach. Therefore, six common immune checkpoints (CTLA4, LAG3, HAVCR2, TIGIT, PD-L1, and PD1) and Tumor Immune Dysfunction and Exclusion (TIDE) (22) scores were included in the cluster analysis.
Building Aging Risk Scoring (ARS) prognostic models and nomogram construction. A risk-scoring model based on clustering was further developed, accounting for the individual heterogeneity and diversity within THCA. Our initial step was to identify genes associated with prognosis by univariate Cox regression using a threshold of p<0.05, followed by a multivariate Cox regression analysis to identify key DE-ARGs. Through this approach, we constructed a predictive model for THCA scores. The specific scoring formula was as follows:
where n represents the number of moduleRNAs (i.e., DEGs); Coef (i) is the coefficient (weight) for each gene, obtained through multivariable Cox regression; X(i) indicates the relative expression level of the z score transformation of each gene determined by multivariate analysis. The survival differences between the two risk stratification groups were determined based on Kaplan–Meier (K–M) survival curves. The training and validation sets were plotted separately to obtain receiver operating characteristic (ROC) curves for model efficacy validation to reveal the stability of the prediction model. Our prediction model and clinically relevant risk factors (T, N, M, age, stage, sex) were then combined and presented as a nomogram to identify the risk factors that could affect prognosis. Calibration plots were used to predict survival probabilities at 1, 3, and 5 years.
Gene set enrichment analysis (GSEA) based on the predictive model. GSEA was used to determine whether certain genes are indeed differentially expressed between the two risk groups. GSEA demonstrated the potential pathway mechanisms between the two sets by determining the enrichment difference between a predefined set of background genes and the sequenced phenotype-associated genes in the risk subgroup. Significantly changed pathways were identified by using 1,000 permutation tests. Adjusted p-values <0.05 were considered significant. Background annotation files ‘c2.cp.kegg.v2022.1.Hs.symbols.gmt’ and ‘c5.go.v2022.1.Hs. symbols.gmt’ were obtained from the MSigDB database (23).
Analysis of high- and low-risk groups with immune landscapes. Tumor progression is known to be inextricably linked to the immune microenvironment (24). Therefore, we explored the immune landscape among different risk groups using the RNA expression matrix of the whole TCGA-THCA cohort. We explored the differences in stromal cells, tumor purity, and 22 immune-related cells in our high- and low-risk individuals using the ESTIMATE and CIBERSORT (25) algorithms. Meanwhile, using the GSVA package (26) in R, both groups were compared in terms of immune cell characteristics and immune function. An in-depth analysis of the differences in HLA expression between high- and low-risk individuals is necessary to understand how HLA exerts an irreplaceable role in immune surveillance and immune escape. In addition, we were curious about the response of THCA to immunotherapy; therefore, 47 common immune checkpoints and TIDE were also included in the analysis.
Gene mutation analysis of the risk score model. BRAF gene mutations are known to contribute to thyroid carcinoma development (27). Thus, we attempted to clarify the predictive value of somatic mutation-associated genes in thyroid cancer. Somatic mutation data were obtained from the TCGA database for 498 samples. We calculated the tumor mutation burden (TMB) differences across samples with the R package ‘maftools’ (28) and selected the top 20 genes with the highest mutation frequencies. Subsequently, survival analyses were plotted using the TMB and risk scores. In addition, we were keenly interested in the mutation information of our model genes; thus, we assessed these genes using the cBioPortal database (29).
Cell culture and real-time polymerase chain reaction (RT-qPCR). To further illustrate the accuracy of our analysis, a normal thyroid cell line (Nthy ori-3-1) and a thyroid cancer cell line (TPC-1) were purchased from Shanghai, PR China. Nthy ori-3-1 cells were cultured in RPMI-1640 medium, while TPC-1 cells were maintained in 10% high glucose DMEM. Both were successfully passaged at 5% CO2 in a 37°C incubator.
We extracted total RNA from Nthy-ori-3-1 and TPC-1 cells using a column extraction reagent to validate the relative expression of the six genes (NGFR, IL6, CLU, CDKN2, HSP-7, and E2F1). GAPDH was used as the reference gene. For the amplification of the mentioned genes, the primer sequences were as following: GAPDH, forward, 5′-GTCTCCTCTGACTTCAACAGCG-3′ and reverse, 5′-ACCACCCTGTTGCTGTAGCCAA-3′; E2F1, forward, 5′-GCTGGACCACCTGATGAA-3′ and reverse, 5′-GGAGGGG CTTTGATCACC-3′; IL6, forward, 5′-ACAGCCACTCACCT CTTCAG-3′ and reverse, 5′-CCATCTTTTTCAG CCATCTTT-3′; CLU, forward, 5′-TGCCTGAAACAGACCTGCAT-3′ and reverse, 5′-CATGACATCCAGCATGTGCG-3′; NGFR, forward, 5′-TCTG ACGTCACAACACCCAG-3′ and reverse, 5′-AGCAGCCAAGAT GGAGCAAT-3′; CDKN2A, forward, 5′-CCCGCTTTCGTAG TTTTCAT-3′ and reverse, 5′-TTATTTGAGCTTTGGTTCTG-3′; HSP-7, forward, 5′-AGCGCAAGTACAAGAAGGAC-3′ and reverse, 5′-TGGTGATGGAGGTGTAGAAG-3′. The fold differences in gene expression levels in the thyroid cancer versus normal cell lines were calculated using normalized CT values.
Results
Identification of DE-ARGs and functional enrichment analysis. We used the limma package in R to obtain 46 differentially expressed genes, 23 up-regulated and 23 down-regulated (Figure 2A). Subsequently, the above differentially expressed genes were correlated in the STRING database and further visualized in Cytoscape, which showed interactions between 41 proteins (Figure 2B). In addition, the above 46 DE-ARGs were assessed using KEGG and GO functional enrichment analyses and found to be mostly concentrated in signaling pathways, such as immune, cell proliferation, and hormonal regulation (Figure 2C and D). The top 10 up- and down-regulated genes were extracted separately to explore their correlations (Figure 2E). The Cytoscape cytoHubba plugin was used to visualize the top 10 hub genes (Figure 2F).
Differentially expressed gene identification and analysis. A: Volcano plot showing the 46 differentially expressed aging-related genes (FDR <0.05 and |logFC| >2). B: PPI network visualization in Cytoscape. C: The top 30 potential pathways enriched in the KEGG analysis. D: BP, CC, and MF were shown to be enriched in the top 10 pathways. E: Correlation analysis of the top 10 up-regulated and down-regulated genes, plotted using the ‘pheatmap’ R package. F: Ten key aging-related genes obtained by the cytoHubba plugin using the cytoscape software.
Identification and evaluation of DE-ARGs molecular subtypes. The “ConsensusClusterPlus” R package with K=2 revealed that the 46 DE-ARGs divided our sample into 2 classes (Cluster1 and Cluster2, Figure 3A and B) separated by a clear boundary, indicating that the two subtypes may differ significantly in some aspects. Survival curves and clinical characteristics were further assessed to highlight the differences between the two subtypes. The survival curves revealed that Cluster1 patients had significantly better overall survival outcomes than those in Cluster2 (Figure 3C). In addition, some similar differences were observed in the clinical characteristic distribution across the two subtypes (Figure 3D). Therefore, immunotherapy for aggressive THCA is a potentially viable approach.
Clustering analysis of DE-ARGs in thyroid carcinoma. A: Clustering matrix with K=2. B: Consensus CDF when K is clustered in the range of 2 to 9. C: Survival curves for Cluster1 and Cluster2. D: The heatmap was generated using the ‘pheatmap’ R package and demonstrates the differences in clinically relevant factors between Cluster1 and Cluster2.
Relationship between the two subtypes and the immune landscape. HLA genes are critical in immune surveillance and responses; thus, we compared the differences in HLA-related genes between the two clusters. We included 39 HLA-associated genes, most of which were significantly differentially expressed (Figure 4A), and further assessed the tumor microenvironment variation between Cluster1 and Cluster2 by using the ESTIMATE algorithm, which typically comprises three scores: ImmuneScore, StromalScore, and TumorPurity (Figure 4B). The results confirmed a higher abundance of immune infiltration in Cluster 1, again in accordance with our prognostic analysis. Finally, to assess the immunotherapy response, we included well-known immune checkpoints, such as CD274 (PD-L1), PDCD1 (PD-1), LAG3, TIGIT, HAVCR2, and CTLA4. The expression levels of CTLA4, LAG3, TIGIT, and HAVCR2 were higher in the Cluster1 subtype (Figure 4C). The findings suggest that immunosuppression occurs more frequently in Cluster1 subtype THCA with high expression of immune checkpoint inhibitors. It is speculated that this may be due to the high level of immune cell infiltration in Cluster1 and uneven distribution among immune cells. The TIDE score also supports this hypothesis (Figure 4D).
Immunological differences between Cluster1 and Cluster2 and response to immunotherapy. A: Significance comparison of HLA genes between Cluster1 and Cluster2. *p<0.05, **p<0.01, ***p<0.001. B: Immune infiltration analysis of the two clusters. C: Differences in the six common immune checkpoints between Cluster1 and Cluster2; D: TIDE score predicting the response to immunotherapy.
Construction of a prognostic model consisting of six DE-ARGs. The 46 DE-ARGs were assessed using univariate and multivariate Cox analyses to accurately risk-stratify PTC patients, resulting in six potential associations with aging predictors, including NGFR, CDKN2A, CLU, E2F1, IL-6, and HSP-7. the risk equation was constructed as follows: (−0.30481 * E2F1 expression) + (−0.08585 * NGFR expression) + (0.061243 * IL6 expression) + (0.00924 * HSPA1A) + (−0.00083 * CLU expression) + (0.139632 * CDKN2A expression). Furthermore, 504 tumor samples in the TCGA-THCA cohort were successfully stratified into two distinct risk groups according to the median risk value. The risk curve demonstrated a gradual increase in patient mortality as the risk score increased (Figure 5A). The THCA-training and THCA-test cohort results supported this observation (Figure 5B and C). The ROC curves’ ability to predict THCA-all cohort PTC patient 1-, 3-, and 5-year survival outcomes were 0.755, 0.872, and 0.902 (Figure 5D) using the R package ‘timeROC’; in the TCGA-training and TCGA-test subsets, these values were 0.850, 0.855 and 0.889 and 0.909, 0.880 and 0.969, respectively (Figure 5E and F). The above results showed that the model’s predictive ability using aging-related genes was very stable. K–M plots for the three cohorts showed that patients at high risk had significantly lower survival rates than those at low risk (Figure 5G and I). Hazard ratios for six key genes were visualized by the ggplot2 R package (Figure 6A). Clinically relevant factors, including age, sex, clinical stage, T stage, N stage, and M stage, were included in the univariate analysis to further clarify whether our risk scores could be used as independent predictors. The results showed that age, clinical stage, T stage, and the risk score were considered significant (p<0.05, Figure 6B). Meaningful clinical correlates in the above univariate analysis results were then included in the multivariate analysis, revealing two independent prognostic factors, age and the risk score (Figure 6C). The risk score was positively correlated with prognosis. In addition, we constructed a nomogram combining age and the risk scores to assess overall survival at 1, 3, and 5 years (Figure 6D). The calibration curves show that the nomogram has superior predictive performance (Figure 6E).
Prognostic model construction and model efficacy evaluation of the six DE-ARGs in thyroid carcinoma. A-C: Risk curves for high- and low-risk groups and the corresponding scatter plots in the TCGA-all, TCGA-training, and TCGA-test cohorts. D-F: ROC curves for prognostic models in the TCGA-all, TCGA-training, and TCGA-test cohorts. G-I: Analysis of survival differences in two different risk groups (left: TCGA-all cohort; middle: TCGA-training cohort; right: TCGA-test cohort).
Analysis of independent prognostic factors and nomogram construction. A: Visualizing the risk ratios of six key genes. B: Univariate analysis of risk scores combined with clinically relevant prognostic factors (age, sex, stage, T, N, and M). C: Results of multivariate Cox regression analysis. D: Nomogram predicting survival at 1, 3, and 5 years. E: Calibration curves showing the difference between the nomogram predicted and actual observations at 1, 3, and 5 years.
GSEA analysis in high- and low-risk groups. To fully understand the differences in potential pathway mechanisms between the two risk groups, based on the risk score prediction models, GSEA to assess the KEGG and GO pathways of the genes of interest to find the top 5 up- and down-regulated pathways in each group. The screening was conducted with a p-value <0.05 and FDR q-value <0.25. GSEA revealed that most of the genes in both groups were concentrated in metabolism-related signaling pathways. In the low-risk group, biosynthesis and ion transport were more abundant, while the cell cycle and DNA repair were more prevalent in the high-risk group (Figure 7A).
GSEA and immune landscape exploration between high- and low-risk individuals. A: GSEA screened the top 5 up-regulated and down-regulated pathways. B-C: Differences between the two groups in terms of immune cells and immune function (left: immune cells). There were no significant differences observed when examining B cells, CD8+ T cells, T helper cells, and Tfh cells. D: Twenty-nine HLA-associated genes showed significant differences between the two risk groups. The differences for the 29 listed HLA-related genes were statistically significant. E: Twenty-eight immune checkpoint-associated genes showed high expression in the low-risk group. All differences were statistically significant. *p<0.05, **p<0.01, ***p<0.001.
Evaluation of immune differences in high and low-risk groups. Our first step was to perform a broad-spectrum study of 22 immune cells in the risk subgroup population to assess potential immune differences In addition, we examined the function of immune cells and found that the low-risk group exhibited a higher level of activity in these cells (Figure 7B and C), consistent with their longer survival times. Furthermore, we found significantly lower expression levels of most MHC molecules in the high-risk than in the lower-risk population (Figure 7D). Additionally, a significant up-regulation of 28 immune checkpoints was present among 47 immune checkpoints examined in our low-risk cohort (Figure 7E), potentially indicating higher immune escape as well as a lower immunotherapy benefit.
Analysis of somatic mutations and TMB in high and low-risk groups. The top 20 most frequently mutated genes in the THCA groups were mapped separately using the R package “maftools” to better understand the intrinsic link between the risk models and gene mutations. In both groups, Braf mutations topped the list (Figure 8A and B). Neither risk group had a significant difference in TMB (p=0.26, Figure 8C). In addition, we explored the relationships between survival and TMB and survival and both scoring models and TMB values (Figure 8D and E). Both results showed a positive correlation between low TMB and better survival outcomes in thyroid carcinoma. Finally, visualization of the broad mutation landscape of our six model genes on the cBioPortal website revealed that our model genes had a lower overall mutation rate (Figure 8F).
Analysis of tumor mutation burden (TMB) and survival differences between the two groups. A: Mutations in the top 20 genes associated with high risk. B: Mutations in the top 20 genes associated with low risk. C: There was no difference in TMB values between the high- and low-risk groups. D: The high-risk group has a poorer prognosis than the low-risk group. E: The combination of TMB and the risk score model predicts survival in thyroid carcinoma. F: Visualization of the overall mutation rate of the six model genes.
Relative expression of NGFR, IL6, CLU, CDKN2A, HSP-7, E2F1 based on RT-qPCR. The above findings were further validated at the RNA level, and the RT–qPCR results demonstrated that four signature genes (NGFR, CDKN2A, CLU, and E2F1) were up-regulated in tumor samples, and the remaining two (IL-6 and HSP-7) were down-regulated (Figure 9), consistent with our previous analysis.
Validation of the expression of the 6 key model genes in normal (Nthy-ori-3-1) and tumor (TPC-1) cell lines by RT–qPCR. The statistical analysis of the data was conducted using a t test for two independent samples (n=3 per group). A: NGFR; B: HSP-7 (HSP1A1); C: IL-6; D E2F; E: CDKN2A; F: CLU. *p<0.05, **p<0.01, ***p<0.001 and ****p<0.0001.
Discussion
In recent years, with the popularity of ultrasound imaging and fine-needle aspiration (FNA) techniques, the incidence of thyroid carcinoma increased linearly and has been detected more frequently in younger patients (30), especially young women (31). According to the current understanding, thyroid cancer is generally associated with a low mortality rate. This is primarily due to a combination of factors (32). Firstly, thyroid cancer is often detected at an early stage, allowing for early initiation of treatment, and increasing the chances of a cure. Secondly, thyroid cancer typically exhibits slow growth and is largely confined to the thyroid gland, making surgical removal the primary treatment approach, and often enabling complete removal of the affected tissue. However, despite the relatively low overall mortality rate, there are still challenges that some patients face, such as disease recurrence, lymph node metastasis, and the development of distant metastases (33, 34). Therefore, early diagnosis and prompt treatment are of critical importance. The current challenge is to stratify thyroid carcinoma patients by risk and identify individuals with poor prognoses early in the hope of achieving individualized and standardized treatment (35, 36).
Aging is a constant research focus in the analysis of the progression and prognosis of various tumors. Aging is a complex, holistic, and stage-specific physiological process involving numerous genes, complex regulatory networks, and influencing factors. Previous studies have noted that senescent cells usually have three main characteristics: cell proliferation arrest, resistance to apoptosis, and the senescence-associated secretory phenotype (SASP) (37). Senescent cells become more active in aged or diseased tissues (38). In addition, senescent cells can secrete senescence-related secretory phenotypes that promote systemic inflammation (39) and increase the side effects of chemotherapy drugs leading to tumor recurrence (40). Senescent cells can also accelerate the continuous accumulation of somatic mutations that eventually lead to the onset of malignant tumors (41). However, in BRAF mutation-driven thyroid cancer, we found that senescent thyroid cancer cells appear at the front of tumor invasion (42). Meanwhile, we also observed that oncogene-induced senescence (OIS) in thyroid cancer promotes the malignant progression of papillary thyroid cancer (43, 44).
Senescence plays a crucial role in various tumors, and its pathogenic mechanism in thyroid carcinoma is unclear. Therefore, in our study, we divided the samples into two classes using 46 DE-ARGs. Compared to Cluster2, Cluster1 exhibited significant immune cell enrichment, and the ESTIMATE algorithm illustrated that Cluster1 samples had lower tumor purity. In contrast, the TIDE score showed a higher rate of immune escape in Cluster1, probably due to immune disorders caused by the higher immune cell abundance in this cluster and the uneven distribution of immune cells. Therefore, Cluster1 THCA patients are less likely to benefit from immunotherapy. We also risk-stratified the sample and constructed a nomogram with risk factors, such as age and the risk score to predict the overall survival rate at 1, 3, and 5 years. The calibration curves demonstrated that the nomogram showed high prediction accuracy. Finally, we performed mutation analyses for both risk groups, including the visualization of mutation information for six key genes. The tumor mutational burden was combined with the risk score model to better guide individualized clinical counseling and treatment.
BRAF gene mutations are prominent in thyroid cancer, with an incidence exceeding 40%, followed closely by N-RAS and H-RAS mutations (45), consistent with the analysis of the top 3 mutations in our high-risk group. Proteins driven by mutations in the Braf gene encode proteins that activate downstream mitogen-activated protein kinase (MAPK) signaling pathways, promoting malignant tumor progression and recurrence risk (46, 47). High tumor mutational burden tends to be positively associated with a poorer prognosis (48), as evidenced by the K–M survival curves based on differences in the TMB values. The combined TMB value risk model predicted that high-risk groups with concomitant high tumor mutation burden have the worst prognosis, which may help clinicians accurately stratify patients and identify individuals with poor prognoses early for individualized treatment.
Among the 6 genes involved in constructing the prognostic model, 4 genes (E2F1, IL6, NGFR, and CDKN2A) were considered to dominate the risk score. E2F1 activates transcription and is essential for regulating the cell cycle, apoptosis, and proliferation by targeting DNA promoters (49). E2F1 has recently been implicated in numerous tumor types, including breast cancer (50) and retinoblastoma (51). Interleukin-6 (IL-6) is remarkably active in immunity, inflammation, and metabolism (52). Studies have confirmed that IL-6 regulates kinase activity by binding to the IL-6 receptor to activate gp130 during inflammation and immune activation (53). Nerve growth factor receptors (NGFR) are involved in cell death and the induction of fetal myogenic progenitor cell differentiation (54). In addition, the up-regulation of NGFR gene expression in melanoma increases the risk of early tumor lymph node metastasis and resistance to immunotherapy (55, 56). CDKN2A is considered a tumor suppressor-associated gene (57) that functions through MDM2 to inhibit P53 degradation, thereby inhibiting cell proliferation and accelerating apoptosis (58). Meanwhile, mutation or deletion of CDKN2A often predicts a poor prognosis (59). As it is widely acknowledged, anaplastic thyroid cancer (ATC) represents a unique histological subtype of thyroid cancer that carries a high risk and frequently leads to distant metastases (60). Several studies have revealed prominent overexpression of enhancer of zeste homolog 2 (EZH2) in ATC cell lines (61), and inhibiting the expression of EZH2 has presented new avenues for treatment. Additionally, the hsp70 gene has the capability to promote tumor cell death by enhancing the PI3K/Akt signaling pathway (62). Consequently, ATC patients can find hope in the potential of restoring their health through these advancements.
Study limitations. Firstly, it was conducted retrospectively, and we lacked an independent external validation set to corroborate our predictive model, which may restrict the generalizability of the model. Secondly, although there may be genes associated with aging that regulate thyroid cancer, further fundamental research is required to confirm this. Lastly, our genetic prognostic analysis was limited to the RNA level, and additional in-depth investigation into the impact on the protein level is necessary. Therefore, in the future, large-scale prospective studies, along with experimental research, will be needed to validate our hypotheses.
Conclusion
In conclusion, a preliminary exploration of the molecular subtypes of thyroid cancer based on 46 DE-ARGs and a prognostic model constructed using these genes were presented. Differences in the immune landscape and gene mutations between molecular subtypes and risk groups were also highlighted. The prognostic models constructed in this study may help to identify thyroid carcinoma patients with poor prognoses earlier in the disease process, allowing clinicians to treat them accordingly.
Acknowledgements
We thank the TCGA database and those who are willing to share information freely.
Footnotes
Conflicts of Interest
The Authors of this project do not have any conflicts of interest.
Authors’ Contributions
Zhi Li and Li jia contributed equally to this article. Zhi Li was responsible for the analysis of the data and the draft of the paper. Li Jia and Lu Zhang helped with the revision of the paper. Zhi-Yong Deng and Chao Liu offered the idea and framework of the paper. ZM, ZHR, BYK, and LJ participated in the discussion. All Authors read and approved the final manuscript.
- Received October 12, 2023.
- Revision received November 21, 2023.
- Accepted November 24, 2023.
- Copyright © 2024, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) 4.0 international license (https://creativecommons.org/licenses/by-nc-nd/4.0).
