Comprehensive analysis of cellular senescence and immune microenvironment in papillary thyroid carcinoma

Materials and Methods

SEER cohort, TCGA - THCA cohort, GTEx and GEO cohorts

We followed the flowchart outlined in Figure 1 for our analysis. The SEER database is an authoritative cancer statistical database in the United States, recording the incidence, mortality, and morbidity of millions of patients with malignant tumors in various states and counties. Our study included patients diagnosed with PTC as their first primary tumor in the SEER cohort between 2004 and 2015. The TCGA-THCA dataset provides comprehensive information, including transcriptome profiling, copy number variation, DNA methylation, and somatic structural variation for a large number of cancer patients. We included PTC cases from the TCGA-THCA cohort with complete follow-up information, and we retrieved clinical data from the University of California at Santa Cruz (UCSC) Xena website (https://xena.ucsc.edu/). Transcriptome data of PTC tissues and normal tissues were obtained from the TCGA (T=510, N=58), GTEx (N = 337) and eight GEO cohorts (http://www.ncbi.nlm.nih.gov/geo): GSE33630 (T=49, N=45), GSE60542 (T=33, N=30), GSE58545 (T=27, N=18), GSE3467 (T=9, N=9), GSE3678 (T=7, N=7), GSE5364 (T=35, N=16), GSE27155 (T=51, N=4), GSE53157 (T=15, N=3). Since we obtained the data from public source websites, no ethical review was required.

Figure 1. Flow chart of data analysis. (A) Prognosis by age group in SEERand TCGA-THCA cohorts; (B) Intersection of differentially expressed cellular senescence genes; (C) Prognosis-related cellular senescence genes; (D) Multi-omics analysis between two clusters; (E) Prognostic analysis between two clusters; (F) Consensus Cluster Analysis; (G) Construction of prognostic model of cellular senescence; (H) Analysis and validation of the performance of prognostic model of cellular senescence; (I) Multi-omics analysis.

Consensus clustering analysis based on genes related to cellular senescence

We extracted 278 genes related to cellular senescence from the CellAge website [21] (https://genomics.senescence.info/cells/) (Supplementary Table 1), which contains manually curated gene information related to cellular senescence. We analyzed the cell senescence genes extracted from PTC and normal tissues in TCGA, GSE33630, GSE60542 and GSE58545 through the “limma” package [22] of R software (R 4.1.0) to obtain differentially expressed genes (DEG) (FDR < 0.05), and then took the intersection [23]. The TCGA-THCA and GTEx data were batch corrected using the “normalizeBetweenArrays” function of the “limma” package. In cases where probe data shared the same gene name, we aggregated them to calculate the average expression. Subsequently, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses on the genes in the intersection. We also constructed a protein-protein interaction (PPI) network using the STRING website (https://string-db.org/) and Cytoscape software (version 3.8.2). Screening for prognostic genes associated with cellular senescence was done with a Cox analysis. Univariate Cox analysis was used to screen cellular senescence genes associated with OS (p < 0.05) and the expression levels of these genes were displayed with a heatmap.

In accordance with the expression of OS-related cellular senescence genes, we conducted a consensus cluster analysis on PTC patients and examined Kaplan-Meier (K-M) curves, clinicopathological characteristics, immune cell infiltration, immune microenvironment, and KEGG enrichment pathways among different clusters.

The abundance of immune cell infiltration was assessed using the CIBERSORT algorithm [24, 25], which can determine the abundance of 22 types of immune cells in the tissue. The tumor microenvironment was evaluated using four indicators (ESTIMATEScore, StromalScore, ImmuneScore and TumorPurity) using the ESTIMATE algorithm [26].

Establishing and validating a prognostic model

To quantitatively analyze the impact of cellular senescence genes on the prognosis of PTC patients, we randomly divided the TCGA-THCA cohort into a training cohort and a validation cohort (1:1). We constructed a CSRS by performing Least absolute shrinkage and selection operator (LASSO) regression analysis [27, 28] on the prognosis-related cellular senescence genes in the training cohort using the “glmnet” package [29]. Based on the prognostic model, we performed risk score calculations for patients in the training cohort, validation cohort, and the entire TCGA-THCA cohort. The calculation formula of the CSRS was as follows [28]:

$$CSRS={\displaystyle \sum _{i=1}^{n}Ex{p}_i\ast Coe{f}_i}$$

where Exp_i is the expression value of the i gene in the model, and Coef_i is the coefficient calculated by LASSO.

We conducted CSRS (Cellular Senescence Risk Score) calculations for patients within the training cohort, validation cohort, and the entire TCGA-THCA cohort. Based on the median risk scores, patients were stratified into high-risk and low-risk groups.

Subsequently, we performed a comprehensive array of analyses, including K-M survival curves, receiver operating characteristic (ROC) curves, decision analysis curves, survival status maps, and heatmaps illustrating the genes associated with CSRS in distinct patient groups.

The expression levels of genes involved in the construction of the CSRS were verified by the GEPIA2 [30] website (http://gepia2.cancer-pku.cn/#index) and in five GEO datasets: GSE3467 (T=9, N=9), GSE3678 (T=7, N=7), GSE5364 (T=35, N=16), GSE27155 (T=51, N=4), GSE53157 (T=15, N=3).

Hierarchical analysis

To further mitigate the potential influence of clinicopathological features on CSRS as a prognostic factor, we conducted a stratified prognostic analysis based on various clinicopathological features within the entire TCGA-THCA cohort. These features include age [≥55 vs. < 55 years], gender [male vs. female], AJCC stage [III - IV vs. I - II], M [M1 vs. M0], N [N1 vs. N0], T [T3 - 4 vs. T1 - 2], BRAF [mutated vs. wild-type], RAS [mutated vs. wild-type], Cell Type [Classical vs. Follicular], Radiation therapy [Yes vs. No].

Multi-omics analysis

In the entire TCGA-THCA cohort, using seven algorithms (TIMER [31], CIBERSORT [24, 25], CIBERSORT-ABS [24], QUANTISEQ [32], MCPCOUNTER [33], XCELL [34], EPIC [35]), and the analyzed algorithms will be employed to elucidate the extent of immune cell infiltration among different groups and identify immune cells playing pivotal roles. We created immune cell heatmaps to compare infiltration rates between high- and low-risk groups.

To further analyze the differences in pathway enrichment between high and low risk groups, we performed Gene Set Variation Analysis (GSVA) with the “GSVA” package [36]. The gene set of “c2.cp.kegg.v7.4.symbols.gmt” was obtained from the MSigDB website [37], and statistics were considered significant was defined as p-values less than 0.05.

Subsequently, we employed clinical correlation heatmaps to illustrate the distribution of clinicopathological features in the high and low-risk groups, along with the expression of genes involved in the CSRS.

The immunophenoscore (IPS) serves as a valuable indicator for predicting the efficacy of anti-CTLA-4 and anti-PD-1 immunotherapy in tumor patients [25]. We obtained the IPS score of PTC patients from The Cancer Immunome Atlas (https://tcia.at/) and used the violin plot to display differences in IPS scores between high and low risk groups.

A model for tumor immune evasion, known as Tumor Immune Dysfunction and Exclusion (TIDE: http://tide.dfci.harvard.edu/) can offer insights into the effects of immunotherapy on two primary mechanisms [38, 39]. We uploaded the gene matrix from TCGA-THCA to the TIDE website and generated violin plots depicting TIDE scores, Dysfunction scores, Exclusion scores and CAF scores between high and low risk groups.

Some advanced cancers have been treated better with immune checkpoint inhibitors (ICIs) [40]. Boxplots, based on 47 immune checkpoints obtained from the literature [41, 42], were used to illustrate the differential expression of immune checkpoints in high-risk and low-risk groups.

Based on information from the Genomics of Drug Sensitivity in Cancer (GDSC) database, an assessment of PTC patients’ response to chemotherapeutic and targeted treatment agents was conducted. The half maximal inhibitory concentration (IC50) of pharmacological agents was estimated with the “pRRophetic” R package, which has been extensively used in medical research [43].

Data availability statement

All the data come from public databases.

Results

Age is a prognostic indicator in SEER and TCGA-THCA cohorts

According to inclusion and exclusion criteria, we included a total of 42,756 patients in the SEER database. We divided patients into two groups based on age (<55 years, ≥55 years). When OS was the outcome event, patients aged 55 or older also had a worse prognosis than the younger group (Figure 2A). Considering the differences in incidence and mortality of thyroid cancer in different gender populations [44–46], we presented the OS data for the younger and older age groups stratified by gender. The results showed that there was a statistically significant difference in overall survival between the young and old groups of women and men (Figure 2B, 2C). Similar results were obtained when we used DSS as the outcome event (Figure 2D–2F). In the TCGA-THCA cohort, when we used OS, DSS, and PFS as outcome events, the older group also had a worse prognosis than the younger group (Figure 2G–2I).

Figure 2. Prognostic differences among different age groups in SEER and TCGA-THCA cohorts. OS between PTC with Age ≥55 years and Age <55 years in all (A), female (B) and male (C) in SEER. DSS between PTC with Age ≥55 years and Age <55 years in all (D), female (E) and male (F) in SEER. OS (G), DSS (H) and PFS (I) between PTC with Age ≥55 years and Age <55 years in TCGA-THCA.

The cluster with shorter overall survival could be distinguished based on cell senescence genes

278 cellular senescence genes expression were differentially analyzed between tumors and normal tissues in the TCGA-THCA, GSE58545, GSE60542 and GSE33630 cohorts. This analysis identified 76 genes that were common across these datasets (Figure 3A and Supplementary Tables 2, 3). GO functional enrichment analysis showed that 76 intersecting genes were mainly enriched in aging; cell aging; cellular senescence; DNA damage response, signal transduction by p53 class mediator; signal transduction in response to DNA damage (Figure 3B). GO functional enrichment analysis indicated that these genes were mainly involved in aging-related functions. KEGG pathway enrichment analysis showed that these genes were enriched in Bladder cancer, Cell cycle, Endocrine resistance, Kaposi sarcoma - associated herpesvirus infection, p53 signaling pathway (Figure 3C). The PPI network diagram of the proteins corresponding to the 76 genes was shown in Figure 3D. It should be noted that P53 protein has many protein interactions, and both GO function enrichment analysis and KEGG pathway enrichment analysis involved p53 signaling mediated regulatory signals. It may be that P53 is involved in cell senescence process in PTC.

Figure 3. Consensus cluster analysis. Intersection of differentially expressed cellular senescence genes across four cohorts (A). GO functional enrichment analysis (B) and KEGG pathway enrichment analysis (C) of 76 intersecting genes. Protein-Protein interactions (PPI) network diagram of 76 intersecting genes (D). Forest plot of prognostic genes among 76 intersecting genes (E). Expression heatmap of prognostic genes (F). Consensus clustering matrix for k = 2 (G). K-M curve of OS probability between cluster1 and cluster2 (H). Heatmap of clinicopathological feature correlations between cluster1 and cluster2 (I). Violin plots of the infiltration of immune cells in cluster1 and cluster2 by CIBERSORT algorithm (J). Boxplot of ESTIMATEScore (K), TumorPurity (L), ImmuneScore (M) and StromalScore (N) on cluster1 and cluster2 by ESTIMATE algorithm. KEGG pathway enrichment analysis in cluster1 (O) and cluster2 (P).

We performed univariate Cox analysis on the OS of PTC patients based on the expression of these 76 genes, and obtained 16 prognosis-related cellular senescence genes (Figure 3E and Supplementary Figure 2). The expression levels of these 16 genes (ASPH, NINJ1, SNAI1, UBTD1, MAP3K6, CDKN2A, E2F1, HDAC4, PLA2R1, PRKCD, WWP1, ID4, NEK1, P3H1, NDRG1, SORBS2) in tumor tissue and normal tissue in the TCGA-THCA cohort were displayed by heat map (Figure 3F). Using consensus clustering analysis based on the expression levels of prognosis-related cellular senescence genes in PTC patients, we investigated the role of those genes in PTC. We divided PTC patients into two clusters (Figure 3G), and details of the consensus cluster analysis are provided in Supplementary Figure 1.

K-M curve shows that cluster 2 had a worse OS (p = 0.021) (Figure 3H), and the clinicopathological heatmap indicated that the number of PTC with pathological type, RAS, BRAF, T-stage, N-stage, M-stage were statistically different between clusters (Figure 3I). These results suggest that these 16 prognostic cellular senescence genes can be used to predict OS in PTC patients and the expression of those genes can provide a new way to differentiate PTC subtypes. The results of the CIBERSORT analysis showed that the content of ten immune cells (T cells CD8, T cells follicular helper, T cells regulatory, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Eosinophils) were statistically different between the two clusters (Figure 3J). Furthermore, the ESTIMATE analysis showed lower ESTIMATEScore, ImmuneScore, StromalScore and higher TumorPurity in cluster 2 (Figure 3K–3N). The results of CIBERSORT and ESTIMATE analysis indicated that cellular senescence genes may affect the prognosis of PTC patients through the immune microenvironment.

To further analyze the reasons for the differing prognosis of the two clusters, we performed KEGG pathway enrichment analysis for clusters 1 and 2, respectively. The results showed that cluster 1 was mainly enriched in allograft rejection, asthma, cell adhesion molecules (CAMs), intestinal immune network for IgA production, leishmania infection, leukocyte transendothelial migration, natural killer cell mediated cytotoxicity, pathogenic Escherichia coli infection, type I diabetes mellitus and viral myocarditis. The results showed that cluster 2 was mainly enriched in arginine and proline metabolism, ascorbate and aldarate metabolism, β alanine metabolism, butanoate metabolism, fatty acid metabolism, glycerolipid metabolism, glycosylphosphatidylinositol (GPI)-anchor biosynthesis, lysine degradation, peroxisome, valine leucine and isoleucine degradation. The results of the analysis in cluster 2 suggest that cellular senescence genes may contribute to the worsening of the prognosis of PTC patients through metabolism-related pathways.

Construction and validation of CSRS

In the training cohort, we constructed CSRS (Supplementary Figure 2) by selecting four genes from 16 prognosis-related cellular senescence genes by LASSO regression analysis for prediction OS. We got the following risk score formula for patients with PTC: Risk score = (0.05013 × expression value of SNAI1) + (0.097623 × expression value of CDKN2A) + (0.54932 × expression value of HDAC4) + (0.02454 × expression value of NDRG1). Based on median the risk score, we divided the patients in the three cohorts into high and low-risk groups.

The K-M survival curve analysis showed that patients with high-risk groups had shorter OS (Figure 4A–4C). The AUC value of PTC patients in 1-10 years was greater than 0.737 that the model has good accuracy (Figure 4D–4F). To compare the predictive power of CSRS and clinicopathological characteristics for OS, we evaluated the AUC value of CSRS and clinical pathological characteristics in the fifth year. The results demonstrated that the risk scores achieved the highest AUC value than other clinical pathological characteristics (Figure 4G–4I). These results indicated that CSRS has good predictive ability and clinical application value. Survival status, a heatmap of 4-gene expression levels, distribution of risk scores in high and low risk groups (Figure 4J–4O and Supplementary Figure 3). We analyzed gene expression data from the GEPIA2 website and five GEO cohorts, and the results further supported the previous conclusion (Figure 5) (the SNAI1 gene expression level showed a trend of lower expression in tumor tissues in GSE3678 and GSE27155, Figure 5A–5X).

Figure 4. Construction and validation of prognostic models. K-M curve of OS probability between high and low risk groups in training cohort (A), validation cohort (B) and all TCGA-THCA (C). The ROC curves of prognostic signature of prognostic model in 1-10 years in training cohort (D), in validation cohort (E) and in all TCGA-THCA cohort (F). The ROC curves of risk score, age, gender, stage, M, N, T, BRAF, and RAS in 5 years in training cohort (G), in validation cohort (H) and in all TCGA-THCA cohort (I). Survival status of patients in the training cohort (J), validation cohort (K) and all TCGA-THCA (L). Expression heatmap of prognostic model genes in the training cohort (M), validation cohort (N) and all TCGA-THCA (O).

Figure 5. Validation of gene expression levels in prognostic models. Expression of SNAI1 in GTEx (A), in GSE3467 (B), in GSE3678 (C), in GSE5364 (D), in GSE27155 (E), in GSE35157 (F). Expression of CDKN2A in GTEx (G), in GSE3467 (H), in GSE3678 (I), in GSE5364 (J), in GSE27155 (K), in GSE35157 (L). Expression of HDAC4 in GTEx (M), in GSE3467 (N), in GSE3678 (O), in GSE5364 (P), in GSE27155 (Q), in GSE35157 (R). Expression of NDRG1 in GTEx (S), in GSE3467 (T), in GSE3678 (U), in GSE5364 (V), in GSE27155 (W), in GSE35157 (X).

The risk score in hierarchical analysis is the predictor for PTC

To further mitigate the effect of clinicopathological features on risk scores predicting OS, we conducted a stratified analysis based on the clinicopathological features of the entire TCGA-THCA cohort. The K-M survival curve results in the stratified analysis showed that the high-risk group had a shorter OS, except for patients with age < 55, M1, RAS-Mutation, Follicular Cell Type, Radiation therapy−NO (Figure 6).

Figure 6. Hierarchical analysis. K-M curve analysis of OS probability was based on risk scores grouped by age (A, B), sex (C, D), stage (E, F), T (G, H), M (I, J), N (K, L), BRAF (M, N), RAS (O, P), pathological subtype (Q, R) and radiation therapy (S, T).

Multiple omics differences between high - and low-risk groups

A heatmap of immune cells with statistically significant differences between high and low risk groups based on the analysis algorithms of immune cell contents (TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, EPIC) was presented in Figure 7A. These results provide a basis for further elucidating the relationship between cell senescence and immune cells. To investigate the biological behavior of the high-risk group, we analyzed by GSVA and found that the high-risk group was mainly enriched pathway in adipocytokine signaling pathway, lysine degradation, glycosaminoglycan biosynthesis-heparan sulfate, glycosaminoglycan biosynthesis-chondroitin sulfate, acute myeloid leukemia, gap junction, axon guidance, focal adhesion, TGF-β signaling pathway, dorso-ventral axis formation, hypertrophic cardiomyopathy (HCM), ECM-receptor interaction, hedgehog signaling pathway, basal cell carcinoma (Figure 7B). Furthermore, the heat map of the correlation between clinicopathological characteristics and risk scores showed that there were statistical differences in radiation therapy, RAS, Age and Cluster in the high and low risk groups (Figure 7C). Additionally, four IPS scores were statistically different between high and low risk groups (Figure 7D–7G). The data provided by the TIDE website showed that TIDE scores, Dysfunction scores, Exclusion scores and CAF scores were statistically different between high and low risk groups (Figure 7H–7K).

Figure 7. Multi-omics analysis. Analysis results of six immune cell infiltration algorithms in high and low risk groups (A). KEGG pathway enrichment analysis in high and low risk groups (B). Heatmap of clinicopathological features in high and low risk groups (C). Comparison of immunophenoscores (IPS) between high and low risk groups who have not been treated with anti-CTLA4 or anti-PD-1 immunotherapy (D), who have been treated with anti-CTLA4 and anti-PD-1 immunotherapy (E), who have been treated with anti-CTLA4 and anti-PD-1 immunotherapy, who have only been treated with anti-PD-1 immunotherapy (F), and who have only been treated with anti-CTLA4 immunotherapy (G). Comparison of TIDE (H), Dysfunction (I), Exclusion (J), CAF (K) between high and low risk groups.

Through an analyzing the expression levels of immune checkpoints, we found that 13 immune checkpoints (ADORA2A, TNFSF9, CD274, NRP1, IDO2, CD40, CD28, IDO1, CD160, TNFSF4, CD80, TNFSF15, CD44) were differentially expressed in high and low risk groups (Figure 7A). Subsequently, we identified ten drugs (BX.795, PF.562271, PF.02341066, Pazopanib, PAC.1, MK.2206, IPA.3, Imatinib, GDC0941, Embelin) that were more sensitive in high-risk patients by drug sensitivity analysis (Figure 8B–8K).

Figure 8. Analysis of drug treatment. There are differences in immune checkpoint expression between groups at high and low risk (A). Ten drugs with lower half-maximum inhibitory concentration in the high-risk group (B–K).

Author Contributions

Conceptualization, J.Z.; methodology, Y.H.; software, H.J.; validation, G.X., X.L. and W.C.; formal analysis, G.X.; investigation, X.L.; resources, W.C.; data curation, G.X.; writing—original draft preparation, Y.H.; writing—review and editing, H.J.; visualization, Y.H.; supervision, Y.L.; project administration, H.J.; funding acquisition, Y.L. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

Thanks to Miss Nan Xiong for editing and proofreading the manuscript.

Conflicts of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This research was funded by the National Natural Science Foundation of China (81600602).

References

• 1. Xia C, Dong X, Li H, Cao M, Sun D, He S, Yang F, Yan X, Zhang S, Li N, Chen W. Cancer statistics in China and United States, 2022: profiles, trends, and determinants. Chin Med J (Engl). 2022; 135:584–90. https://doi.org/10.1097/CM9.0000000000002108 [PubMed]
• 2. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022; 72:7–33. https://doi.org/10.3322/caac.21708 [PubMed]
• 3. Udelsman R, Zhang Y. The epidemic of thyroid cancer in the United States: the role of endocrinologists and ultrasounds. Thyroid. 2014; 24:472–9. https://doi.org/10.1089/thy.2013.0257 [PubMed]
• 4. Ahn HS, Welch HG. South Korea’s Thyroid-Cancer “Epidemic”--Turning the Tide. N Engl J Med. 2015; 373:2389–90. https://doi.org/10.1056/NEJMc1507622 [PubMed]
• 5. Brito JP, Morris JC, Montori VM. Thyroid cancer: zealous imaging has increased detection and treatment of low risk tumours. BMJ. 2013; 347:f4706. https://doi.org/10.1136/bmj.f4706 [PubMed]
• 6. Caulley L, Eskander A, Yang W, Auh E, Cairncross L, Cho NL, Golbon B, Iyer S, Liu JC, Lee PJ, Lindeman B, Meltzer C, Molin N, et al. Trends in Diagnosis of Noninvasive Follicular Thyroid Neoplasm With Papillarylike Nuclear Features and Total Thyroidectomies for Patients With Papillary Thyroid Neoplasms. JAMA Otolaryngol Head Neck Surg. 2022; 148:99–106. https://doi.org/10.1001/jamaoto.2021.3277 [PubMed]
• 7. Kim SY, Kim HJ, Kim SM, Chang H, Lee YS, Chang HS, Park CS. Thyroid Hormone Supplementation Therapy for Differentiated Thyroid Cancer After Lobectomy: 5 Years of Follow-Up. Front Endocrinol (Lausanne). 2020; 11:520. https://doi.org/10.3389/fendo.2020.00520 [PubMed]
• 8. Wilson M, Patel A, Goldner W, Baker J, Sayed Z, Fingeret AL. Postoperative thyroid hormone supplementation rates following thyroid lobectomy. Am J Surg. 2020; 220:1169–73. https://doi.org/10.1016/j.amjsurg.2020.06.052 [PubMed]
• 9. Ullmann TM, Papaleontiou M, Sosa JA. Current Controversies in Low-Risk Differentiated Thyroid Cancer: Reducing Overtreatment in an Era of Overdiagnosis. J Clin Endocrinol Metab. 2023; 108:271–80. https://doi.org/10.1210/clinem/dgac646 [PubMed]
• 10. Perrier ND, Brierley JD, Tuttle RM. Differentiated and anaplastic thyroid carcinoma: Major changes in the American Joint Committee on Cancer eighth edition cancer staging manual. CA Cancer J Clin. 2018; 68:55–63. https://doi.org/10.3322/caac.21439 [PubMed]
• 11. Cui Y, Mubarik S, Li R, Na, Yu C. Trend dynamics of thyroid cancer incidence among China and the U.S. adult population from 1990 to 2017: a joinpoint and age-period-cohort analysis. BMC Public Health. 2021; 21:624. https://doi.org/10.1186/s12889-021-10635-w [PubMed]
• 12. Shah Y, Verma A, Marderstein AR, White J, Bhinder B, Garcia Medina JS, Elemento O. Pan-cancer analysis reveals molecular patterns associated with age. Cell Rep. 2021; 37:110100. https://doi.org/10.1016/j.celrep.2021.110100 [PubMed]
• 13. López-Otín C, Blasco MA, Partridge L, Serrano M, Kroemer G. The hallmarks of aging. Cell. 2013; 153:1194–217. https://doi.org/10.1016/j.cell.2013.05.039 [PubMed]
• 14. Xu M, Pirtskhalava T, Farr JN, Weigand BM, Palmer AK, Weivoda MM, Inman CL, Ogrodnik MB, Hachfeld CM, Fraser DG, Onken JL, Johnson KO, Verzosa GC, et al. Senolytics improve physical function and increase lifespan in old age. Nat Med. 2018; 24:1246–56. https://doi.org/10.1038/s41591-018-0092-9 [PubMed]
• 15. van Deursen JM. Senolytic therapies for healthy longevity. Science. 2019; 364:636–7. https://doi.org/10.1126/science.aaw1299 [PubMed]
• 16. Di Micco R, Krizhanovsky V, Baker D, d’Adda di Fagagna F. Cellular senescence in ageing: from mechanisms to therapeutic opportunities. Nat Rev Mol Cell Biol. 2021; 22:75–95. https://doi.org/10.1038/s41580-020-00314-w [PubMed]
• 17. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov. 2022; 12:31–46. https://doi.org/10.1158/2159-8290.CD-21-1059 [PubMed]
• 18. Romaniello D, Gelfo V, Pagano F, Sgarzi M, Morselli A, Girone C, Filippini DM, D’Uva G, Lauriola M. IL-1 and senescence: Friends and foe of EGFR neutralization and immunotherapy. Front Cell Dev Biol. 2023; 10:1083743. https://doi.org/10.3389/fcell.2022.1083743 [PubMed]
• 19. Bleyer A, Barr R, Hayes-Lattin B, Thomas D, Ellis C, Anderson B, and Biology and Clinical Trials Subgroups of the US National Cancer Institute Progress Review Group in Adolescent and Young Adult Oncology. The distinctive biology of cancer in adolescents and young adults. Nat Rev Cancer. 2008; 8:288–98. https://doi.org/10.1038/nrc2349 [PubMed]
• 20. Keegan TH, Ries LA, Barr RD, Geiger AM, Dahlke DV, Pollock BH, Bleyer WA, and National Cancer Institute Next Steps for Adolescent and Young Adult Oncology Epidemiology Working Group. Comparison of cancer survival trends in the United States of adolescents and young adults with those in children and older adults. Cancer. 2016; 122:1009–16. https://doi.org/10.1002/cncr.29869 [PubMed]
• 21. Avelar RA, Ortega JG, Tacutu R, Tyler EJ, Bennett D, Binetti P, Budovsky A, Chatsirisupachai K, Johnson E, Murray A, Shields S, Tejada-Martinez D, Thornton D, et al. A multidimensional systems biology analysis of cellular senescence in aging and disease. Genome Biol. 2020; 21:91. https://doi.org/10.1186/s13059-020-01990-9 [PubMed]
• 22. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
• 23. Chen T, Zhang H, Liu Y, Liu YX, Huang L. EVenn: Easy to create repeatable and editable Venn diagrams and Venn networks online. J Genet Genomics. 2021; 48:863–6. https://doi.org/10.1016/j.jgg.2021.07.007 [PubMed]
• 24. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–7. https://doi.org/10.1038/nmeth.3337 [PubMed]
• 25. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
• 26. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
• 27. Huang R, Li G, Wang Z, Hu H, Zeng F, Zhang K, Wang K, Wu F. Identification of an ATP metabolism-related signature associated with prognosis and immune microenvironment in gliomas. Cancer Sci. 2020; 111:2325–35. https://doi.org/10.1111/cas.14484 [PubMed]
• 28. Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med. 1997; 16:385–95. https://doi.org/10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3 [PubMed]
• 29. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010; 33:1–22. https://doi.org/10.18637/jss.v033.i01 [PubMed]
• 30. Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019; 47:W556–60. https://doi.org/10.1093/nar/gkz430 [PubMed]
• 31. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017; 77:e108–10. https://doi.org/10.1158/0008-5472.CAN-17-0307 [PubMed]
• 32. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, Krogsdam A, Loncova Z, Posch W, Wilflingseder D, Sopper S, Ijsselsteijn M, Brouwer TP, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med. 2019; 11:34. https://doi.org/10.1186/s13073-019-0638-6 [PubMed]
• 33. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, de Reyniès A. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016; 17:218. https://doi.org/10.1186/s13059-016-1070-5 [PubMed]
• 34. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017; 18:220. https://doi.org/10.1186/s13059-017-1349-1 [PubMed]
• 35. Racle J, Gfeller D. EPIC: A Tool to Estimate the Proportions of Different Cell Types from Bulk Gene Expression Data. Methods Mol Biol. 2020; 2120:233–48. https://doi.org/10.1007/978-1-0716-0327-7_17 [PubMed]
• 36. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
• 37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]
• 38. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–8. https://doi.org/10.1038/s41591-018-0136-1 [PubMed]
• 39. Chen Y, Li ZY, Zhou GQ, Sun Y. An Immune-Related Gene Prognostic Index for Head and Neck Squamous Cell Carcinoma. Clin Cancer Res. 2021; 27:330–41. https://doi.org/10.1158/1078-0432.CCR-20-2166 [PubMed]
• 40. Wang J, Ma X, Ma Z, Ma Y, Wang J, Cao B. Research Progress of Biomarkers for Immune Checkpoint Inhibitors on Digestive System Cancers. Front Immunol. 2022; 13:810539. https://doi.org/10.3389/fimmu.2022.810539 [PubMed]
• 41. Xu K, He J, Zhang J, Liu T, Yang F, Ren T. A novel prognostic risk score model based on immune-related genes in patients with stage IV colorectal cancer. Biosci Rep. 2020; 40:BSR20201725. https://doi.org/10.1042/BSR20201725 [PubMed]
• 42. Xu Y, Chen Y, Niu Z, Xing J, Yang Z, Yin X, Guo L, Zhang Q, Qiu H, Han Y. A Novel Pyroptotic and Inflammatory Gene Signature Predicts the Prognosis of Cutaneous Melanoma and the Effect of Anticancer Therapies. Front Med (Lausanne). 2022; 9:841568. https://doi.org/10.3389/fmed.2022.841568 [PubMed]
• 43. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014; 9:e107468. https://doi.org/10.1371/journal.pone.0107468 [PubMed]
• 44. Chen W, Zheng R, Zuo T, Zeng H, Zhang S, He J. National cancer incidence and mortality in China, 2012. Chin J Cancer Res. 2016; 28:1–11. https://doi.org/10.3978/j.issn.1000-9604.2016.02.08 [PubMed]
• 45. Gilliland FD, Hunt WC, Morris DM, Key CR. Prognostic factors for thyroid carcinoma. A population-based study of 15,698 cases from the Surveillance, Epidemiology and End Results (SEER) program 1973-1991. Cancer. 1997; 79:564–73. https://doi.org/10.1002/(sici)1097-0142(19970201)79:3<564::aid-cncr20>3.0.co;2-0 [PubMed]
• 46. Lee YH, Lee YM, Sung TY, Yoon JH, Song DE, Kim TY, Baek JH, Ryu JS, Chung KW, Hong SJ. Is Male Gender a Prognostic Factor for Papillary Thyroid Microcarcinoma? Ann Surg Oncol. 2017; 24:1958–64. https://doi.org/10.1245/s10434-017-5788-4 [PubMed]
• 47. Erbe R, Wang Z, Wu S, Xiu J, Zaidi N, La J, Tuck D, Fillmore N, Giraldo NA, Topper M, Baylin S, Lippman M, Isaacs C, et al. Evaluating the impact of age on immune checkpoint therapy biomarkers. Cell Rep. 2021; 36:109599. https://doi.org/10.1016/j.celrep.2021.109599 [PubMed]
• 48. Chatsirisupachai K, Lesluyes T, Paraoan L, Van Loo P, de Magalhães JP. An integrative analysis of the age-associated multi-omic landscape across cancers. Nat Commun. 2021; 12:2345. https://doi.org/10.1038/s41467-021-22560-y [PubMed]
• 49. Chatsirisupachai K, Palmer D, Ferreira S, de Magalhães JP. A human tissue-specific transcriptomic analysis reveals a complex relationship between aging, cancer, and cellular senescence. Aging Cell. 2019; 18:e13041. https://doi.org/10.1111/acel.13041 [PubMed]
• 50. Kim YH, Choi YW, Han JH, Lee J, Soh EY, Park SH, Kim JH, Park TJ. TSH signaling overcomes B-RafV600E-induced senescence in papillary thyroid carcinogenesis through regulation of DUSP6. Neoplasia. 2014; 16:1107–20. https://doi.org/10.1016/j.neo.2014.10.005 [PubMed]
• 51. Kim YH, Choi YW, Lee J, Soh EY, Kim JH, Park TJ. Senescent tumor cells lead the collective invasion in thyroid cancer. Nat Commun. 2017; 8:15208. https://doi.org/10.1038/ncomms15208 [PubMed]
• 52. Lin W, Wang X, Wang Z, Shao F, Yang Y, Cao Z, Feng X, Gao Y, He J. Comprehensive Analysis Uncovers Prognostic and Immunogenic Characteristics of Cellular Senescence for Lung Adenocarcinoma. Front Cell Dev Biol. 2021; 9:780461. https://doi.org/10.3389/fcell.2021.780461 [PubMed]
• 53. Lin W, Wang X, Xu Z, Wang Z, Liu T, Cao Z, Feng X, Gao Y, He J. Identification and validation of cellular senescence patterns to predict clinical outcomes and immunotherapeutic responses in lung adenocarcinoma. Cancer Cell Int. 2021; 21:652. https://doi.org/10.1186/s12935-021-02358-0 [PubMed]
• 54. Guan L, Crasta KC, Maier AB. Assessment of cell cycle regulators in human peripheral blood cells as markers of cellular senescence. Ageing Res Rev. 2022; 78:101634. https://doi.org/10.1016/j.arr.2022.101634 [PubMed]
• 55. Furuya S, Endo K, Takahashi A, Miyazawa K, Saitoh M. Snail suppresses cellular senescence and promotes fibroblast-led cancer cell invasion. FEBS Open Bio. 2017; 7:1586–97. https://doi.org/10.1002/2211-5463.12300 [PubMed]
• 56. Xu S, Zhang L, Cheng X, Yu H, Bao J, Lu R. Capsaicin inhibits the metastasis of human papillary thyroid carcinoma BCPAP cells through the modulation of the TRPV1 channel. Food Funct. 2018; 9:344–54. https://doi.org/10.1039/c7fo01295k [PubMed]
• 57. Vosgha H, Ariana A, Smith RA, Lam AK. miR-205 targets angiogenesis and EMT concurrently in anaplastic thyroid carcinoma. Endocr Relat Cancer. 2018; 25:323–37. https://doi.org/10.1530/ERC-17-0497 [PubMed]
• 58. Ma S, Jia W, Ni S. miR-199a-5p inhibits the progression of papillary thyroid carcinoma by targeting SNAI1. Biochem Biophys Res Commun. 2018; 497:181–6. https://doi.org/10.1016/j.bbrc.2018.02.051 [PubMed]
• 59. Emadi Baygi M, Soheili ZS, Schmitz I, Sameie S, Schulz WA. Snail regulates cell survival and inhibits cellular senescence in human metastatic prostate cancer cell lines. Cell Biol Toxicol. 2010; 26:553–67. https://doi.org/10.1007/s10565-010-9163-5 [PubMed]
• 60. Ferru A, Fromont G, Gibelin H, Guilhot J, Savagner F, Tourani JM, Kraimps JL, Larsen CJ, Karayan-Tapon L. The status of CDKN2A alpha (p16INK4A) and beta (p14ARF) transcripts in thyroid tumour progression. Br J Cancer. 2006; 95:1670–7. https://doi.org/10.1038/sj.bjc.6603479 [PubMed]
• 61. Wang HQ, Li Y, Song X, Ma YQ, Li JL, Li YX, Wang GF, Liu P, Liu PL, Cao S, Shi HY. Significance of interstitial fibrosis and p16 in papillary thyroid carcinoma. Endocr J. 2022; 69:1253–9. https://doi.org/10.1507/endocrj.EJ22-0010 [PubMed]
• 62. Jiang JL, Tian GL, Chen SJ, Xu LI, Wang HQ. Promoter methylation of p16 and RASSF1A genes may contribute to the risk of papillary thyroid cancer: A meta-analysis. Exp Ther Med. 2015; 10:1549–55. https://doi.org/10.3892/etm.2015.2656 [PubMed]
• 63. Zou M, Baitei EY, Al-Rijjal RA, Parhar RS, Al-Mohanna FA, Kimura S, Pritchard C, Binessa HA, Alzahrani AS, Al-Khalaf HH, Hawwari A, Akhtar M, Assiri AM, et al. TSH overcomes Braf(V600E)-induced senescence to promote tumor progression via downregulation of p53 expression in papillary thyroid cancer. Oncogene. 2016; 35:1909–18. https://doi.org/10.1038/onc.2015.253 [PubMed]
• 64. Hou F, Wei W, Qin X, Liang J, Han S, Han A, Kong Q. The posttranslational modification of HDAC4 in cell biology: Mechanisms and potential targets. J Cell Biochem. 2020; 121:930–7. https://doi.org/10.1002/jcb.29365 [PubMed]
• 65. Di Giorgio E, Paluvai H, Dalla E, Ranzino L, Renzini A, Moresi V, Minisini M, Picco R, Brancolini C. HDAC4 degradation during senescence unleashes an epigenetic program driven by AP-1/p300 at selected enhancers and super-enhancers. Genome Biol. 2021; 22:129. https://doi.org/10.1186/s13059-021-02340-z [PubMed]
• 66. Feng Y, Wang X, Xu L, Pan H, Zhu S, Liang Q, Huang B, Lu J. The transcription factor ZBP-89 suppresses p16 expression through a histone modification mechanism to affect cell senescence. FEBS J. 2009; 276:4197–206. https://doi.org/10.1111/j.1742-4658.2009.07128.x [PubMed]
• 67. Lu WJ, Chua MS, So SK. Suppressing N-Myc downstream regulated gene 1 reactivates senescence signaling and inhibits tumor growth in hepatocellular carcinoma. Carcinogenesis. 2014; 35:915–22. https://doi.org/10.1093/carcin/bgt401 [PubMed]
• 68. Xie Z, Li X, He Y, Wu S, Wang S, Sun J, He Y, Lun Y, Zhang J. Immune Cell Confrontation in the Papillary Thyroid Carcinoma Microenvironment. Front Endocrinol (Lausanne). 2020; 11:570604. https://doi.org/10.3389/fendo.2020.570604 [PubMed]