Association of a Novel Prognosis Model with Tumor Mutation Burden and Tumor-Infiltrating Immune Cells in Thyroid Carcinoma

Although immunotherapy has recently demonstrated a substantial promise in treating advanced thyroid carcinoma (THCA), it is not appropriate for all THCA patients. As a result, this study aims to identify biomarkers for predicting immunotherapy efficacy and prognosis in THCA patients based on a constructed prognostic model. The transcriptomic and corresponding clinical data of THCA patients were obtained from the Cancer Genome Atlas (TCGA) database. We identified differentially expressed genes (DEGs) between THCA and normal samples and performed an intersection analysis of DEGs with immune-related genes (IRGs) downloaded from the ImmPort database. Functional enrichment analysis was performed on the chosen immune-related DEGs. Subsequently, Cox and LASSO regression analyses were conducted to obtain three hub immune-related DEGs, including PPBP, SEMA6B, and GCGR. Following that, a prognostic risk model was established and validated based on PPBP, SEMA6B, and GCGR genes to predict immunotherapy efficacy and THCA prognosis. Finally, we investigated the association between the constructed risk model and tumor mutational burden (TMB), abundance of tumor-infiltrating immune cells (TICs) as well as immunotherapeutic targets (PDL-1, PD-1, and CTLA4) in THCA. THCA patients in the high-risk score (RS) group showed higher TMB levels and worse prognosis than the low RS group. Patients in the high-RS group had higher proportions of monocytes, M2 macrophages, and activated dendritic cells, whereas those in the low-RS group exhibited higher numbers of M1 macrophages and dendritic resting cells. Our data implied that the constructed THCA prognostic model was sound and we concluded that the THCA patients having high TMB and low PD-L1 expression levels might respond poorly to immunotherapy. Taken together, we constructed a novel prognostic model for THCA patients to predict their prognosis and immunotherapy efficacy, providing a viable option for the future management of THCA patients in the clinic.

Introduction

Thyroid carcinoma (THCA) is the fifth most prevalent malignancy affecting women (Zhang et al., 2019b) and a major cause of annual endocrine malignancies death (Aboelnaga and Ahmed., 2015). Over the past decades, the incidence of THCA has been escalating worldwide (Cabanillas et al., 2016). THCA is generally classified into four pathological subtypes: papillary thyroid cancer (PTC, 80–85%), follicular thyroid cancer (FTC, 10–15%), medullary thyroid carcinoma (MTC, less than 2%) and anaplastic thyroid carcinoma (ATC, less than 2%) (Laha et al., 2020). Although most THCA is well-differentiated PTC with a 10-years survival rate of over 95%, some variants of PTC may demonstrate increased aggressive behavior, particularly in older patients, contributing to significant mortality (Nath and Erickson., 2018). Therefore, surgical resection and radioiodine (RAI) therapy were considered standard treatments for most THCA patients. Nonetheless, survival rates for advanced thyroid cancer patients remain low. For these reasons, novel therapeutic strategies, such as immunotherapies and targeted molecular therapy, are under investigation for treating advanced or metastatic THCA patients.

In recent years, immune checkpoint inhibitor (ICI) has achieved remarkable progress in treating breast cancer (Santa-Maria and Nanda., 2018), lung cancer (Anagnostou et al., 2017), melanoma (Amaria et al., 2018; Barrios et al., 2020) and squamous cell carcinoma of the head and neck (HNSCC) (Ferris et al., 2016; Ferris et al., 2018). It has been shown that numerous immune cells and their mediators exist in the tumor microenvironment (TME) of THCA with various interactions between them (Varricchi et al., 2019). In addition, it was demonstrated that increased frequency of regulatory and PD-1⁺ T cells was associated with the recurrent or aggressive PTC (French et al., 2012), and that the frequency of PD-1⁺ T cells was higher in patients with extranodal invasion than in those without lymph nodes metastases, indicating that it may be associated with the THCA prognosis. Furthermore, Gunda and his colleagues have recently revealed that anti-PD-1/PD-L1 therapy could beneficially modulate the immune microenvironment in orthotopic ATC murine model while simultaneously enhancing the efficacy of lenvatinib, a multi-targeted tyrosine kinase inhibitor (Gunda et al., 2019). The efficacy of ICIs has also yielded encouraging results in clinical trials, manifesting potent antitumor effects and improved tolerability in advanced FTC patients (Mehnert et al., 2019). However, the expensive cost of ICIs, which averages around $150,000 per year (Oiseth and Aziz., 2017), and their unpredictable efficacies (Giannone et al., 2020), prelude their widespread clinical use. Considering these factors, screening for appropriate molecular markers to improve treatment precision is highly demanding.

Tumor mutational burden (TMB), or the number of nonsynonymous mutations in a genomic region of somatic cells, is a biomarker employed to predict immunotherapy efficacy in various cancers (Samstein et al., 2019). As is widely known, tumor-infiltrating immune cells (TICs), an essential component of TME, are critical for in tumor initiation and progression (Bissell and Hines., 2011). Furthermore, gene locus mutations are observed in many histological subtypes of THCA, such as anaplastic lymphoma kinase (ALK), neurotrophic receptor tyrosine kinase 1 (NTRK1) genes rearrangements, BRAF and GTPase RAS family genes (Bos, 1989; Nikiforov and Nikiforova., 2011; Cancer Genome Atlas Research Network, 2014; Arndt et al., 2018; Varricchi et al., 2019). As high-throughput sequencing advances, large-scale acquisition of relevant cancer genomic data has become possible. However, few studies have conducted in-depth investigations to evaluate immunotherapy efficacy and prognosis in THCA.

