Prognostic Protein Biomarker Screening for Thyroid Carcinoma Based on Cancer Proteomics Profiles

Article type: Research Article

Keywords: thyroid carcinoma, TCPA, overall survival, TMB

Authors: Pu Xie ⚬, Qinglei Yin, Shu Wang ⚬, Dalong Song

Affiliations: Department of Endocrine and Metabolic Diseases, Shanghai Institute of Endocrine and Metabolic Diseases, Ruijin Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai 200025, China; xiepu1997@163.com; Shanghai National Clinical Research Center for Metabolic Diseases, Key Laboratory for Endocrine and Metabolic Diseases of the National Health Commission of the PR China, Shanghai Key Laboratory for Endocrine Tumor, Ruijin Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai 200025, China; Guangdong Geriatric Institute, Guangdong Provincial People’s Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou 510000, China; yinqinglei128@163.com; Reproductive Medicine Center, Department of Obstetrics and Gynecology, Guangdong Provincial People’s Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou 510000, China

License: © 2024 by the authors. CC BY 4.0 Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Article links: DOI: 10.3390/biomedicines12092066  |  PubMed: 39335579  |  PMC: PMC11428938

Relevance: Relevant: mentioned in keywords or abstract

Full text: PDF (4.4 MB)

1. Introduction

Thyroid carcinoma (THCA) is one of the most common neoplastic diseases [ref. 1] and is expected to become the fourth leading type of cancer in the world [ref. 2]. In the USA, the incidence of THCA has increased rapidly from 4.56 cases per 100,000 person-years to 14.41 cases in the past four decades [ref. 3]. Risk factors of THCA include radiation exposure at a young age, excess body weight, hormonal exposures, certain environmental pollutants, and family history [ref. 4]. Most patients with THCA can be treated by surgery or radioactive iodine therapy followed by endocrine suppression [ref. 5]. Nevertheless, the mortality rate of advanced thyroid cancer remains a significant concern. Reports indicate an increased incidence of more advanced forms and higher morbidity associated with thyroid carcinoma [ref. 6,ref. 7]. Thus, the identification of novel biomarkers for prognosis and drug screening in THCA, providing additional strategies for survival and optimizing treatments, is of great significance.

The pathogenesis of THCA is not fully understood. With the development of high-throughput technologies, recent decades have witnessed great improvements in cancer genomics [ref. 8]. Using a multi-omics approach in recent years, the characteristics of THCA have been extensively studied in the aspects of genomics, transcriptomics, and proteomics [ref. 9]. With the application of proteomics, mass spectrometry-based proteomics enable the reclassification of cancer, which may affect prognosis and guide clinical decision-making [ref. 10]. Recently, prognostic signatures from public databases have garnered significant attention and demonstrated potentials for predicting outcomes in THCA patients. However, these studies were mainly focused on transcriptomic profiles with little known about proteome modulation and protein function. Previous studies have shown that, in many cases, transcript levels are insufficient to predict protein levels [ref. 11] and the average correlation between mRNA expression and protein expression is below 0.5 [ref. 12,ref. 13]. Therefore, prognostic signatures based on proteomic level will be crucial.

The application of advanced proteomics technologies and bioinformatics approaches to THCA research could offer unprecedented insight into cancer biology and treatment [ref. 14]. In our study, we determined a protein-associated prognostic signature from The Cancer Proteome Atlas (TCPA) THCA cohort. A risk model was thus established for predicting THCA patients’ prognosis. To validate the prognostic signature, we investigated its efficiency and accuracy using both training and testing sets. We also explored the association of this signature with immune cell infiltration and the tumor mutation burden (TMB) based on results from The Cancer Genome Atlas (TCGA). Importantly, our study proved a new approach for predicting the response to treatment, including chemotherapy and targeted therapy in THCA patients.

2. Materials and Methods

2.1. Data Collection