Herein, based on Cancer Genome Atlas (TCGA) and ImmPort databases, we screened immune-related differentially expressed genes (IRDEGs), explored their functional enrichment, and constructed a prognostic prediction model. Additionally, we further investigated the association of TMB with immune infiltration and prognosis in THCA.

Materials and Methods

Flow Diagram of Analysis

We designed a flow chart of analysis for the construction, validation and evaluation of the prognostic model in THCA. The analysis process was performed strictly according to the flow chart (Figure 1).

FIGURE 1

FIGURE 1. Flow chart for construction and validation the THCA prognostic model. THCA, thyroid carcinoma.

Data Downloading and Processing

We downloaded transcriptomic data with HTSeq-FPKM workflow type from TCGA of THCA project database (https://portal.gdc.cancer.gov/), including 510 tumor and 58 normal samples. Besides, masked somatic mutation data processed by VarScan software was also acquired from the GDC (Genomic Data Commons) data portal of TCGA database. Next, we downloaded clinical data of corresponding patients, including their age, gender, tumor grade, tumor stage, survival time and survival status.

Screening of Immune-Related Differentially Expressed Genes

Firstly, differentially expressed genes (DEGs) between normal and tumor tissues of THCA were analyzed using R software “limma” package⁴, and the screening criteria were false discovery rate (FDR) < 0.05 and |log2FC| > 1. Immune-related genes (IRGs) were obtained from ImmPort database (http://www.immport.org/) and then performed an intersection analysis on DEGs and IRGs to identify IRDEGs. R “venn” package was used to visualize to result of intersection analysis and R “pheatmap” package was performed to manifest the expression levels of IRDEGs of the THCA samples in the TCGA database.

Functional Enrichment Analysis

We used R package “org.Hs.eg.db” to obtain Entrez-ID of each IRDEG, and then performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses on IRDEGs using “clusterProfiler,” “enrichplot,” and “ggplot2” R packages. p < 0.05 and q < 0.05 were utilized as statistical significance thresholds.

Construction and Validation of a Prediction Model

For cross-validation, we randomly divided all samples into two groups with a 2:1 ratio using R package “caret”, referred to as the training cohort and the test cohort, containing 2/3 and 1/3 of THCA cases, respectively. First, we determined IRDEGs associated with THCA survival in the training cohort using univariate Cox regression analysis with a threshold value of p < 0.001 via R package “survival”. Then, key IRDEGs were identified using LASSO regression using R package “glmnet”, and multivariate Cox regression analyses, thus constructing a THCA prediction model. Following that, we divided all THCA samples into high- and low-risk score (RS) groups according to their median RS. For high- and low-RS groups, Kaplan-Meier survival analysis and time-related receiver operating characteristic (ROC) curve analysis through R package “timeROC”, were employed to cross-validate the predictive power of the model in the training test and the combined cohorts (containing all THCA samples).

Analysis of Tumor-Infiltrating Immune Cells

We utilized CIBERSORT algorithm (R script v 1.03), a deconvolution algorithm to quantify immune cell proportions based on transcriptomic expression profiles (Newman et al., 2015; Yao et al., 2020), to calculate the relative frequencies of 22 tumor-infiltrating immune cells (TIC) subtypes in tumor samples with high and low RS groups. p < 0.05 was utilized as statistical significance thresholds, and Wilcoxon test was performed to compare differences in the relative frequencies of various type of immune cells.

Analysis of High- and Low-RS Groups in Terms of their Clinical Information

After establishing the prediction model, we conducted the difference analysis to present the diversity of the clinical information between the high- and low-RS groups through R packages “ggpubr” and “limma”. Next, Kaplan-Meier survival analysis was used to compare the survival outcome of patients between the two groups, and performed the difference analysis of TMB levels and immune checkpoint genes (ICGs) expression levels between high- and low-RS groups to predict their treatment responses when receiving immunotherapy.

Statistical Analysis

Statistical analyses were performed using R language (version 4.0.3). Kaplan-Meier survival analysis was used to assess the differences of the survival outcome between different groups. Univariate and multivariate Cox regression analyses were used to identify independent prognostic factors. p value < 0.05 was considered statistically significant in all tests.

Results

Identification of DEGs and IRDEGs

We first used R software “limma” package⁴ with false discovery rate (FDR) < 0.05 and |log2FC| >1 to identify DEGs between normal and THCA tissues, resulting in 2567 DEGs, as displayed in the volcano plot (Figure 2A). Then 294 IRDEGs were obtained from the intersection analysis of 567 DEGs and 1793 IRGs from ImmPort database and displayed in a Venn diagram (Figure 2B). expression levels of 294 IRDEGs between tumor and normal samples were visualized via a heatmap (Figure 2C).

FIGURE 2

FIGURE 2. Identification of immune-related DEGs in THCA. (A) Volcano plot for DEGs between THCA and normal samples from TCGA database. (B) Intersection of 2567 DEGs and 1793 IRGs in Venn plot (C) Heatmap showing immune-related IRDEGs. DEGs, differentially expressed genes; IRGs, immune-related genes.

Functional Enrichment Analysis for IRDEGs

To further investigate relevant pathway, biological process, cellular component, and molecular function of 294 IRDEGs, we performed GO and KEGG pathway enrichment analysis for these IRDEGs, as shown in Figure 3. GO enrichment analysis indicated that these IRDEGs were mostly involved in chemotaxis-related activities, cell adhesion regulation and other immune-related responses (Figure 3A). Besides, KEGG pathway analysis revealed that these IRDEGs were chiefly enriched in immune-related pathways, such as cytokine−cytokine receptor interaction, Viral protein with cytokine receptor, T cell receptor signaling pathway, Chemokine signaling pathway, Natural killer cell mediated cytotoxicity, PD−L1 expression and PD−1 checkpoint pathway in cancer (Figure 3B).

FIGURE 3

FIGURE 3. Enrichment analysis of IRDEGs. (A) GO enrichment analysis. (B) KEGG pathway enrichment analysis. IRDEGs, immune-related differentially expressed gene; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function.

Construction and Validation of THCA Prognosis Model

We first performed univariate Cox regression analysis on 294 IRDEGs to find the key IRDEGs that affect patients’ survival outcome with a threshold value of p < 0.001, which five genes (PPBP, RBP4, SEMA6B, VGF and GCGR) were selected as being significantly associated with THCA patients’ prognosis in this step (Table 1). Subsequently, LASSO regression was used for the following analyses (Figures 4A–B), thereby obtaining four candidate genes (PPBP, RBP4, SEMA6B, and GCGR). Eventually, three prognostic hub genes (PPBP, SEMA6B and GCGR) were screened from five these candidate genes by multivariate Cox regression analysis (Table 1). Next, we randomly assigned all THCA cases into training and test cohorts at a 2: 1 ratio for cross-validation, referred to as the training cohort and the test cohort, containing 2/3 and 1/3 of THCA cases, respectively. Additionally, no significant difference was observed in the clinical characteristics between the training and test cohorts (p > 0.01, Table 2). The THCA prognostic model was constructed based on the three hub genes and the RS of each patient in the prognostic model was calculated using the following formula: RS = (1.167391 × expression of PPBP) + (0.900831 × expression of SEMA6B) + (0.471683 × expression of GCGR), thereby dividing the patients into low- and high-RS groups according to their median value of RS. The heatmap for THCA tissues in the combined, test, and training three sets (Figures 5A–C) revealed that expression levels of three prognostic hub genes were downregulated in the low-RS group. RS distribution (Figures 5D–F) and RS-related survival status among patients (Figures 5G–I) indicated that higher RS corresponded to higher mortality risk in THCA patients.

TABLE 1

TABLE 1. Cox regression analysis for screening of IRDEGs influencing the THCA patients’ prognosis.

FIGURE 4

FIGURE 4. Identification of IRDEGs associated with THCA prognosis. (A–B) LASSO coefficient profiles of prognostic-related IRDEGs. LASSO, least absolute shrinkage and selection operator. IRDEGs, immune-related differentially expressed genes.

TABLE 2

TABLE 2. Clinical characteristics between test and training cohorts.

FIGURE 5

FIGURE 5. Construction of THCA prognostic model (A–C) Expression levels of three prognostic hub IRDEGs among THCA patients in the combined test, and training cohorts, respectively. (D–F) The distribution of RS among THCA patients in three cohorts. (G–I) The survival status correlated with RS among THCA patients. THCA, thyroid carcinoma; RS, risk score.

To further validate the reliability of THCA prognostic model, we further performed survival and ROC curves analyses. Kaplan–Meier survival plots for THCA patients in the combined, test and training three sets indicated that patients in the high-RS group had a substantially worse survival outcome than those in the low-RS group (Figures 6A–C), with p < 0.001, p = 0.023, and p = 0.006, respectively. Besides, the combined set showed that area under the curve (AUC) values for 1-, 3-, 5-yearsurvival were 0.847, 0.722, and 0.781, respectively (Figure 6D). Similarly, each AUC for 1-, 3-, 5-years survival was more than 0.6 in the test and training sets (Figures 6E, F). In brief, these findings demonstrated predictive accuracy and stability of the constructed model.

FIGURE 6

FIGURE 6. Cross-validation of a prognostic model in THCA. (A–C) OS analysis of THCA patients with high- and low-RS in the combined test, and training cohorts, with p < 0.001, p = 0.023, and p = 0.006, respectively. (D–F) ROC curve analysis for 1−, 3−, 5-years validation of RS model in three different cohorts. THCA, thyroid carcinoma; OS, overall survival; RS, risk score; ROC, receiver operating characteristic.

Identification of Independent Risk Factors for Prognosis in THCA Patients

We conducted univariate and multivariate Cox regression analyses to determine the correlation between RS and clinical features in THCA patients. The results revealed that prognosis of THCA patients was linked to three significant risk factors, including the age at diagnosis, pathological stage and RS (Table 3). Boxplots of the relationship between patients’ clinical characteristics and RS revealed that patients aged over 65 years old (Figure 7A), had advanced pathological stages (Figure 7C), or had a higher T stage (Figure 7D) showed a higher RS, with p = 0.003, p = 0.014, and p = 0.002, respectively. However, no significant difference was observed in the correlations of RS with gender (Figure 7B), N (Figure 7E) and M (Figure 7F)stages, with p > 0.05. These results revealed that RS was associated with progression and development of THCA patients.

TABLE 3

TABLE 3. Cox regression analysis for clinical characteristics and RS influencing prognosis in THCA patients.

FIGURE 7

FIGURE 7. Correlation between RS and clinical characteristics in THCA. (A) age > 65 years old and < 65 years old (B) female and male (C) Stage I-II and stage III−IV, (D) T stage (E) Nstage (F) Mstage. Age over 65 years old, higher pathological stage, and higher T stage correlated with higher RS, with p = 0.003, p = 0.014, and p = 0.002, respectively. RS, risk score; THCA, thyroid carcinoma.

Associations of TMB With Prognosis and RS

To examine the value of TMB in THCA prognosis, we analyzed the associations between TMB and overall survival (OS) time or RS. The Kaplan–Meier survival plot revealed that THCA patients with high TMB exhibited shorter OS time than those with low TMB (Figure 8A, p = 0.033). Besides, THCA patients in the high-RS group presented higher TMB levels than those in the low-RS group (Figure 8B, p = 0.0041). Such results indicated a negative relationship between THCA prognosis model and TMB. Waterfall plots displayed a landscape of the top 20 somatic gene mutations in 475 THCA samples from TCGA database, with different colors signifying different mutation types (Figures 8C,D).

FIGURE 8

FIGURE 8. Association of TMB with OS and RS (A) Survival analysis of THCA patients with high- and low-TMB (B) Comparisons of TMB in low- and high-RS groups (C) (D) Waterfall plot for mutation profiles of the top 20 genes in THCA samples of high- and low-RS groups, respectively. Annotations with different colors at the bottom referred to the various mutation types and bar chart above presented mutation burden. The right showed name of mutated genes and the right displayed percent of gene mutation. TMB, tumor mutation burden; OS, overall survival; RS, risk score.

Evaluation of TICs and Immune Checkpoint Genes Expression

To investigate the effect of constructed prognostic model on TICs we used CIBERSORT algorithm to calculate the distribution of 22 TIC subtypes in THCA samples. The Wilcoxon rank-sum test was utilized to compare the proportions of different TIC subtypes in high- and low-RS groups, and the findings were then visualized using a violin plot (Figure 9A). The results revealed that the relative abundance of monocytes, M2 macrophages, activated dendritic cells was significantly higher in high-RS group than those in low-RS group (p = 0.017, p = 0.042 and p < 0.001, respectively), whereas M1 macrophages and dendritic resting cells in high-RS group exhibited lower relative abundance (p = 0.036 and p = 0.002, respectively). According to their median TICs values, we classified THCA patients into high- and low- RS groups and then performed Kaplan–Meier survival analysis. The results displayed that high memory B cells was associated with poor OS (p = 0.002, Figure 9B).

FIGURE 9

FIGURE 9. TICs profile of THCA cases (A) Violin plot for quantification of 22 TIC subtypes between low- and high- RS groups. (B) Memory B cells influencing survival outcome of THCA patients. TIC, tumor-infiltrating immune cell; THCA, thyroid carcinoma; RS, risk score.

To further evaluate the effectiveness of immunotherapy on THCA, we analyzed expression levels of immunotherapy targets (PD-L1, PD-1 and CTLA4) in the high- and low-RS groups (Figure 10). The results revealed that the low-risk group significantly increased PD−L1 expression levels compared with the high-risk group (p = 0.044, Figure 10A), implying that low-RS THCA patients might be more susceptible to immunotherapy.

FIGURE 10

FIGURE 10. Expression levels of ICGs in high- and low-RS groups (A) PD-L1 expression of THCA patients in high- and low-RS groups. (B) PD-1 expression of THCA patients in high- and low-RS groups. (C) CTLA4 expression of THCA patients in high- and low-RS groups. RS, risk score; PD-L1, programmed cell death ligand 1; THCA, thyroid carcinoma; PD-1, programmed death 1; CTLA4, cytotoxic T lymphocyte-associated protein 4.

Discussion

THCA is the fifth most prevalent cancer in the United States, with an annual incidence increase of ∼5% (Varricchi et al., 2019). Although most THCA patients have a favorable prognosis, approximately 15–20% of DTC patients, most ATC patients display resistance to standard treatment methods, such as RAI. Sorafenib and lenvatinib, two recently approved multikinase inhibitors (MKIs), have demonstrated encouraging outcomes in progressive, RAI-refractory DTC. However, adverse consequences have been identified, restricting their use (Cabanillas and Habra., 2016; Krajewska et al., 2017). Notably, the development of ICIs, such as anti–CTLA-4 and anti–PD-1 agents, have revolutionized THCA treatment (Antonelli et al., 2018; Varricchi et al., 2019). However, not all patients benefit from them due to individual variances. In the current study, we constructed a prognostic prediction model by screening IRDEGs to predict immunotherapy efficacy and survival outcome of THCA.

In the present study, we first identified 294 IRDEGs between THCA and normal samples and then explored potential functional enrichment pathways of these IRDEGs via KEGG and GO functional enrichment analyses. We discovered that the IRDEGs were primally enriched in some immune-related pathways, suggesting that IRDEGs might play a critical role in alteration of THCA immune microenvironment. Subsequently, three hub IRDEGs associated with THCA prognosis were identified, including PPBP, SEMA6B and GCGR. PPBP, alternatively known as C-X-C chemokine ligand 7, stimulates various cellular processes, such as DNA synthesis, glucose metabolism, intracellular cAMP accumulation and fibrinogen activator synthesis. Some studies have indicated that PPBP and its encoded proteins might be linked to the progression of Wilms tumor (Guo et al., 2017) and gastric cancer (Yamamoto et al., 2019). SEMA6B a member of the semaphorin family, is mainly involved in developing peripheral and central nervous systems (Andermatt et al., 2014). SEMA6B has been demonstrated to have a critical role in the prognosis of various tumors, such as gliomas (Sun et al., 2019), breast cancer (Murad et al., 2006) and testicular cancer (Ji et al., 2020). A recent study has shown that SEMA6B promoted occurrence and development of THCA via regulating Notch signaling pathway (Lv et al., 2021). GCGR, a member of G protein-coupled receptor family, is critical in regulating glucose homeostasis. Previous studies indicated that GCGR aberrant expression might lead to adverse survival in endometrial stromal sarcoma (Davidson et al., 2013), renal cell carcinoma (Schrödter et al., 2016), and non-small cell lung cancer (NSCLC) (Li et al., 2020). In addition, it has been revealed that GCGR overexpression results in poor survival of PTC by activating epithelial-mesenchymal transition (EMT) and P38/ERK pathway (Jiang et al., 2020). Although we speculated that these genes might be a potential therapeutic target and/or prognostic biomarker for treating THCA patients, this hypothesis still requires additional validation in future studies.

Next, we firstly built a THCA prognostic model using the three screened hub genes (PPBP, SEMA6B and GCGR), and validated the reliability of the prognostic model in THCA. The results in the combined, test, and training sets revealed that THCA patients in the high-RS group exhibited poor survival outcomes than those in the low-RS group, and AUC values were over 0.6, implying that the constructed prognostic model accurately predicted OS in THCA. Besides, we found that older patients in the more advanced stage had significantly greater RS levels than younger ones in the earlier stage, consistent with the conclusion drawn by Ibrahimpasic et al. (Ibrahimpasic et al., 2019).

TMB, as a novel biomarker, has been recently known to forecast the clinical efficacy of immunotherapy in many cancers (Samstein et al., 2019), since it is closely associated with tumor immune infiltration and microenvironment alteration (Kang et al., 2020). Zhou et al. (Zhou et al., 2021) discovered that patients in the low-TMB group exhibited a worse survival outcome and immune response than those in the high-TMB group in. However, few studies focused on the predictive value of TMB for immunotherapy in THCA patients. In this current study we revealed that THCA patients in the high-RS group possessed higher TMB than those in the low-RS group. Herein, we speculated that higher TMB might be correlated with a worse prognosis in THCA, consistent with the conclusion drawn by Zhang et al. (Zhang et al., 2019a) in their work on clear cell renal cell carcinoma. Besides considering the important roles of TICs in TME on prognosis of numerous malignancies (Zhang et al., 2003; Tran Janco et al., 2015), we further assessed the distribution of 22 TIC subtypes in THCA samples from high- and low-RS groups and the relationship with survival outcome. The analysis results revealed that patients in the high-RS group chiefly possessed a higher level of monocytes, M2 macrophages and activated dendritic cells (DCs). In contrast, those in the low-RS group exhibited higher proportions of M1 macrophages and resting DCs, manifesting that tumor-associated macrophages were associated with tumor progression (Jackaman et al., 2017). Nevertheless, the underlying mechanism remains yet unclear. Similarly, Travers et al. (Travers et al., 2019) found that increasing tumor-killing M1 macrophages and decreasing M2 macrophages in TME contributed to reduced TMB and improved survival in mice with ovarian cancer. Additionally, a study indicated a strong correlation between DCs and advanced patients with PTC (Bergdorf et al., 2019), which might explain why THCA samples in the high-RS group exhibited advanced pathological stages. Besides, it has been reported that tumor-infiltrating memory B cells are linked to superior clinical outcomes in breast cancer (Garaud et al., 2019). In contrast, in the present study, we revealed that a high percent of memory B cells was highly correlated with poor survival in THCA.

Immune checkpoint molecules, such as PD-L1, PD-1 and CTLA4, have been demonstrated to connect with the efficacy of immunotherapy (Kalbasi and Ribas., 2020) Hence, we further investigated these biomarkers, including PD-L1, PD-1 and CTLA4, expression levels of THCA patients between in the high- and low-RS group. Results showed that PD-L1 expression was significantly upregulated in the low-RS group compared to the high-RS group, suggesting that this prognostic model might have ability to determine THCA patients’ response to immunotherapy. Based on these findings, we supposed that high TMB and low PD-L1 expression in THCA patients might respond poorly to immunotherapy.

Conclusions

In conclusion, this study constructed and validated a THCA prognostic prediction model based on TCGA database, displaying good predictability and reliability for THCA prognosis. To our knowledge, our group was the first to screen out three potential therapeutic target genes and elucidate the association of TICs with RS and OS in THCA. Additionally, we figured out that THCA patients in the high-RS group had high TMB and low PD-L1 expression, establishing a baseline and reference for predicting THCA immunotherapy efficacy in clinical trials. However, future investigations require relevant basic experiments and large sample clinical trials.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Author Contributions

SZ and SC wrote the manuscript. SC, YW, YZ, and JL analyzed data. BG and XN were responsible for data acquisition and interpretation SZ and SC revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the project of improving scientific ability of young college teachers in Guangxi (2021KY0110), Nanning Qingxiu District Science and Technology project (2020027), and Guangxi Health Commission project (Z20201043).

Conflict 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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Aboelnaga, E. M., and Ahmed, R. A. (2015). Difference between Papillary and Follicular Thyroid Carcinoma Outcomes: an Experience from Egyptian Institution. Cancer Biol. Med. 12 (1), 53–59. doi:10.7497/j.issn.2095-3941.2015.0005

PubMed Abstract | CrossRef Full Text | Google Scholar

Amaria, R. N., Reddy, S. M., Tawbi, H. A., Davies, M. A., Ross, M. I., Glitza, I. C., et al. (2018). Neoadjuvant Immune Checkpoint Blockade in High-Risk Resectable Melanoma. Nat. Med. 24 (11), 1649–1654. doi:10.1038/s41591-018-0197-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Anagnostou, V., Smith, K. N., Forde, P. M., Niknafs, N., Bhattacharya, R., White, J., et al. (2017). Evolution of Neoantigen Landscape during Immune Checkpoint Blockade in Non-small Cell Lung Cancer. Cancer Discov. 7 (3), 264–276. doi:10.1158/2159-8290.cd-16-0828

PubMed Abstract | CrossRef Full Text | Google Scholar

Andermatt, I., Wilson, N. H., Bergmann, T., Mauti, O., Gesemann, M., Sockanathan, S., et al. (2014). Semaphorin 6B Acts as a Receptor in post-crossing Commissural Axon Guidance. Development 141 (19), 3709–3720. doi:10.1242/dev.112185

PubMed Abstract | CrossRef Full Text | Google Scholar

Antonelli, A., Ferrari, S. M., and Fallahi, P. (2018). Current and Future Immunotherapies for Thyroid Cancer. Expert Rev. Anticancer Ther. 18 (2), 149–159. doi:10.1080/14737140.2018.1417845

PubMed Abstract | CrossRef Full Text | Google Scholar

Arndt, A., Steinestel, K., Rump, A., Sroya, M., Bogdanova, T., Kovgan, L., et al. (2018). Anaplastic Lymphoma Kinase (ALK ) Gene Rearrangements in Radiation-Related Human Papillary Thyroid Carcinoma after the Chernobyl Accident. J. Path: Clin. Res. 4 (3), 175–183. doi:10.1002/cjp2.102

CrossRef Full Text | Google Scholar

Barrios, D. M., Do, M. H., Phillips, G. S., Postow, M. A., Akaike, T., Nghiem, P., et al. (2020). Immune Checkpoint Inhibitors to Treat Cutaneous Malignancies. J. Am. Acad. Dermatol. 83 (5), 1239–1253. doi:10.1016/j.jaad.2020.03.131

CrossRef Full Text | Google Scholar

Bergdorf, K., Ferguson, D. C., Mehrad, M., Ely, K., Stricker, T., and Weiss, V. L. (2019). Papillary Thyroid Carcinoma Behavior: Clues in the Tumor Microenvironment. Endocrine-Related Cancer 26 (6), 601–614. doi:10.1530/erc-19-0074

PubMed Abstract | CrossRef Full Text | Google Scholar

Bissell, M. J., and Hines, W. C. (2011). Why Don't We Get More Cancer? A Proposed Role of the Microenvironment in Restraining Cancer Progression. Nat. Med. 17 (3), 320–329. doi:10.1038/nm.2328

PubMed Abstract | CrossRef Full Text | Google Scholar

Bos, J. L. (1989). Ras Oncogenes in Human Cancer: a Review. Cancer Res. 49 (17), 4682–4689.

PubMed Abstract | Google Scholar

Cabanillas, M. E., and Habra, M. A. (2016). Lenvatinib: Role in Thyroid Cancer and Other Solid Tumors. Cancer Treat. Rev. 42, 47–55. doi:10.1016/j.ctrv.2015.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabanillas, M. E., McFadden, D. G., and Durante, C. (2016). Thyroid Cancer. The Lancet 388 (10061), 2783–2795. doi:10.1016/s0140-6736(16)30172-6

CrossRef Full Text | Google Scholar

Cancer Genome Atlas Research Network (2014). Integrated Genomic Characterization of Papillary Thyroid Carcinoma. Cell 159 (3), 676–690. doi:10.1016/j.cell.2014.09.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, B., Abeler, V. M., Hellesylt, E., Holth, A., Shih, I.-M., Skeie-Jensen, T., et al. (2013). Gene Expression Signatures Differentiate Uterine Endometrial Stromal Sarcoma from Leiomyosarcoma. Gynecol. Oncol. 128 (2), 349–355. doi:10.1016/j.ygyno.2012.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferris, R. L., Blumenschein, G., Fayette, J., Guigay, J., Colevas, A. D., Licitra, L., et al. (2018). Nivolumab vs Investigator's Choice in Recurrent or Metastatic Squamous Cell Carcinoma of the Head and Neck: 2-year Long-Term Survival Update of CheckMate 141 with Analyses by Tumor PD-L1 Expression. Oral Oncol. 81, 45–51. doi:10.1016/j.oraloncology.2018.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferris, R. L., Blumenschein, G., Fayette, J., Guigay, J., Colevas, A. D., Licitra, L., et al. (2016). Nivolumab for Recurrent Squamous-Cell Carcinoma of the Head and Neck. N. Engl. J. Med. 375 (19), 1856–1867. doi:10.1056/NEJMoa1602252

CrossRef Full Text | Google Scholar

French, J. D., Kotnis, G. R., Said, S., Raeburn, C. D., McIntyre, R. C., Klopper, J. P., et al. (2012). Programmed Death-1+ T Cells and Regulatory T Cells Are Enriched in Tumor-Involved Lymph Nodes and Associated with Aggressive Features in Papillary Thyroid Cancer. J. Clin. Endocrinol. Metab. 97 (6), E934–E943. doi:10.1210/jc.2011-3428

CrossRef Full Text | Google Scholar

Garaud, S., Buisseret, L., Solinas, C., Gu-Trantien, C., de Wind, A., Van den Eynden, G., et al. (2019). Tumor-infiltrating B Cells Signal Functional Humoral Immune Responses in Breast Cancer. JCI Insight 4 (18). doi:10.1172/jci.insight.129641

PubMed Abstract | CrossRef Full Text | Google Scholar

Giannone, G., Ghisoni, E., Genta, S., Scotto, G., Tuninetti, V., Turinetto, M., et al. (2020). Immuno-Metabolism and Microenvironment in Cancer: Key Players for Immunotherapy. Ijms 21 (12), 4414. doi:10.3390/ijms21124414

PubMed Abstract | CrossRef Full Text | Google Scholar

Gunda, V., Gigliotti, B., Ashry, T., Ndishabandi, D., McCarthy, M., Zhou, Z., et al. (2019). Anti-PD-1/PD-L1 Therapy Augments Lenvatinib's Efficacy by Favorably Altering the Immune Microenvironment of Murine Anaplastic Thyroid Cancer. Int. J. Cancer 144 (9), 2266–2278. doi:10.1002/ijc.32041

CrossRef Full Text | Google Scholar

Guo, F., Zhang, J., Wang, L., Zhao, W., Yu, J., Zheng, S., et al. (2017). Identification of Differentially Expressed Inflammatory Factors in Wilms Tumors and Their Association with Patient Outcomes. Oncol. Lett. 14 (1), 687–694. doi:10.3892/ol.2017.6261

PubMed Abstract | CrossRef Full Text | Google Scholar

Ibrahimpasic, T., Ghossein, R., Shah, J. P., and Ganly, I. (2019). Poorly Differentiated Carcinoma of the Thyroid Gland: Current Status and Future Prospects. Thyroid 29 (3), 311–321. doi:10.1089/thy.2018.0509

PubMed Abstract | CrossRef Full Text | Google Scholar

Jackaman, C., Tomay, F., Duong, L., Razak, N. B. A., Pixley, F. J., Metharom, P., et al. (2017). Aging and Cancer: The Role of Macrophages and Neutrophils. Ageing Res. Rev. 36, 105–116. doi:10.1016/j.arr.2017.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, C., Wang, Y., Wang, Y., Luan, J., Yao, L., Wang, Y., et al. (2020). Immune-related Genes Play an Important Role in the Prognosis of Patients with Testicular Germ Cell Tumor. Ann. Transl Med. 8 (14), 866. doi:10.21037/atm-20-654

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, H.-C., Chen, X.-R., Sun, H.-F., and Nie, Y.-W. (2020). Tumor Promoting Effects of Glucagon Receptor: a Promising Biomarker of Papillary Thyroid Carcinoma via Regulating EMT and P38/ERK Pathways. Hum. Cell 33 (1), 175–184. doi:10.1007/s13577-019-00284-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalbasi, A., and Ribas, A. (2020). Tumour-intrinsic Resistance to Immune Checkpoint Blockade. Nat. Rev. Immunol. 20 (1), 25–39. doi:10.1038/s41577-019-0218-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, K., Xie, F., Mao, J., Bai, Y., and Wang, X. (2020). Significance of Tumor Mutation Burden in Immune Infiltration and Prognosis in Cutaneous Melanoma. Front. Oncol. 10, 573141. doi:10.3389/fonc.2020.573141

PubMed Abstract | CrossRef Full Text | Google Scholar

Krajewska, J., Gawlik, T., and Jarzab, B. (2017). Advances in Small Molecule Therapy for Treating Metastatic Thyroid Cancer. Expert Opin. Pharmacother. 18 (11), 1049–1060. doi:10.1080/14656566.2017.1340939

PubMed Abstract | CrossRef Full Text | Google Scholar

Laha, D., Nilubol, N., and Boufraqech, M. (2020). New Therapies for Advanced Thyroid Cancer. Front. Endocrinol. 11, 82. doi:10.3389/fendo.2020.00082

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Liu, X., Zhou, X. J., Chen, X., Li, J. P., Yin, Y. H., et al. (2020). Identification and Validation of the Prognostic Value of Immune-Related Genes in Non-small Cell Lung Cancer. Am. J. Transl Res. 12 (9), 5844–5865.

Google Scholar

Lv, X. J., Chen, X., Wang, Y., Yu, S., Pang, L., and Huang, C. (2021). Aberrant Expression of Semaphorin 6B Affects Cell Phenotypes in Thyroid Carcinoma by Activating the Notch Signalling Pathway. Endokrynol Pol. 72 (1), 29–36. doi:10.5603/EP.a2020.0072

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehnert, J. M., Varga, A., Brose, M. S., Aggarwal, R. R., Lin, C.-C., Prawira, A., et al. (2019). Safety and Antitumor Activity of the Anti-PD-1 Antibody Pembrolizumab in Patients with Advanced, PD-L1-Positive Papillary or Follicular Thyroid Cancer. BMC Cancer 19 (1), 196. doi:10.1186/s12885-019-5380-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Murad, H., Collet, P., Huin-Schohn, C., Al-Makdissy, N., Kerjan, G., Chedotal, A., et al. (2006). Effects of PPAR and RXR Ligands in Semaphorin 6B Gene Expression of Human MCF-7 Breast Cancer Cells. Int. J. Oncol. 28 (4), 977–984. doi:10.3892/ijo.28.4.977

CrossRef Full Text | Google Scholar

Nath, M. C., and Erickson, L. A. (2018). Aggressive Variants of Papillary Thyroid Carcinoma: Hobnail, Tall Cell, Columnar, and Solid. Adv. Anat. Pathol. 25 (3), 172–179. doi:10.1097/pap.0000000000000184

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12 (5), 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikiforov, Y. E., and Nikiforova, M. N. (2011). Molecular Genetics and Diagnosis of Thyroid Cancer. Nat. Rev. Endocrinol. 7 (10), 569–580. doi:10.1038/nrendo.2011.142

PubMed Abstract | CrossRef Full Text | Google Scholar

Oiseth, S. J., and Aziz, M. S. (2017). Cancer Immunotherapy: a Brief Review of the History, Possibilities, and Challenges Ahead. Jcmt 3, 250–261. doi:10.20517/2394-4722.2017.41

CrossRef Full Text | Google Scholar

Samstein, R. M., Lee, C.-H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor Mutational Load Predicts Survival after Immunotherapy across Multiple Cancer Types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Santa-Maria, C. A., and Nanda, R. (2018). Immune Checkpoint Inhibitor Therapy in Breast Cancer. J. Natl. Compr. Canc Netw. 16 (10), 1259–1268. doi:10.6004/jnccn.2018.7046

PubMed Abstract | CrossRef Full Text | Google Scholar

Schrödter, S., Braun, M., Syring, I., Klümper, N., Deng, M., Schmidt, D., et al. (2016). Identification of the Dopamine Transporter SLC6A3 as a Biomarker for Patients with Renal Cell Carcinoma. Mol. Cancer 15, 10. doi:10.1186/s12943-016-0495-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Liu, X., Xia, M., Shao, Y., and Zhang, X. D. (2019). Multicellular Gene Network Analysis Identifies a Macrophage-Related Gene Signature Predictive of Therapeutic Response and Prognosis of Gliomas. J. Transl Med. 17 (1), 159. doi:10.1186/s12967-019-1908-1

CrossRef Full Text | Google Scholar

Tran Janco, J. M., Lamichhane, P., Karyampudi, L., and Knutson, K. L. (2015). Tumor-infiltrating Dendritic Cells in Cancer Pathogenesis. J.I. 194 (7), 2985–2991. doi:10.4049/jimmunol.1403134

CrossRef Full Text | Google Scholar

Travers, M., Brown, S. M., Dunworth, M., Holbert, C. E., Wiehagen, K. R., Bachman, K. E., et al. (2019). DFMO and 5-Azacytidine Increase M1 Macrophages in the Tumor Microenvironment of Murine Ovarian Cancer. Cancer Res. 79 (13), 3445–3454. doi:10.1158/0008-5472.can-18-4018

PubMed Abstract | CrossRef Full Text | Google Scholar

Varricchi, G., Loffredo, S., Marone, G., Modestino, L., Fallahi, P., Ferrari, S. M., et al. (2019). The Immune Landscape of Thyroid Cancer in the Context of Immune Checkpoint Inhibition. Ijms 20 (16), 3934. doi:10.3390/ijms20163934

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamamoto, Y., Kuroda, K., Sera, T., Sugimoto, A., Kushiyama, S., Nishimura, S., et al. (2019). The Clinicopathological Significance of the CXCR2 Ligands, CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, CXCL7, and CXCL8 in Gastric Cancer. Anticancer Res. 39 (12), 6645–6652. doi:10.21873/anticanres.13879

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, Y., Yan, Z., Lian, S., Wei, L., Zhou, C., Feng, D., et al. (2020). Prognostic Value of Novel Immune-Related Genomic Biomarkers Identified in Head and Neck Squamous Cell Carcinoma. J. Immunother. Cancer 8 (2), e000444. doi:10.1136/jitc-2019-000444

CrossRef Full Text | Google Scholar

Zhang, C., Li, Z., Qi, F., Hu, X., and Luo, J. (2019a). Exploration of the Relationships between Tumor Mutation burden with Immune Infiltrates in clear Cell Renal Cell Carcinoma. Ann. Transl. Med. 7 (22), 648. doi:10.21037/atm.2019.10.84

CrossRef Full Text | Google Scholar

Zhang, L., Conejo-Garcia, J. R., Katsaros, D., Gimotty, P. A., Massobrio, M., Regnani, G., et al. (2003). Intratumoral T Cells, Recurrence, and Survival in Epithelial Ovarian Cancer. N. Engl. J. Med. 348 (3), 203–213. doi:10.1056/NEJMoa020177

CrossRef Full Text | Google Scholar

Zhang, W., Sun, W., Qin, Y., Wu, C., He, L., Zhang, T., et al. (2019b). Knockdown of KDM1A Suppresses Tumour Migration and Invasion by Epigenetically Regulating the TIMP1/MMP9 Pathway in Papillary Thyroid Cancer. J. Cell Mol Med 23 (8), 4933–4944. doi:10.1111/jcmm.14311

CrossRef Full Text | Google Scholar

Zhou, H., Chen, L., Lei, Y., Li, T., Li, H., and Cheng, X. (2021). Integrated Analysis of Tumor Mutation burden and Immune Infiltrates in Endometrial Cancer. Curr. Probl. Cancer 45 (2), 100660. doi:10.1016/j.currproblcancer.2020.100660

PubMed Abstract | CrossRef Full Text | Google Scholar