The protein expression data of THCA were retrieved from the TCPA (https://www.tcpaportal.org/tcpa/ (accessed on 18 July 2022)). The RNA-seq expression profile data (HTSeq-FPKM) of THCA patients were obtained from the TCGA data portal (http://portal.gdc.cancer.gov/cart (accessed on 18 July 2022)). Clinical information such as age, gender, survival time, survival status, and clinical stage was also downloaded from the TCGA data portal. Patients were randomly split into a training set and a testing set at a ratio of 5:5 using the “caret” packages. The distribution of age, gender, and clinical stage was similar between the two datasets (Supplementary Table S1). The training set was used to identify the prognostic signature and the testing set was used to validate its prognostic capability.

2.2. Construction of the Protein Prognostic Signature

To explore the possible protein in relation to prognosis for THCA patients, univariate Cox proportional hazard regression analysis was used to identify the relationship between prognostic proteins and OS in the training set, and p < 0.05 was considered to be statistically significant. After that, we used the “glmnet” packages [ref. 15,ref. 16], prognostic and selection operation (LASSO) to minimize overfitting and identify the most significant prognosis-related proteins. p < 0.05 was considered statistically significant and p < 0.01 was considered a significant difference. Risk score calculation formula: expression protein 1 × coefficient protein 1 + expression protein 2× coefficient protein 2 + … + expression protein n × coefficient protein n. A Kaplan–Meier plot was performed to compare the overall survival between two risk groups using the R package survival. The time-dependent receiver operating characteristic (ROC) curve was used to assess all the predictive values of the prognostic protein signature for overall survival using R package survivalROC. The AUC (Area Under ROC Curve) value of the ROC indicates that 0.5–0.7 is acceptable, 0.7–0.9 is good, and >0.9 is excellent.

2.3. Nomogram Construction and Calibration

The risk score and relevant clinical parameters such as age and stage were incorporated into the construction of a prognostic nomogram via an “rms” R package (version 6.2.0) to predict 1-, 3-, and 5-year OS of THCA patients in the TCPA cohort. We used a calibration plot comparing predicted and observed overall survival to evaluate the performance of the prognostic nomogram (method = “boot”, B = 1000).

2.4. DEPs Identification

We conducted a differential protein expression analysis in proteomic data. The R package limma (version 4.1.3) and Wilcox test were used to identify differentially expressed proteins (DEPs) between the high-risk and the low-risk group. Proteins with p < 0.05 were considered as DEPs.

2.5. Principal Component Analysis (PCA), GO and KEGG Analysis

PCA analysis was performed using the limma and scatterplot3d packages to explore the distribution of patients with different risk scores. To assess the potential biologic functions of differentially expressed proteins, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed by the cluster Profiler package in R. Functional categories with an adjusted p value < 0.05 were considered as significant pathways. To investigate the functional biological differences of the hub proteins, GO analysis was automatically completed and visualized by the Metascape [ref. 17].

2.6. Protein–Protein Interaction Network Construction

The common significant DEPs were mapped to the STRING database (http://www.string-db.org/ (accessed on 20 July 2022)), then set confidence > 0.7 to conduct the PPI network analysis, which were visualized by Cytoscape software (v3.9.1) [ref. 18]. Module analysis was performed to identify significant modules and explore the hub proteins using the Molecular Complex Detection (MCODE v2.0.0) [ref. 19], which is a plug-in of Cytoscape (MCODE score ≥ 6).

2.7. Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was utilized to elucidate the molecular mechanisms. We divided the samples in the entire set into the high-risk and the low-risk groups based on scores from the signature and compared the KEGG pathway enrichment between the two groups. A p-value of <0.05 was considered significant.

2.8. Analysis of Tumor Infiltrating Immune Cells

We employed the CIBERSORT to evaluate the immune cell infiltration in each sample based on the RNA-seq data. With “CIBERSORT” (R package), we used the CIBERSORT algorithm to analyze gene expression data. Using the standard p < 0.05, we screened samples and calculated the percentage of 22 immune cells.

2.9. Analysis of Tumor Mutation Burden

The mutation data collected from TCGA were analyzed by R package maftools. The tumor mutation burden (TMB) was calculated using the formula: TMB = (total mutation/total covered bases) × 10⁶.

2.10. Immunotherapy Analysis

The Tumor Immune Dysfunction and Exclusion (TIDE) is a data-driven Web platform that combines large-scale omics data from 33,000 cases across 188 cohorts, 998 tumor samples from 12 immune checkpoint blockade (ICB) clinical studies, and eight clustered regularly interspaced short palindromic repeat (CRISPR) screens. TIDE aids in hypothesis generation and the optimization of immunological biomarkers. In this study, we used TIDE to evaluate immune-suppressive rejection scores and therapeutic responses to ICB.

2.11. Potential Chemotherapeutic Response

The response to common chemo drugs, such as temsirolimus, roscovitine, AZD6244, AKT inhibitor, and metformin, were predicted by the Genomics of Drug Sensitivity in Cancer (GDSC) (http://www.cancerrxgene.org/ (accessed on 20 July 2022)) to analyze the relationship between the signature and chemotherapeutic response. We used the R package pRRophetic to estimate and compare the half-maximal inhibitory concentration (IC50) between different risk groups.

2.12. Statistical Analysis

The differences in immune cells between the risk score were compared with the Wilcox test. The correlations between the expression of every protein in the signature and the expression of immune checkpoints were also calculated with Spearman rank correlation. All statistical analyses were performed with R software (version 4.1.3), and p < 0.05 was selected as statistically significant.

3. Results

3.1. Establishment of a Proteomic Prognostic Signature in Thyroid Carcinoma

The whole dataset was split into a training set (n = 188) and a testing set (n = 188). The training set was used to develop the protein prognostic model, while the testing set and the entire dataset were used to validate the prognostic model. To screen prognosis-related proteins, univariate Cox regression analysis was conducted to identify proteins associated with prognosis. After univariate Cox analysis, 28 proteins met the criteria of p < 0.05 and were retained for further analysis. Among them, 18 proteins were associated with increased risk (HRs > 1), while the other 10 proteins were protective proteins (HRs < 1) (Figure 1A). Next, least absolute shrinkage and selection operator (LASSO) penalized Cox regression were applied to reduce the prognostic proteinsand a 13-protein signature was constructed according to the optimum λ value (Figure 1B,C, Supplementary Figure S1). The risk score could be calculated using the data in Table 1. Principal component analysis (PCA) demonstrated that the prognostic signature could effectively distinguish THCA patients into two groups (Supplementary Figure S2). We evaluated the risk model’s survival prediction performance using validation data. Samples from the training, testing, and entire datasets were categorized into high- and low-risk groups based on the median risk score. The individual risk score and survival status and the key protein expression profiles of patients are shown in Figure 1. Kaplan–Meier survival curves showed significantly worse overall survival (OS) for the high-risk group compared to the low-risk group (Figure 1D–F). Survival ROC curve analysis was also performed, and the area under the curve (AUC) at 1, 3, and 5 years was 1.000, 0.940, and 0.932, respectively, which verified the predictive value of the signature, and the AUC value of the test and the whole set also indicated a similar potential for predicting survival (Figure 1G–I). According to the risk graph, patient mortality appears to increase as the risk score increases (Figure 1J–L). In the high-risk group, S100A4, PAI1, IGFBP2, RICTOR, B7-H3, COLLAGENVI, PAR, SNAIL, FAK, and Connexin-43 proteins were upregulated, while Rheb, EVI1, and P90RSK_pT359S363 proteins were downregulated (Figure 1J–L). The risk curves and protein profile heatmap displayed features consistent with both the test set and the whole set (Figure 1J–L). These findings indicate that the protein profiles may serve as effective prognostic biomarkers for THCA patients.

3.2. Construction of Nomogram Based on the Protein Signature and Clinical Data

To assess the independent prognostic force of the prognostic signature, both univariable and multivariable Cox proportional hazard regression models were applied in the training, the testing, and the entire sets (Supplementary Figure S3). The results demonstrated that the prognostic signature could serve as an independent predictor. A nomogram was then created using age, gender, TMN status, clinical stage, and risk-score levels to predict 1-, 3-, and 5-year overall survival in THCA patients (Figure 2A). The calibration curve presented desirable prediction of the nomogram for the 1-, 3-, and 5-year clinical outcomes (Figure 2B–D). Furthermore, it is notable that the AUC of the nomogram exceeded that of other characteristics (Figure 2E–G). These results suggest that the predictive nomogram was highly accurate for forecasting survival in THCA patients.

3.3. Clinical Relevance Assessment and Construction of the Protein Coexpression Network

In order to explore whether the prognostic signature participated in the development and progression of THCA, risk scores with clinical characteristics were compared. There were significant differences in the risk scores between different groups in age (p < 0.05) and TNM stage (p = 0.002) (Figure 3A,B), but no significant differences in gender, T stage, N stage, and M stage (Figure 3C–F). In addition, stratification analysis was further conducted to investigate the prognostic significance of THCA patients in subgroups. Our analyses suggested that a protein-based signature showed excellent performance in predicting outcome in age > 65 (p < 0.05), male (p < 0.05), female (p < 0.01), T3–4 stage (p < 0.01), N1 stage (p < 0.05), M0 stage (p < 0.01), Stage I–II (p < 0.05), and Stage III–IV (p < 0.01) (Supplementary Figure S4). Protein co-expression analysis for all 13 proteins was performed, and 69 proteins with correlation coefficients > 0.4 and a p value < 0.001 were identified (Figure 3G). The correlation of the proteins included in the prognostic signature is displayed in Figure 3H, in which Rheb and EVI1 showed the strongest positive correlation, while S100A4 and P90RSK_pT359S363 showed the strongest negative correlation.

3.4. Functional Enrichment Analysis Based on Risk Model

To identify the involved biological processes, we used GSEA to analyze protein expressions in THCA patients across different risk groups. In the high-risk group, representative KEGG pathways were chemokine signaling, cytokine–cytokine receptor interaction, ECM receptor interaction, focal adhesion, and the hematopoietic cell lineage pathway (Figure 4A). In the low-risk group, prominent KEGG pathways were fructose and mannose metabolism, linoleic acid metabolism, steroid hormone biosynthesis, and so on (Figure 4B). The differentially expressed proteins were further analyzed by GO and KEGG analysis. GO analysis revealed that primary functional categories in the biological processes (BP) were positive regulation of cell activation, regulation of DNA metabolic process and regulation of cell–cell adhesion. For cellular components (CCs), the major enriched GO terms were the chromosomal region, mismatch repair complex, and DNA repair complex. In terms of molecular functions (MFs), catalytic activity acting on DNA, phosphatidylinositol 3-kinase binding and phosphoprotein binding were enriched (Figure 4C,D, Supplementary Figure S5A,B). The KEGG pathway indicated that the differentially expressed proteins were mainly involved in EGFR tyrosine kinase inhibitor resistance, FoxO signaling central carbon metabolism in cancer, and the mTOR signaling pathway (Figure 4E,F, Supplementary Figure S5C,D).

To uncover the potential relationship between the DEPs, a PPI network that included 161 nodes and 256 edges was identified through the STRING tool. Then, the most significant modules were recognized by the MCODE plug-in of cytoscape. RBBP8 (CTIP), PMS2 (ATM), TP53BP1 (P53), MLH1 (MLH1), and MSH2 (MSH2) occupied the dominant position (Supplementary Figure S6A). GO analysis showed that the hub proteins were mainly engaged in response to stimulus and regulation of the biological process, metabolic process and immune system process (Supplementary Figure S6B).

3.5. Difference in Tumor-Infiltrating Immune Cells in Different Risk Groups

As immune micro-environmental abnormality plays an important role in oncogenesis, invasion, and metastasis, to further explore the relationship between the protein-based prognostic signature and immunity, we further evaluated the tumor micro-environmental (TME) characteristics between different risk groups. The result of TME showed that the stromal score (p < 0.05), immune score (p < 0.05), and ESTIMATE score (p < 0.05) were higher in the high-risk group than in the low-risk group (Figure 5A). CD8⁺ T cells, plasma cells, and activated NK cells were significantly elevated in the low-risk group (p < 0.05). In contrast, resting dendritic cells, activated dendritic cells, and neutrophils were obviously increased in the high-risk group (p < 0.05) (Figure 5B). The protein-associated signature was positively correlated with resting dendritic cells and activated dendritic cells (p < 0.001), but negatively correlated with activated NK cells (p < 0.01), plasma cells (p = 0.01), and CD8⁺ T cells (p < 0.001, Figure 5C). Differential and correlated analyses demonstrated consistency in the immunologic characteristics. Weak to moderate correlations were observed among various tumor-infiltrating immune cells (Figure 5D). Significant differences were found in the enrichment scores of different immune cell subpopulations (Figure 5E) and associated functions or pathways (Figure 5F) using ssGSEA, including HLA and checkpoint, across different risk groups. Thus, the heterogeneity of immune cell infiltration in THCA could be a novel indicator with a potential clinical significance for immunotherapy.

3.6. Protein-Based Signature Is Associated with Immunization Checkpoint Block

Cancer cells must have evaded the anti-tumor immune response to grow progressively, which relies in part on the expression on their surface of proteins with immunosuppressive functions, such as programmed death ligand 1 (PD-L1) [ref. 20]. The ability of cancer cells to evade immune detection often leads to malignant progression and poor outcome in THCA patients. To understand how the protein-based signature affects the prognosis, we assessed the differences in immune checkpoint blocks (ICBs) in THCA, including the levels of immune checkpoint genes and HLA component expression. To avoid bias from different THCA histological types, we counted all types in the dataset and found four cases of nonencapsulated sclerosing carcinomas, one case of oxyphilic adenocarcinoma, one case of follicular carcinoma, and the remaining were papillary thyroid carcinomas (PTCs). Since the sample size of THCA types other than PTCs was too small to be representative, we only evaluated the tumor immune dysfunction and exclusion (TIDE) scores in PTCs. We found that in the high-risk group, levels of immune checkpoint genes like HAVCR2, TNFSF4, CD274, TNFRSF18, CD28, CD276, CD80, CD86, CD70, CTLA4, BTLA, TNFSF9, ICOS, VTCN1, TNFSF18, TNFRSF9, TIGIT, and PDCD1LG2 were elevated, whereas KIR3DL1 and CD44 declined (Figure 6A). The correlation between immune checkpoint genes and risk scores also confirmed these results when we enrolled more immune checkpoint genes. Several genes in the signature including Rheb, EVI1, and SNAIL, were found to have negative correlations with nearly all checkpoints (Figure 6B). We further compared the differences in HLA expression. HLA-A, HLA-B, HLA-C, HLA-DMA, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DQB2, HLA-DRA, HLA-DRB1, HLA-DRB5, HLA-G, and HLA-H were evidently increased in the high-risk group (Figure 6C). Finally, we evaluated the potential immunotherapy response in each patient by the TIDE algorithm. However, the results showed that both PTC groups responded consistently to an immune checkpoint blockade with no difference (Figure 6D). The risk score tended to be higher in the PTC responder group than in the PTC non-responder group, yet not to a significant degree (Figure 6E).

3.7. Mutational Landscape Based on the Protein-Associated Signature

To further elucidate the mechanism underlying tumorigenesis and progression in THCA patients, we compared the tumor mutational burden (TMB) across different risk groups based on somatic mutation. We found that the high-risk group had a significantly higher TMB than the low-risk group (p = 0.016) (Supplementary Figure S7A). Additionally, the mutational landscape revealed a marked increase in BRAF mutation in the high-risk group (Supplementary Figure S7B). The overall survival of the low-TMB group was better compared to the high-TMB group (p < 0.001) (Supplementary Figure S7C). In a stratified survival analysis, we examined the combined impact of the TMB and prognosis signature, which confirmed that the TMB would not affect the predictive power of the signature for prognosis (p < 0.001) (Supplementary Figure S7D).

3.8. Potential Predictive Biomarker for Chemotherapy and Targeted Therapy

We investigated potential predictive biomarkers for chemotherapy and targeted therapy in all PTC samples and found that the high-risk group had a higher sensitivity to AGI-6780, AT13148, AZD4547, AZD5582, AZD5991, GDC0810, GSK591, GSK269962A, I-BET-762, I-BRD9, KRAS inhibitor, LCL161, Mirin, OSI-027, Rapamycin, Dactinomycin, Carmustine, Zoledronate, Pevonedistat, Tozasertib, Uprosertib, Axitinib, Dabrafenib, Crizotinib, Entospletinib, and Ibrutinib than the low-risk group (p < 0.001) (Figure 7, Supplementary Figure S8). Therefore, the established protein-associated signature might aid in choosing effective chemotherapeutic agents and developing personalized target drugs for PTC.

4. Discussion

Thyroid carcinoma, as the most common endocrine system malignant tumor, poses a threat to health and well-being, especially in patients who cannot be operated or who have a recurrence after an operation and who have no response to iodine therapy [ref. 21,ref. 22]. The mechanisms of tumorigenesis and progression of THCA are still unclear. There is an urgent demand to establish a novel prognosis signature and clarify potential mechanisms in THCA. Although the application of next-generation sequencing makes the analysis of RNA and DNA levels popular in oncology research, it should be noted that the majority of genes play their roles at the protein level. The development of proteomic technology and bioinformatics has provided new insights for comprehensive exploration of the molecular mechanisms of cancer pathogenesis and valuable information for prognosis and treatment decisions [ref. 23,ref. 24,ref. 25,ref. 26,ref. 27]. However, to date, there are relatively few studies on the role of proteomic models in THCA patients.

In our study, we randomly divided the entire set into a training and a testing set, and developed the model using the training set. After univariate Cox analysis and LASSO Cox regression analysis, we created a prognostic model composed of 13 proteins. Patients in the high-risk group showed significantly poor OS. The candidate proteins involved in the prognostic signature have been proven to engage in cancer development. S100A4 is a small calcium-binding protein that can affect multiple biological processes on its binding partners [ref. 28,ref. 29]. S100A4 is associated with poor survival in THCA as well as in glioma [ref. 30], bladder [ref. 31], pancreatic [ref. 32], and breast cancers [ref. 33]. PAI1 (SERPINE1) can not only inhibit caspase 3 [ref. 34], but also inhibit the cleavage of FasL and its abscission by the plasmin on the cell surface [ref. 35] to resist tumor cell apoptosis. IGFBP2 promotes tumor progression by inducing the alternative polarization of macrophages through the STAT3 pathway [ref. 36]. Emerging evidence has demonstrated that COL6A1 has diverse biological functions, including cell migration, differentiation, embryonic development, and maintenance of cell stemness in human malignancies [ref. 37]. Protease-activated receptors (PARs) are a unique family of G-protein coupled receptors, which play important roles in promoting cancer metastasis [ref. 38,ref. 39]. It has been widely reported that Snail played crucial roles in hepatocarcinogenesis and metastasis [ref. 40]. Several lines of evidence have indicated that the gap junction protein connexin 43 (Cx43) controls the response of glioblastoma to temozolomide (TMZ) [ref. 41], and blocking Cx43 using different approaches restores TMZ sensitivity [ref. 42,ref. 43]. Rictor correlates with poor prognosis of gastric cancer patients mainly through sensing growth factor concentration, regulating cell proliferation, survival, metabolism, and cytoskeletal remodeling [ref. 44]. In addition, our results suggest that B7-H3 was associated with increased risk in THCA patients, which was inconsistent with a previous report [ref. 45]. FAK expression is crucial for the development and growth of thyroid tumors. Recent research indicates that periostin, derived from cancer-associated fibroblasts, enhances thyroid tumor growth by activating FAK-STAT3 signaling [ref. 46]. Rheb is an important positive regulator of the mTORC1 pathway. Recent research indicates that Rheb degradation and subsequent mTORC1 inactivation play roles in cancer cell survival under glucose deprivation [ref. 47]. To date, the prognostic impact of EVI1 has not been clearly defined in real-world datasets. Some studies have found an association between high EVI1 expression and improved survival rates [ref. 48,ref. 49], which is consistent with our findings. High levels of phosphorylated p90RSK expression in ER-positive breast cancer tissues were linked to tumor shrinkage and decreased tumor volume following surgery. This association was particularly significant in ER-positive tumors [ref. 50]. The accuracy of the signature prediction for 1, 3, and 5 years of the three sets was more than 0.7. The combination of multiple proteins provided superior predictive efficiency compared to a single biomarker. The protein signature outperformed previously published signatures in predicting the prognosis of THCA [ref. 22,ref. 51,ref. 52]. Subgroup survival analysis, incorporating the clinical features, found that a high risk was associated with poorer prognosis in age > 65, female, T3–4 stage, N1 stage, M0 stage, and Stage III–IV subgroups. Additionally, multi-ROC curve analyses and the evaluation of the nomogram support the potential clinical application of this signature for prognosis prediction.

Expression analysis revealed that 245 proteins were differently expressed between the two risk groups. The results of GO and KEGG revealed that these differently expressed proteins were associated with the regulation of cell–cell adhesion, chromosomal region, PI3K binding, regulation of DNA metabolic process, EGFR tyrosine kinase inhibitor resistance, FoxO signaling pathway, NF-κB signaling pathway, PD-L1 expression, and PD-1 checkpoint pathway in cancer. Interestingly, these pathways were engaged in the progression of cancer. During malignant progression, epithelial cancer cells dissolve their cell–cell adhesion and gain invasive features [ref. 53]. Chromosomal rearrangements often result in active regulatory regions juxtaposed upstream of an oncogene to generate an expressed gene fusion [ref. 54]. PI3Ks are enzymes that belong to a family of plasma membrane-associated lipid kinases [ref. 55], implicated in the pathogenesis of many cancers [ref. 56]. NF-κB regulates the expression of genes involved in many processes that play a key role in the development and progression of cancer such as proliferation, migration, and apoptosis [ref. 57]. Epidermal growth factor receptor (EGFR), one of the most studied receptor kinases, is a drug target for cancer therapy, because its kinase activity correlates with tumorigenicity [ref. 58]. Emerging evidence has revealed that the expression of PD-L1 on tumor cells leads to immunosuppression and consequently enhances aggressiveness [ref. 59]. Recent studies have reported that PD-L1 is highly expressed in a subset of patients with advanced thyroid cancer and correlates with a higher risk of recurrence and shortened disease-free survival [ref. 60,ref. 61]. These differences may contribute to the different prognosis between the two groups.

The tumor immune microenvironment plays a key role in regulating the processes of tumor development, invasion, and metastasis [ref. 62]. To gain deeper insights into the immune cell infiltration in THCA, we explored the relationship between the protein-based prognostic signature and immune cell infiltration. We found that resting dendritic cells, activated dendritic cells, and neutrophils were significantly elevated in the high-risk group, while CD8⁺ T cells, plasma cells, and activated NK cells were remarkably elevated in the low-risk group. CD8⁺ T cells are typically regarded as a uniform population of cells that secrete large amounts of IFN-γ and the protease granzyme B, which act synergistically to kill infected or tumorigenic cells [ref. 63]. Plasma cells play a beneficial role in the majority of cancer types [ref. 64], suggesting the low-risk group THCA patients benefit from the infiltrating plasma cells. There is emerging evidence of increased NK cell-mediated tumor cell killing [ref. 65]. In THCA patients, tumor-associated neutrophils, which are capable of promoting tumor progression by secreting a wide range of cytokines [ref. 66], correlated with a larger tumor size [ref. 67]. Several studies have shown that tumor-infiltrating dendritic cells are related to more advanced tumor T stage and lymph node metastasis in THCA [ref. 68,ref. 69]. Our results indicate increased anti-tumor immune activity in the low-risk group, partially explaining the predictive value of the prognostic signature. The immune cell subpopulation enrichment analysis and the ssGSEA-related functions or pathway analysis revealed that scores for nearly all the items were higher in the high-risk group, which indicated that THCA patients in the high-risk group seemed to be more immunogenic. Moreover, the high-risk group had a higher immune score, stromal score, and ESTIMATE score, which further suggests that patients with a high-risk score appeared to be more immunogenic and might benefit from immunotherapy.

We further evaluated the differences in the ICB between the two groups. The high-risk groups featured higher levels of HLA family genes and higher levels of immune checkpoint genes. Lymphocyte-activation gene 3 (LAG3) is an immune inhibitory checkpoint and is expressed on the surface of lymphocytes [ref. 70], such as CD4+ T cells, CD8+ T cells, natural killer (NK) cells, NK T (NKT) cells, and regulatory T (Treg) cells. LAG3 inhibits the tumor immune microenvironment by accelerating T cell exhaustion and blocking T cell proliferation [ref. 71]. Programmed death-ligand 1 (PD-L1, CD274), an immune checkpoint protein expressed on cancer cells and immune cells, interacts with its binding partner programmed cell death-1 (PD-1) to inhibit T cell proliferation and cytokine production [ref. 72]. TIGIT in CD4+ T cells induce immunosuppression through inhibiting T cell proliferation directly by inducing the down-expression of T-bet, GATA3, IRF4, and RORc, which reduce the level of pro-inflammatory IFN-γ while increasing the level of anti-inflammatory IL-10 [ref. 73]. PDCD1LG2 (PD-L2) is a second ligand for PD-1 and inhibits T cell activation [ref. 74]. The immune inhibitory checkpoint contributed to the worse survival in the high-risk groups. On the other hand, MHC molecules were found to be overexpressed in the high-risk groups. Tumor recognition by immune cells requires presentation of non-self peptides (neoantigens) by tumor cells through MHC Class I or II complexes. Loss or reduced expression of MHC or their subunits abrogates T cell-mediated anti-tumor immunity [ref. 75]. Tumor-specific MHC-II is associated with good outcomes in cancer patients, including those with immunotherapies [ref. 76]. The higher expression of MHC molecules in THCA patients of the high-risk group may explain the better prognosis of THCA than other tumors. Our results suggest that patients in the high-risk group have a higher ICB, which may have a possible good anti-tumor immune response.

To delve into the cause of tumorigenesis, the TMB of two groups was compared and a higher TMB was found in the high-risk group. In our present study, the most frequent mutation was the BRAF missense mutation. Previous studies have found the BRAF T1799A mutation in approximately 45% of PTC and 25% of apparently PTC-derived anaplastic thyroid cancers, but not in follicular thyroid cancer (FTC) and benign thyroid tumors [ref. 77]. BRAF T1799A mutation contributes to poor clinicopathologic outcomes of PTC [ref. 78] and this mutation is often correlated with a loss of radiotherapy response and recurrence of PTC [ref. 77]. The BRAF V600E mutant has increased kinase activity, resulting in promoting cell growth, differentiation, and survival [ref. 79]. Kaplan–Meier analysis shows that the low-TMB group had better OS than the high-TMB group; in addition, the combined Kaplan–Meier analysis demonstrated that the risk score was a reliable predictive indicator independent of the TMB. Patients with a high TMB have been reported to have a more sensitive response to immune checkpoint blockade, such as anti-PD-1 agents [ref. 80]. The risk score tended to be higher in the ICB responder group than in the non-responder group, although not to a significant degree.

Effectiveness of chemotherapy and targeted treatment between the two risk groups were predicted by the GDSC dataset. In the high-risk group, the THCA patients were more sensitive to AGI-6780, AT13148, AZD4547, AZD5582, AZD5991, GDC0810, GSK591, GSK269962A, I-BET-762, I-BRD9, KRAS inhibitor, LCL161, Mirin, OSI-027, Rapamycin, Dactinomycin, Carmustine, Zoledronate, Pevonedistat, Tozasertib, Uprosertib, Axitinib, Dabrafenib, Crizotinib, Entospletinib, and Ibrutinib, which might be a potential treatment option for THCA. Our results suggest that THCA patients in the high-risk groups had more ICB expressions, a higher TMB, and PTC patients in the high-risk groups responded more to sensitivity to chemotherapy and targeted treatment.

However, this study has some limitations that should be considered. To begin with, it is a preliminary study based on bioinformatics tools, incorporating real-world samples; other cohort validation will allow for a more rigorous design. In addition, our molecular results lack experimental evidence in vivo and in vitro.

5. Conclusions

In summary, we constructed a valid prognostic protein signature to predict outcomes for THCA patients, demonstrating a strong predictive ability. We also evaluated the differences in immunotherapy response and chemotherapeutic drug sensitivity between the two risk groups. These results could advance our understanding of THCA pathogenesis and provide new strategies for personalized therapy.

References

1. R.L. Siegel, K.D. Miller, H.E. Fuchs, A. Jemal. Cancer Statistics, 2021. CA Cancer J. Clin., 2021. [DOI | PubMed]
2. J. Kim, J.E. Gosnell, S.A. Roman. Geographic influences in the global rise of thyroid cancer. Nat. Rev. Endocrinol., 2020. [DOI | PubMed]
3. H. Lim, S.S. Devesa, J.A. Sosa, D. Check, C.M. Kitahara. Trends in Thyroid Cancer Incidence and Mortality in the United States, 1974–2013. JAMA, 2017. [DOI | PubMed]
4. V. Drozd, R. Schneider, T. Platonova, G. Panasiuk, T. Leonova, N. Oculevich, I. Shimanskaja, I. Vershenya, T. Dedovich, T. Mitjukova. Feasibility Study Shows Multicenter, Observational Case-Control Study Is Practicable to Determine Risk of Secondary Breast Cancer in Females with Differentiated Thyroid Carcinoma Given Radioiodine Therapy in Their Childhood or Adolescence; Findings Also Suggest Possible Fertility Impairment in Such Patients. Front. Endocrinol., 2020
5. C.L. Jillard, R.P. Scheri, J.A. Sosa. What Is the Optimal Treatment of Papillary Thyroid Cancer?. Adv. Surg., 2015. [DOI | PubMed]
6. L.G.T. Morris, R.M. Tuttle, L. Davies. Changing Trends in the Incidence of Thyroid Cancer in the United States. JAMA Otolaryngol. Head. Neck Surg., 2016. [DOI | PubMed]
7. C.D. Seib, J.A. Sosa. Evolving Understanding of the Epidemiology of Thyroid Cancer. Endocrinol. Metab. Clin. North. Am., 2019. [DOI | PubMed]
8. R. Stark, M. Grzelak, J. Hadfield. RNA sequencing: The teenage years. Nat. Rev. Genet., 2019. [DOI | PubMed]
9. V. Subbiah, D. Yang, V. Velcheti, A. Drilon, F. Meric-Bernstam. State-of-the-Art Strategies for Targeting RET-Dependent Cancers. J. Clin. Oncol., 2020. [DOI | PubMed]
10. S.F. Idris, S.S. Ahmad, M.A. Scott, G.S. Vassiliou, J. Hadfield. The role of high-throughput technologies in clinical cancer genomics. Expert. Rev. Mol. Diagn., 2013. [DOI | PubMed]
11. Y. Liu, A. Beyer, R. Aebersold. On the Dependency of Cellular Protein Levels on mRNA Abundance. Cell, 2016. [DOI | PubMed]
12. F. Edfors, F. Danielsson, B.M. Hallström, L. Käll, E. Lundberg, F. Pontén, B. Forsström, M. Uhlén. Gene-specific correlation of RNA and protein levels in human cells and tissues. Mol. Syst. Biol., 2016. [DOI | PubMed]
13. D.P. Nusinow, J. Szpyt, M. Ghandi, C.M. Rose, E.R. McDonald, M. Kalocsay, J. Jané-Valbuena, E. Gelfand, D.K. Schweppe, M. Jedrychowski. Quantitative Proteomics of the Cancer Cell Line Encyclopedia. Cell, 2020. [DOI | PubMed]
14. F. Finotello, D. Rieder, H. Hackl, Z. Trajanoski. Next-generation computational tools for interrogating cancer immunity. Nat. Rev. Genet., 2019. [DOI | PubMed]
15. H. Wang, B.J. Lengerich, B. Aragam, E.P. Xing. Precision Lasso: Accounting for correlations and linear dependencies in high-dimensional genomic data. Bioinformatics, 2019. [DOI | PubMed]
16. J. Friedman, T. Hastie, R. Tibshirani. Regularization Paths for Generalized Linear Models via Coordinate Descent. J. Stat. Softw., 2010. [DOI | PubMed]
17. Y. Zhou, B. Zhou, L. Pache, M. Chang, A.H. Khodabakhshi, O. Tanaseichuk, C. Benner, S.K. Chanda. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun., 2019. [DOI | PubMed]
18. P. Shannon, A. Markiel, O. Ozier, N.S. Baliga, J.T. Wang, D. Ramage, N. Amin, B. Schwikowski, T. Ideker. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res., 2003. [DOI | PubMed]
19. G.D. Bader, C.W.V. Hogue. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinform., 2003. [DOI]
20. T. Reisländer, F.J. Groelly, M. Tarsounas. DNA Damage and Cancer Immunotherapy: A STING in the Tale. Mol. Cell, 2020. [DOI | PubMed]
21. N. George, A. Agarwal, N. Kumari, S. Agarwal, N. Krisnani, S.K. Gupta. Mutational Profile of Papillary Thyroid Carcinoma in an Endemic Goiter Region of North India. Indian. J. Endocrinol. Metab., 2018. [PubMed]
22. R. Qin, C. Li, X. Wang, Z. Zhong, C. Sun. Identification and validation of an immune-related prognostic signature and key gene in papillary thyroid carcinoma. Cancer Cell Int., 2021. [DOI | PubMed]
23. Y.-A. Chen, C.-Y. Lu, W.-F. Cheng, K.-T. Kuo, C.-W. Yu, H.-N. Ho, H.-F. Chen, S.-H. Pan. An experimental model for ovarian cancer: Propagation of ovarian cancer initiating cells and generation of ovarian cancer organoids. BMC Cancer, 2022. [DOI | PubMed]
24. Y. Qu, J. Feng, X. Wu, L. Bai, W. Xu, L. Zhu, Y. Liu, F. Xu, X. Zhang, G. Yang. A proteogenomic analysis of clear cell renal cell carcinoma in a Chinese population. Nat. Commun., 2022. [DOI | PubMed]
25. W. Chong, X. Zhu, H. Ren, C. Ye, K. Xu, Z. Wang, S. Jia, L. Shang, L. Li, H. Chen. Integrated multi-omics characterization of KRAS mutant colorectal cancer. Theranostics, 2022. [DOI | PubMed]
26. S. Saha, R. Verma, C. Kumar, B. Kumar, A.K. Dey, M. Surjit, S.V.S. Mylavarapu, T.K. Maiti. Proteomic analysis reveals USP7 as a novel regulator of palmitic acid-induced hepatocellular carcinoma cell death. Cell Death Dis., 2022. [DOI | PubMed]
27. C.K.Y. Ng, E. Dazert, T. Boldanova, M. Coto-Llerena, S. Nuciforo, C. Ercan, A. Suslov, M.-A. Meier, T. Bock, A. Schmidt. Integrative proteogenomic characterization of hepatocellular carcinoma across etiologies and stages. Nat. Commun., 2022. [DOI | PubMed]
28. N. Ambartsumian, J. Klingelhöfer, M. Grigorian. The Multifaceted S100A4 Protein in Cancer and Inflammation. Methods Mol. Biol., 2019. [PubMed]
29. F. Fei, J. Qu, M. Zhang, Y. Li, S. Zhang. S100A4 in cancer progression and metastasis: A systematic review. Oncotarget, 2017. [DOI | PubMed]
30. N. Abdelfattah, P. Kumar, C. Wang, J.-S. Leu, W.F. Flynn, R. Gao, D.S. Baskin, K. Pichumani, O.B. Ijare, S.L. Wood. Single-cell analysis of human glioma and immune cells identifies S100A4 as an immunotherapy target. Nat. Commun., 2022. [DOI | PubMed]
31. B.R. Davies, M. O’Donnell, G.C. Durkan, P.S. Rudland, R. Barraclough, D.E. Neal, J.K. Mellon. Expression of S100A4 protein is associated with metastasis and reduced survival in human bladder cancer. J. Pathol., 2002. [DOI | PubMed]
32. S.H. Lee, H. Kim, J.-H. Hwang, E. Shin, H.S. Lee, D.W. Hwang, J.Y. Cho, Y.-S. Yoon, H.-S. Han, B.H. Cha. CD24 and S100A4 expression in resectable pancreatic cancers with earlier disease recurrence and poor survival. Pancreas, 2014. [DOI | PubMed]
33. S. de Silva Rudland, L. Martin, C. Roshanlall, J. Winstanley, S. Leinster, A. Platt-Higgins, J. Carroll, C. West, R. Barraclough, P. Rudland. Association of S100A4 and osteopontin with specific prognostic factors and survival of patients with minimally invasive breast cancer. Clin. Cancer Res., 2006. [DOI | PubMed]
34. U.A. Lademann, M.U. Rømer. Regulation of programmed cell death by plasminogen activator inhibitor type 1 (PAI-1). Thromb. Haemost., 2008. [PubMed]
35. M. Valiente, A.C. Obenauf, X. Jin, Q. Chen, X.H.F. Zhang, D.J. Lee, J.E. Chaft, M.G. Kris, J.T. Huse, E. Brogi. Serpins promote cancer cell survival and vascular co-option in brain metastasis. Cell, 2014. [DOI | PubMed]
36. L. Sun, X. Zhang, Q. Song, L. Liu, E. Forbes, W. Tian, Z. Zhang, Y.A. Kang, H. Wang, J.B. Fleming. IGFBP2 promotes tumor progression by inducing alternative polarization of macrophages in pancreatic ductal adenocarcinoma through the STAT3 pathway. Cancer Lett., 2021. [DOI | PubMed]
37. M. Cescon, F. Gattazzo, P. Chen, P. Bonaldo. Collagen VI at a glance. J. Cell Sci., 2015. [DOI | PubMed]
38. D. Fujimoto, Y. Hirono, T. Goi, K. Katayama, A. Yamaguchi. Prognostic value of protease-activated receptor-1 (PAR-1) and matrix metalloproteinase-1 (MMP-1) in gastric cancer. Anticancer. Res., 2008. [PubMed]
39. J.Y. Suen, A. Cotterell, R.J. Lohman, J. Lim, A. Han, M.K. Yau, L. Liu, M.A. Cooper, D.A. Vesey, D.P. Fairlie. Pathway-selective antagonism of proteinase activated receptor 2. Br. J. Pharmacol., 2014. [DOI | PubMed]
40. X.-Y. Huang, P.-F. Zhang, C.-Y. Wei, R. Peng, J.-C. Lu, C. Gao, J.-B. Cai, X. Yang, J. Fan, A.-W. Ke. Circular RNA circMET drives immunosuppression and anti-PD1 therapy resistance in hepatocellular carcinoma via the miR-30-5p/snail/DPP4 axis. Mol. Cancer, 2020. [DOI | PubMed]
41. J.L. Munoz, V. Rodriguez-Cruz, S.J. Greco, S.H. Ramkissoon, K.L. Ligon, P. Rameshwar. Temozolomide resistance in glioblastoma cells occurs partly through epidermal growth factor receptor-mediated induction of connexin 43. Cell Death Dis., 2014. [DOI | PubMed]
42. L. Wang, Y. Wang, A. Chen, A. Jalali, S. Liu, Y. Guo, S. Na, H. Nakshatri, B.-Y. Li, H. Yokota. Effects of a checkpoint kinase inhibitor, AZD7762, on tumor suppression and bone remodeling. Int. J. Oncol., 2018. [DOI | PubMed]
43. X.-H. Zhang, Y. Qian, Z. Li, N.-N. Zhang, Y.-J. Xie. Let-7g-5p inhibits epithelial-mesenchymal transition consistent with reduction of glioma stem cell phenotypes by targeting VSIG4 in glioblastoma. Oncol. Rep., 2016. [DOI | PubMed]
44. Y. Bian, Z. Wang, J. Xu, W. Zhao, H. Cao, Z. Zhang. Elevated Rictor expression is associated with tumor progression and poor prognosis in patients with gastric cancer. Biochem. Biophys. Res. Commun., 2015. [DOI | PubMed]
45. B. Zhao, Z. Huang, X. Zhu, H. Cai, Y. Huang, X. Zhang, Z. Zhang, H. Lu, C. An, L. Niu. Clinical Significance of the Expression of Co-Stimulatory Molecule B7-H3 in Papillary Thyroid Carcinoma. Front. Cell Dev. Biol., 2022. [DOI | PubMed]
46. X. Jin, Q. Deng, S. Ye, S. Liu, Y. Fu, Y. Liu, G. Wu, G. Ouyang, T. Wu. Cancer-associated fibroblast-derived periostin promotes papillary thyroid tumor growth through integrin-FAK-STAT3 signaling. Theranostics, 2024. [DOI | PubMed]
47. M. Li, W. Huang, Y. Zhang, Y. Du, S. Zhao, L. Wang, Y. Sun, B. Sha, J. Yan, Y. Ma. Glucose deprivation triggers DCAF1-mediated inactivation of Rheb-mTORC1 and promotes cancer cell survival. Cell Death Dis., 2024. [DOI | PubMed]
48. M. Nanjundan, Y. Nakayama, K.W. Cheng, J. Lahad, J. Liu, K. Lu, W.-L. Kuo, K. Smith-McCune, D. Fishman, J.W. Gray. Amplification of MDS1/EVI1 and EVI1, located in the 3q26.2 amplicon, is associated with favorable patient prognosis in ovarian cancer. Cancer Res., 2007. [DOI | PubMed]
49. P. Jank, J. Leichsenring, S. Kolb, I. Hoffmann, P. Bischoff, C.A. Kunze, M.P. Dragomir, M. Gleitsmann, M. Jesinghaus, W.D. Schmitt. High EVI1 and PARP1 expression as favourable prognostic markers in high-grade serous ovarian carcinoma. J. Ovarian Res., 2023. [DOI | PubMed]
50. H.-G. Moon, J.K. Yi, H.S. Kim, H.Y. Lee, K.-M. Lee, M. Yi, S. Ahn, H.-C. Shin, J.-h. Ju, I. Shin. Phosphorylation of p90RSK is associated with increased response to neoadjuvant chemotherapy in ER-positive breast cancer. BMC Cancer, 2012. [DOI | PubMed]
51. B. Wang, W. Jiang, X. Zheng, Y. Han, R. Liu. Research on a Weighted Gene Co-expression Network Analysis method for mining pathogenic genes in thyroid cancer. PLoS ONE, 2022. [DOI | PubMed]
52. Y. Xia, X. Jiang, Y. Huang, Q. Liu, Y. Huang, B. Zhang, Z. Mei, D. Xu, Y. Shi, W. Tu. Construction of a Tumor Immune Microenvironment-Related Prognostic Model in BRAF-Mutated Papillary Thyroid Cancer. Front. Endocrinol., 2022. [DOI | PubMed]
53. D. Buechel, N. Sugiyama, N. Rubinstein, M. Saxena, R.K.R. Kalathur, F. Lüönd, V. Vafaizadeh, T. Valenta, G. Hausmann, C. Cantù. Parsing β-catenin’s cell adhesion and Wnt signaling functions in malignant mammary tumor progression. Proc. Natl. Acad. Sci. USA, 2021. [DOI | PubMed]
54. A. Nickless, J. Zhang, G. Othoum, J. Webster, M.J. Inkman, E. Coonrod, S. Fontes, E.B. Rozycki, C.A. Maher, N.M. White. Pan-Cancer Analysis Reveals Recurrent BCAR4 Gene Fusions across Solid Tumors. Mol. Cancer Res., 2022. [DOI | PubMed]
55. T.R. Donahue, L.M. Tran, R. Hill, Y. Li, A. Kovochich, J.H. Calvopina, S.G. Patel, N. Wu, A. Hindoyan, J.J. Farrell. Integrative survival-based molecular profiling of human pancreatic cancer. Clin. Cancer Res., 2012. [DOI | PubMed]
56. T.L. Yuan, L.C. Cantley. PI3K pathway alterations in cancer: Variations on a theme. Oncogene, 2008. [DOI | PubMed]
57. X. Dolcet, D. Llobet, J. Pallares, X. Matias-Guiu. NF-kB in development and progression of human cancer. Virchows Arch., 2005. [DOI | PubMed]
58. A. De Luca, A. Carotenuto, A. Rachiglio, M. Gallo, M.R. Maiello, D. Aldinucci, A. Pinto, N. Normanno. The role of the EGFR signaling in tumor microenvironment. J. Cell Physiol., 2008. [DOI | PubMed]
59. B. Sasi, P. Ethiraj, J. Myers, A.-P. Lin, S. Jiang, Z. Qiu, K.N. Holder, R.C.T. Aguiar. Regulation of PD-L1 expression is a novel facet of cyclic-AMP-mediated immunosuppression. Leukemia, 2021. [DOI | PubMed]
60. S. Ahn, T.H. Kim, S.W. Kim, C.S. Ki, H.W. Jang, J.S. Kim, J.H. Kim, J.-H. Choe, J.H. Shin, S.Y. Hahn. Comprehensive screening for PD-L1 expression in thyroid cancer. Endocr. Relat. Cancer, 2017. [DOI | PubMed]
61. R.-L. Shi, N. Qu, T.-X. Luo, J. Xiang, T. Liao, G.-H. Sun, Y. Wang, Y.-L. Wang, C.-P. Huang, Q.-H. Ji. Programmed Death-Ligand 1 Expression in Papillary Thyroid Cancer and Its Correlation with Clinicopathologic Factors and Recurrence. Thyroid., 2017. [DOI | PubMed]
62. Y. Kurebayashi, H. Ojima, H. Tsujikawa, N. Kubota, J. Maehara, Y. Abe, M. Kitago, M. Shinoda, Y. Kitagawa, M. Sakamoto. Landscape of immune microenvironment in hepatocellular carcinoma and its additional impact on histological and molecular classification. Hepatology, 2018. [DOI | PubMed]
63. M. St Paul, P.S. Ohashi. The Roles of CD8+ T Cell Subsets in Antitumor Immunity. Trends Cell Biol., 2020. [DOI | PubMed]
64. M.C.A. Wouters, B.H. Nelson. Prognostic Significance of Tumor-Infiltrating B Cells and Plasma Cells in Human Cancer. Clin. Cancer Res., 2018. [DOI | PubMed]
65. N. Shimasaki, A. Jain, D. Campana. NK cells for cancer immunotherapy. Nat. Rev. Drug Discov., 2020. [DOI | PubMed]
66. C. Lu, D. Rong, B. Zhang, W. Zheng, X. Wang, Z. Chen, W. Tang. Current perspectives on the immunosuppressive tumor microenvironment in hepatocellular carcinoma: Challenges and opportunities. Mol. Cancer, 2019. [DOI | PubMed]
67. M.R. Galdiero, G. Varricchi, S. Loffredo, C. Bellevicine, T. Lansione, A.L. Ferrara, R. Iannone, S. di Somma, F. Borriello, E. Clery. Potential involvement of neutrophils in human thyroid cancer. PLoS ONE, 2018. [DOI | PubMed]
68. K. Bergdorf, D.C. Ferguson, M. Mehrad, K. Ely, T. Stricker, V.L. Weiss. Papillary thyroid carcinoma behavior: Clues in the tumor microenvironment. Endocr. Relat. Cancer, 2019. [DOI | PubMed]
69. Z. Xie, X. Li, Y. He, S. Wu, S. Wang, J. Sun, Y. He, Y. Lun, J. Zhang. Immune Cell Confrontation in the Papillary Thyroid Carcinoma Microenvironment. Front. Endocrinol., 2020. [DOI]
70. C. Solinas, E. Migliori, P. De Silva, K. Willard-Gallo. LAG3: The Biological Processes That Motivate Targeting This Immune Checkpoint Molecule in Human Cancer. Cancers, 2019. [DOI | PubMed]
71. S.C. Wei, C.R. Duffy, J.P. Allison. Fundamental Mechanisms of Immune Checkpoint Blockade Therapy. Cancer Discov., 2018. [DOI | PubMed]
72. L. Ai, A. Xu, J. Xu. Roles of PD-1/PD-L1 Pathway: Signaling, Cancer, and Beyond. Adv. Exp. Med. Biol., 2020. [PubMed]
73. E. Lozano, M. Dominguez-Villar, V. Kuchroo, D.A. Hafler. The TIGIT/CD226 axis regulates human T cell function. J. Immunol., 2012. [DOI | PubMed]
74. Y. Latchman, C.R. Wood, T. Chernova, D. Chaudhary, M. Borde, I. Chernova, Y. Iwai, A.J. Long, J.A. Brown, R. Nunes. PD-L2 is a second ligand for PD-1 and inhibits T cell activation. Nat. Immunol., 2001. [DOI | PubMed]
75. A. Balasubramanian, T. John, M.-L. Asselin-Labat. Regulation of the antigen presentation machinery in cancer and its implication for immune surveillance. Biochem. Soc. Trans., 2022. [DOI | PubMed]
76. M.L. Axelrod, R.S. Cook, D.B. Johnson, J.M. Balko. Biological Consequences of MHC-II Expression by Tumor Cells in Cancer. Clin. Cancer Res., 2019. [DOI | PubMed]
77. M. Xing, W.H. Westra, R.P. Tufano, Y. Cohen, E. Rosenbaum, K.J. Rhoden, K.A. Carson, V. Vasko, A. Larin, G. Tallini. BRAF mutation predicts a poorer clinical prognosis for papillary thyroid cancer. J. Clin. Endocrinol. Metab., 2005. [DOI | PubMed]
78. M. Xing. BRAF mutation in papillary thyroid cancer: Pathogenic role, molecular bases, and clinical implications. Endocr. Rev., 2007. [DOI | PubMed]
79. H. Davies, G.R. Bignell, C. Cox, P. Stephens, S. Edkins, S. Clegg, J. Teague, H. Woffendin, M.J. Garnett, W. Bottomley. Mutations of the BRAF gene in human cancer. Nature, 2002. [DOI | PubMed]
80. L. Fancello, S. Gandini, P.G. Pelicci, L. Mazzarella. Tumor mutational burden quantification from targeted gene panels: Major advancements and challenges. J. Immunother. Cancer, 2019. [DOI | PubMed]