Characterization of the Immune Cell Infiltration Landscape of Thyroid Cancer for Improved Immunotherapy

Within the endocrine system, thyroid cancer (THCA) is the most typical malignant tumor. Tumor-infiltrating immune cells play vital roles in tumor progression, recurrence, metastasis as well as response to immunotherapy. However, THCA’s immune infiltrative landscape is still not clarified. Therefore, we utilized two statistical algorithms to investigate the immune cell infiltration (ICI) landscape of 505 THCA samples and defined three ICI immune subtypes. The ICI scores were calculated using principal-component analysis. Increased tumor mutation burden (TMB) and immune-related signaling pathways were associated to a high ICI score. The high ICI score group indicated a relatively longer overall survival (OS) than the low ICI score group. Most immune checkpoint-related and immune activation-related genes were considerably upregulated in the ICI high group, which indicates stronger immunogenicity and a greater likelihood of benefiting from immunotherapy. In two cohort studies of patients receiving immunotherapy, high-ICI-score group showed notable therapeutic effects and clinical advantages compared to those with lower ICI scores. These results demonstrate that ICI score acts as an effective prognostic indicator and predictor of response to immunotherapy.

Introduction

The most common malignant tumor of the endocrine system is thyroid cancer (THCA). As tumor detection has advanced, the global incidence of THCA has risen rapidly. Younger people and women are more likely to suffer from it (La Vecchia et al., 2015; Du et al., 2018; Siegel et al., 2020). There are four major pathological subtypes of thyroid cancer: papillary thyroid carcinoma, follicular thyroid carcinoma, medullary thyroid carcinoma, and undifferentiated carcinoma. The prognoses vary greatly (Cabanillas et al., 2016). Traditional treatments including chemotherapy and radiotherapy cannot improve the therapeutic effect of locally advanced or metastatic thyroid cancer, which makes the development of effective treatments crucial.

The tumor microenvironment (TME) of thyroid cancer is rich in immune cells, which indicates that it is an ideal candidate for immunotherapy. And certain immune properties have been shown to affect the prognosis of thyroid cancer (Ferrari et al., 2019). Extensive researches on the immune microenvironment of thyroid cancer have shown that components of an individual’s immune system were closely related to the occurrence, invasion, and metastasis as well as the therapeutic response to immunotherapy (Xie et al., 2020; Yin et al., 2020). A new clinical experiment showed that tumor immune cell infiltration (ICI) in thyroid cancer is related to sensitivity to immunotherapy and a better prognosis (Na and Choi, 2018). In papillary thyroid carcinoma, tumor-associated macrophages (TAMs) are correlated to lymph node metastasis, increased tumor size, and a diminished survival rate (Kim et al., 2013; Fang et al., 2014). Increased regulatory T cell (Treg) infiltration suggests a positive correlation with disease stage, whereas natural killer (NK) cell infiltration is adversely associated with disease stage (Gogali et al., 2012). Additionally, thyroid cancer cells can produce a variety of cytokines and chemokines, which can promote the tumor progression. Reducing the concentration of cytokines and chemokines in the tumor microenvironment produces therapeutic benefits (Coperchini et al., 2019). An improved understanding of the molecular and immunological properties of the THCA tumor microenvironment can be used to aid immunotherapy.

We collected The Cancer Genome Atlas (TCGA)-THCA datasets, evaluating the immune microenvironment and immune cell infiltration in THCA using two calculation tools, CIBERSORT, and ESTIMATE. According to the TCGA-THCA immunophenoscores (IPSs), we identified the relevant genes to establish an immune cell infiltration (ICI) scoring model and verified it using the TCGA-THCA dataset. In summary, this ICI score system can be an effective prognostic biomarker and predictor for the valuation of THCA immunotherapy response.

Methods

Obtaining Expression Profile Data and Clinical Data

The overall analysis ideas of this research are shown as follows. First, download THCA expression profile data as well as follow-up clinical evidences based on the TCGA database (https://portal.gdc.cancer.gov/). The RNA sequence data of TCGA-THCA is processed in the following steps: 1) Remove samples which did not include subsequent clinical information; 2) Remove samples with no survival period and survival status; 3) Convert the probe to Gene Symbol; 4) One probe corresponds to several genes, remove the probe; 5) Use the median value for the expression of multiple Gene Symbols. The pre-processed TCGA-THCA data has a total of 505 tumor samples, and the clinical statistics of the samples are shown in Table 1.

TABLE 1

TABLE 1. Clinical information of TCGA-THCA datasets.

Tumor Immunophenoscore (IPS Database)

Tumor immunophenoscore (IPS) comes from the website of The Cancer Immunome database (TCIA) (https://tcia.at/patients). It is in view of the characteristics of tumor-infiltrating immune cells, acting as a bridge between immune cell infiltrating subtypes and immune gene subtypes. (Charoentong et al., 2016).

Consensus Clustering Analysis of Immune Cell Infiltration

22 different immune cells infiltration level in THCA (naïve B cells, memory B cells, plasma cells, CD8⁺ T cells, naïve CD4⁺ T cells, inactivated and activated memory CD4⁺ T cells, T follicular helper cells, Tregs, gamma delta T cells, resting and activated NK cells, monocytes, M0, M1, and M2 macrophages, resting and activated dendritic cells, resting and activated mast cells, eosinophils, and neutrophils) was quantified utilizing the CIBERSORT R package, based on the LM22 signature and 1,000 permutations. We use ESTIMATE R package to estimate the degree of immune infiltration and matrix purity scores in each THCA sample. “Pam” method according to Euclid and Ward’s linkages was utilized, and executed by “ConsensuClusterPlus package” R package. Unsupervised clustering was performed, and researchers repeated the clustering 1,000 times to ensure the classification’s stability.

Differential Expression of Genes Associated With the Tumor IPS

The tumor samples were divided into two groups according to the optimal density gradient threshold of tumor immunophenoscore (IPS) related to survival with the method of R Survminer package. The differential expression of genes was analyzed between the TCGA-THCA tumor sample with different score groups through specifying the cutoff value 0.05 (adjusted) and | log2 (Fold Change) | >1 using the limma R package. Furthermore, the genes which are positively related to the consistent classification result are called immune genotype-related feature A, and the remaining genes are called for feature B.

Dimension Reduction of Gene Features and Construction of ICI Score

In order to measure the immune cell infiltration in tumors through gene expression, this study constructed a tumor immune cell infiltration score (ICI scores) model based on feature A and feature B gene sets related to immune gene subtypes. First of all, we conduct the Boruta algorithm to reduce dimension of the feature A and B genes. After reduction, ICI score A and B which were defined by ICI feature gene A and B were evaluated by the principal component analysis (PCA). Thirdly we construct an immune cell infiltration score (ICI scores) model, and the calculation formula is as follows (Sotiriou et al., 2006):

$\mathrm{I}CI\_scores={\displaystyle \sum PC1\left(A\right)}-{\displaystyle \sum PC1\left(B\right)}$

Collection of Somatic Alteration Data

Mutation data of TCGA-THCA were downloaded from the TCGA (https://www.cancer.gov/tcga/). We used the Survminer R package to calculate the optimal density gradient threshold related to the tumor mutation burden (TMB) and survival and categorized the samples into high and low tumor TMB groups. Implementing the maftool R package, the mutation frequencies of the top 30 driver genes in different immune cell infiltration score (ICI scores) groups were compared.

Acquisition of Immunotherapy Data Sets

To explore the effectiveness of ICI score in predicting the benefit of immunotherapy treatment, we downloaded the expression profile data and clinical materials of the IMvigor210 cohort (http://research-pub.gene.com/IMvigor210CoreBiologies/), the ICI scoring model is used to divide all samples into a high and a low scoring group. Similarly, the GSE78220 data set was downloaded from GEO and analyzed accordingly.

Statistical Analyses

All the statistical comparisons involved in this study and the hypothesis testing of the significance of differences between groups are based on the statistical analysis method of R 3.6.

Results

The Immune Cell Infiltration Landscape

For each sample in the TCGA-THCA dataset, 22 types of immune cell’s infiltration statuses (naïve B cells, memory B cells, plasma cells, CD8⁺ T cells, naïve CD4⁺ T cells, resting and activated memory CD4⁺ T cells, T follicular helper cells, Tregs, gamma delta T cells, resting and activated NK cells, monocytes, M0, M1, and M2 macrophages, resting and activated dendritic cells, resting and mast cells, eosinophils, and neutrophils) were recorded (Supplementary Table S1). From this, a cluster correlation and cumulative distribution function (CDF) curve were calculated (k = 3) and used for stable Immune Cell Infiltration (ICI) subtype classification (Figures 1A,B). Of the three main subtypes, ICI3 had the worst prognosis and a median OS of 937 days, whereas ICI1/2 had a mean OS of 1,017 days (p = 0.0065, Figure 1C).

FIGURE 1

FIGURE 1. Identification of tumor immune subtypes and characteristics of immune infiltration. (A) Clustering result when k = 3; (B) CDF curve distribution of consistent clustering; (C) Survival curve when K = 3; (D) Heat map of immune cell infiltration characteristics; (E) Box diagram of differences in immune cell infiltration characteristics among the three ICI subtypes; (F) The expression difference of PD1 among the three ICI subtypes; (G) the expression difference of PD-L1 among the three ICI subtypes.

To better clarify the inherent biological causes of the different clinical phenotypes, the TMEs of the three molecular subtypes were compared and visualized using heat maps (Figure 1D,E). The ICI3 samples were characterized by high stromal scores, low infiltration of CD8⁺ T cells, and greater infiltration of static memory CD4⁺ T cells and active dendritic cells. Subtypes ICI1 and ICI2, which had better prognoses, were characterized by lower stromal scores, greater CD8⁺ T cell infiltration, lower infiltration of Tregs, inactive memory CD4⁺ T cells and active dendritic cells.

In addition to analyzing immune cell infiltration, the expression levels of two important immune checkpoints (PD1 and PD-L1) in each ICI subtype were further analyzed. One of the characteristics of the ICI3 subtype was significantly lower PD1/PD-L1 expression, while the ICI2 subtype had higher PD1/PD-L1 expression. To assess the statistical significance of the differences in immune infiltration and PD1/PD-L1 expression among the ICI subtypes, the Kruskal-Wallis test and a hypothesis test were used. As shown in Figures 1F,G, the differences in PD1/PD-L1 expression and the infiltration of most immune cell types were statistically significant.

Immune Gene Expression Subtypes

In order to display different biological characteristics of immune phenotype, the tumor immunophenoscores (IPSs) samples in the TCGA-THCA dataset were collected from The Cancer Immunome database (TCIA), as shown in Supplementary Table S2. The Survminer R package was used to calculate the optimal density gradient threshold for the IPSs. A high and a low group of IPS with a threshold of seven were presented from TCGA-THCA tumor samples. The survival of high IPS group was much longer than the other group significantly (Supplementary Figure S1).

We analyzed differential gene expression (DEG) of different IPS group utilizing the limma R package with the set screening threshold changed (p < 0.05 and | log2 (Fold Change) | >1). One thousand two-hundred and ten differentially expressed genes (Supplementary Table S3) were identified, of which 1,030 genes were actively expressed in the high IPS group and 180 genes were greatly expressed in the low IPS group. We then performed unsupervised clustering of 1,210 differentially expressed genes related to IPS (IPS_DEGs) with the ConsensuClusterPlus R package. Finally, according to the cluster correlation and CDF curve, the TCGA-THCA tumor samples were divided into four immune gene expression subtypes, IPS_DEGs Clusters 1-4 (Figures 2A,B). There were clear differences in survival between the different immune gene expression subtypes (log rank test, p = 0.03, Figure 2C). Gene Type A includes 191 gene signatures which had a positive association with immune gene subtypes, while the remaining 1019 IPS_DEGs were represented by Gene Type B (Supplementary Table S4)

FIGURE 2

FIGURE 2. Consistent clustering of differentially expressed genes related to tumor immune antigenicity. (A) Clustering results when k = 4; (B) CDF curve distribution of consistent clustering; (C) Survival curve of molecular subtypes when k = 4; (D) Heat map of immune gene subtypes; (E) Functional enrichment of immune gene subtype related features A gene; (F) Functional enrichment of immune gene subtype related features B gene; (G) Differences in infiltration characteristics of immune cells between immune gene subtypes.

We revealed the heatmap of gene expression for each ICI subtype, IPS_DEGs Cluster, and Gene Type (Figure 2D) and implemented gene ontology (GO) functional enrichment study on Gene Type A and B genes using the clusterProfiler R package. A bubble map was used to display the first 10 pathways enriched in three groups of different functions (Biological Process, Cellular Component, Molecular Function), as shown in Figures 2E,F. Most of the enriched pathways were related to immunobiological processes.

We examined the 22 immune infiltration cells in the four IPS_DEG subtypes and found considerable differences between different subtypes (Figure 2G). We also observed great differences in PD1/PD-L1 expression. IPS_DEGs clusters 1, 2, and 4 were related to greater PD1/PD-L1 expression, while IPS_DEGs cluster 3 had the lowest PD1/PD-L1 expression (Supplementary Figure S2A, B). This low PD1/PD-L1 expression indicates a lack of sensitivity to immunotherapy, which is related to a poor prognosis. The consistent associations between immune cell infiltration characteristics and the prognosis of different gene expression subtypes indicate that the classification strategy for the immune cell subtype is scientific and reasonable.

ICI Score Construction

To measure the immune cell infiltration in thyroid tumors using differential expression genes, we constructed a tumor immune cell infiltration (ICI) score model based on the type A and type B gene sets associated with immune gene subtypes. The Survminer R package was utilized to calculate the most accurate density gradient threshold related to ICI scores and survival. The cutoff value was -27.33, which divided the TCGA-THCA tumor samples into high and low ICI score groups. The high ICI score group manifested longer OS than the low ICI score group in Figures 3A,B (p = 0.014). A Sangji diagram was adopted to visualize the relations between the immune gene expression subtypes, ICI score groups, and survival statuses (Figure 3C).

FIGURE 3

FIGURE 3. Grouping and model construction of tumor immune cell infiltration scores (ICI scores). (A) Distribution (above) and the selection of the best density gradient points (below) of immune cell infiltration score (ICI scores); (B) survival curve between high and low ICI scores; (C) Sankey diagram of immune gene subtypes, ICI scores grouping and survival status; (D) Differences in immune checkpoints and immune activation gene expression between ICI scores groups; (E) Gene set function enrichment of differentially expressed genes between ICI scores groups; (F) Survival curve between radiotherapy and ICI scores groups; (G) Survival curve received therapy in different ICI scores.

Before the prognostic value of ICI score could be assessed for the TCGA-THCA cohort, it was necessary to analyze the immune activity and tolerance conditions within different groups. We chose CD274, CTLA4, HAVCR2, IDO1, LAG3, and PDCD1 as immune checkpoint features, and CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF as immune activation-related genes (Zhang et al., 2020a). Except for IDO1, TBX2, and TNF, expression of most genes related to immune checkpoints and immune activity was significantly increased in the high ICI score group (Figure 3D). We conducted gene set enrichment analysis (GSEA) to make an in-depth exploration of the biological differences between the two ICI score groups. The top ten enriched pathways were cell adhesion molecules (CAMs), chemokine signaling pathways, complement and coagulation cascades, cytokine receptor interactions, hematopoietic cell lineage, leishmania infection, NK cell mediated cytotoxicity, T cell receptor signaling, and viral myocarditis (Figure 3E). We also evaluated the impact of radiotherapy and targeted therapy on the prognosis of each ICI subgroup. The patients of high ICI score group had a better survival compared with the other group no matter what kind of treatment they received (Figures 3F,G).

Somatic Variants in Different ICI Groups

A large body of data suggests that the tumor mutation burden (TMB) may forecast the response to immunotherapy. Thus, researching on the relationship between TMB and ICI score to clarify the genetic characteristics is crucial. We first used the Survminer R package to calculate the optimal density gradient threshold for TMB and survival. A score of 0.42 was the threshold that divided the TCGA-THCA tumor samples into two groups: high TMB group and low TMB group. There were significant differences in survival between the two groups, as shown in Supplementary Figure S3.

We first conducted a comparison of the TMB in the high and low ICI score groups. According to Figure 4A, TMB were significantly higher in the high ICI score group (Wilcoxon test, p = 0.024). Further correlation analysis showed a significant positive relationship between ICI score and TMB score (Spearman coefficient: R = 0.018, p = 0.047; Figure 4B). However, the prognosis of patients with high TMB performed worse prognosis compared with low TMB patients, which was conflicted with the prognoses based on ICI scores. Therefore, in the prognostic stratification of THCA, we further analyzed the synergy of TMB and ICI score. According to the TMB, the stratified survival analysis showed that either in the high or low TMB subgroups, the prognosis of patients with the high ICI score were always better than the low ICI score group (p <0.05; Figure 4C). This indicates that TMB status does not impede predictive effects of ICI scores. Accordingly, the above-mentioned results demonstrated that the ICI score may be a potential predictor of patient prognosis independent of TMB and that it is an effective predictor of response to immunotherapy.

FIGURE 4

FIGURE 4. Correlation between ICI scores and Somatic Variants. (A) Tumor mutation burden (TMB) between high and low ICI scores groups; (B) Association between ICI scores and tumor mutation burden (TMB); (C) Survival curve between tumor mutation burden (TMB) and ICI scores groups; (D) waterfall chart of gene mutation distribution in tumors of high ICI scores group; (E) waterfall chart of gene mutation distribution in tumors of low ICI scores group.

In addition to overall TMB, the distribution of somatic variation in THCA driver genes between different ICI subgroups was further analyzed. And the 30 drive genes with the most significant change frequencies were compared (Figures 4D,E). We analyzed the mutation annotation files from the TCGA-THCA data and found huge variations in the mutation profiles of the high and low ICI subgroups. These findings may provide new insight for studying the tumor ICI composition and understanding mutational mechanisms affecting immune checkpoint inhibitor therapy.

The Value of ICI Scores in Predicting Response to Immunotherapy

Immune checkpoint blockade therapy has emerged to treat cancer by blocking T cell inhibition pathways. To explore the relationship between ICI scores and response to immunotherapy, we analyzed expression profile data and clinical materials from the IMvigor210 cohort. All samples were defined within high and low ICI score groups using the ICI scoring model. Our results showed that in the IMvigor210 cohort, high ICI scores were associated with greater objective response to anti-PD-L1 treatment (Wilcoxon test, p = 0.036, Figure 5A) and longer survival (log rank test, p = 0.046, Figure 5B). Patients in the high ICI score group had higher objective responses to anti-PD-L1 treatment than (Wilcoxon test, p < 0.01; Figure 5C). The GSE78220 cohort in GEO, which received different types of immunotherapies, including cytokines, vaccines, and checkpoint blockers, showed similar findings (Wilcoxon test, p = 0.017; Figure 5D; log rank test, p = 0.023, Figure 5E; chi-square test, p0.01, Figure 5F). Overall, these findings pointed out that the ICI score may predict response to immunotherapy.

FIGURE 5

FIGURE 5. The role of tumor immune cell infiltration scores in the prediction of immunotherapy. (A) The difference of ICI scores between different treatment response groups in the IMvigor210 cohort; (B) the survival curve between the high and low ICI scores groups in the IMvigor210 cohort; (C) the difference of response between different ICI scores in the IMvigor210 cohort; (D) Difference of ICI scores between different treatment response groups in the GSE78220 cohort; (E) survival curve between different ICI score groups in the GSE78220 cohort; (F) difference of proportion of response between different ICI score groups in the GSE78220 cohort.

Discussion

Immunotherapy has shown amazing efficacy against a variety of malignant tumors, but for thyroid cancer it is still in exploratory stages. An early clinical trial showed that the PD-1 inhibitor pembrolizumab may be safe and effective in treating PD-L1-positive thyroid cancer. Studies have also shown a positive therapeutic effect of immunotherapy for patients with undifferentiated thyroid cancer (Capdevila et al., 2020; Ma et al., 2020). Combination therapies, including immune checkpoint inhibitors and tyrosine kinase (TK) or serine/threonine protein kinase B-Raf (BRAF) inhibitors, have shown significant potential for the treatment of advanced thyroid cancer (Iyer et al., 2018; Varricchi et al., 2019). Finding predictive indicators for evaluating the efficacy and prognosis of immunotherapy in thyroid cancer has now become a problem that needs to be solved urgently.

We analyzed the ICI of 505 THCA samples and divided the samples into three immune subtypes. Activated CD8⁺ T lymphocytes are critical for adaptive immune defense and kill cancer cells through a variety of mechanisms; they are generally considered the main cytotoxic lymphocytes exerting anti-tumor effects (Farhood et al., 2019; van der Leun et al., 2020). CD8⁺ T cell infiltration into the tumor indicates better prognosis (Bai et al., 2017; Ye et al., 2019; E et al., 2018). Tumor-associated dendritic cells (DCs) have functional defects, which can cause cancer immunosuppression (Veglia and Gabrilovich, 2017). In the TME, there are DCs with impaired antigen cross-presentation, and this can lead to changes in the activation and maintenance of anti-tumor immunity and promote tumor progression (Liu and Cao, 2015). In the present study, the ICI3 subtype with poor prognosis was characterized by low CD8⁺ T cell infiltration, high infiltration of static memory CD4⁺ T cells and active DCs, and a high stromal score. The ICI1 and ICI2 subtypes, which had good prognoses, were characterized by high infiltration of CD8⁺ T cells, low infiltration of regulatory T cells and activated DCs, and a low stromal score; this was consistent with previous results. Our survival analysis showed a higher OS of the high ICI score group, which is in agreement with previous studies on TME and patient prognosis (Li et al., 2019). We also found that most immune checkpoint-related and immune activation-related genes were upregulated in the high ICI score group, which indicates better immunogenicity and greater benefit from immunotherapy.

The gene set enrichment analysis (GSEA) was applied to study the biological distinctions between the two ICI score subgroups and to obtain the following immune-related pathways: chemokine signaling, cytokine receptor interactions, NK cell-mediated cytotoxicity, and T cell receptor signaling. Chemokines play important immunomodulatory roles by regulating cell proliferation, migration, activation, differentiation, and homing (Nagarsheth et al., 2017). NK cells mediate anti-tumor immunity through cytotoxicity and cytokine secretion (Müller-Hermelink et al., 2008; Gross et al., 2013). The TME of thyroid cancer produces soluble modulators that negatively regulate the maturation, proliferation, and effector functions of NK cells (Melaiu et al., 2020). It has also been shown that TCR signaling shapes the TCR composition of Foxp3+ regulatory T cells (Burchill et al., 2008). Research on cytokine receptor signaling pathways in thyroid cancer has thus far focused on the role of IL-1 family members (Mantovani et al., 2019).

Previous work has suggested that TMB is a prognostic marker and effective predictor of response towards immunotherapy in many cancers (Marabelle et al., 2020). The correlation between TMB and immune infiltration can be used to further assess the response to immunotherapy for a variety of tumors more effectively, including melanoma, bladder urothelial cancer, head and neck cancer, renal clear cell carcinoma, endometrial cancer, and others. Low TMB is associated with reduced immune infiltration, suggesting a poor prognosis for patients (Zhang et al., 2019; Zhang et al., 2020b; Jiang et al., 2020; Kang et al., 2020; Jiang et al., 2021; Zhou et al., 2021). In this study, ICI scores showed significant differences in survival either in high TMB or low TMB subgroups. Moreover, regardless of the expression in TMB, the prognosis of patients in the high ICI score group is always better than that of the low ICI score group, which shows that the TMB status will not interfere with the prediction based on the ICI score. These findings indicate that ICI score and TMB represent different aspects of tumor immunobiology, and ICI score may be a potential biomarker independent of TMB and effective in predicting response to immunotherapy.

There is prior evidence for gene mutations associated with response or tolerance to immunotherapy. For example, POLE/POLD1 gene mutations are associated with microsatellite instability (MSI), mismatch repair (MMR), and TMB, and CMTM6 is related to PD-L1 expression and regulation of anti-tumor immunity (Burr et al., 2017; Yao et al., 2019). In the present study, there were significant differences in gene mutation profiles between the different ICI score groups; this may provide new insight for studying the ICI composition of thyroid tumors and its relationship to gene mutations.

To further evaluate the value of the ICI score for estimating patient responses to immunotherapy, the IMvigor210, and GSE78220 cohorts of patients receiving immunotherapy were analyzed. As seen from the results, patients of the high ICI score group possessed higher objective responses to anti-PD-L1 therapy and greater overall survival. These results suggested that the ICI score relates to immunotherapy response and is an effective predictor of immunotherapy efficacy. Patients with high ICI scores may benefit more from immunotherapy.

Considering the enriched chemokine signaling pathways in the high ICI score group, the patients may also benefit from chemokine suppression and immune checkpoint blockade therapies. Previous studies have shown that the CXCR1/2 inhibitor Reparixin affects the viability of various thyroid cancer cells (Liotti et al., 2017). A CXCR4 antagonist (AMD3100) has also been shown to have significant anti-tumor effects on BHP10-3 papillary thyroid carcinoma cells in vivo and in vitro (Jung et al., 2016). Chemokine antagonists are expected to be adjuvant in the thyroid cancer’s treatment.

In the present study, we adapted retrospective research and used public databases and bioinformatic approaches to discover potential biomarkers related to the development of THCA and the efficacy of immunotherapy. Due to tumor heterogeneity, future studies should include more clinicopathological characteristics to improve the accuracy of predictions. Further validation of our approach is still needed in a potential large-scale cohort of THCA patients treated immunotherapy.

In summary, we revealed the characteristics of ICI in THCA and established an ICI score which could accurately predict the patient’s prognosis and response to immunotherapy. The results exhibited a comprehensive outlook of the internal immune landscape of THCA tumors.

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.

Ethics Statement

Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author Contributions

BJ designed the article, made important revisions to the paper and provided financial support; JG was responsible for data analysis and wrote the manuscript; LS was responsible for data collection; NL provided financial support and approved the final version of the paper to be published.

Funding

This research was supported by National Natural Science Foundation of China (Grant No. 81301767), Science and Technology Plan Project of Liaoning Province (Grant No. 20180551230).

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.714053/full#supplementary-material

References

Bai, M., Zheng, Y., Liu, H., Su, B., Zhan, Y., and He, H. (2017). CXCR5+ CD8+ T Cells Potently Infiltrate Pancreatic Tumors and Present High Functionality. Exp. Cel Res. 361 (1), 39–45. doi:10.1016/j.yexcr.2017.09.039

CrossRef Full Text | Google Scholar

Burchill, M. A., Yang, J., Vang, K. B., Moon, J. J., Chu, H. H., Lio, C.-W. J., et al. (2008). Linked T Cell Receptor and Cytokine Signaling Govern the Development of the Regulatory T Cell Repertoire. Immunity 28 (1), 112–121. doi:10.1016/j.immuni.2007.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Burr, M. L., Sparbier, C. E., Chan, Y.-C., Williamson, J. C., Woods, K., Beavis, P. A., et al. (2017). CMTM6 Maintains the Expression of PD-L1 and Regulates Anti-tumour Immunity. Nature 549 (7670), 101–105. doi:10.1038/nature23643

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

Capdevila, J., Wirth, L. J., Ernst, T., Ponce Aix, S., Lin, C.-C., Ramlau, R., et al. (2020). PD-1 Blockade in Anaplastic Thyroid Carcinoma. Jco 38 (23), 2620–2627. doi:10.1200/JCO.19.02727

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2016). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 18 (1), 248–262. doi:10.1101/056101

PubMed Abstract | CrossRef Full Text | Google Scholar

Coperchini, F., Croce, L., Marinò, M., Chiovato, L., and Rotondi, M. (2019). Role of Chemokine Receptors in Thyroid Cancer and Immunotherapy. Endocr. Relat. Cancer 26 (8), R465–R478. doi:10.1530/ERC-19-0163

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, L., Wang, Y., Sun, X., Li, H., Geng, X., Ge, M., et al. (2018). Thyroid Cancer: Trends in Incidence, Mortality and Clinical-Pathological Patterns in Zhejiang Province, Southeast China. BMC Cancer 18 (1), 291. doi:10.1186/s12885-018-4081-7

PubMed Abstract | CrossRef Full Text | Google Scholar

E, J., Kang, Z., Zhu, L., Xing, J., and Yu, E. (2018). CD8+CXCR5+ T Cells in Tumor-Draining Lymph Nodes Are Highly Activated and Predict Better Prognosis in Colorectal Cancer. Hum. Immunol. 79 (6), 446–452. doi:10.1016/j.humimm.2018.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, W., Ye, L., Shen, L., Cai, J., Huang, F., Wei, Q., et al. (2014). Tumor-associated Macrophages Promote the Metastatic Potential of Thyroid Papillary Cancer by Releasing CXCL8. Carcinogenesis 35 (8), 1780–1787. doi:10.1093/carcin/bgu060

PubMed Abstract | CrossRef Full Text | Google Scholar

Farhood, B., Najafi, M., and Mortezaee, K. (2019). CD8 + Cytotoxic T Lymphocytes in Cancer Immunotherapy: A Review. J. Cel Physiol 234 (6), 8509–8521. doi:10.1002/jcp.27782

CrossRef Full Text | Google Scholar

Ferrari, S. M., Fallahi, P., Galdiero, M. R., Ruffilli, I, Elia, G., Ragusa, F., et al. (2019). Immune and Inflammatory Cells in Thyroid Cancer Microenvironment. Ijms 20 (18), 4413. doi:10.3390/ijms20184413

PubMed Abstract | CrossRef Full Text | Google Scholar

Gogali, F., Paterakis, G., Rassidakis, G. Z., Kaltsas, G., Liakou, C. I., Gousis, P., et al. (2012). Phenotypical Analysis of Lymphocytes with Suppressive and Regulatory Properties (Tregs) and NK Cells in the Papillary Carcinoma of Thyroid. J. Clin. Endocrinol. Metab. 97 (5), 1474–1482. doi:10.1210/jc.2011-1838

CrossRef Full Text | Google Scholar

Gross, E., Sunwoo, J. B., and Bui, J. D. (2013). Cancer Immunosurveillance and Immunoediting by Natural Killer Cells. Cancer J. 19 (6), 483–489. doi:10.1097/PPO.0000000000000005

PubMed Abstract | CrossRef Full Text | Google Scholar

Iyer, P. C., Dadu, R., Gule-Monroe, M., Busaidy, N. L., Ferrarotto, R., Habra, M. A., et al. (2018). Salvage Pembrolizumab Added to Kinase Inhibitor Therapy for the Treatment of Anaplastic Thyroid Carcinoma. J. Immunotherapy Cancer 6 (1), 68. doi:10.1186/s40425-018-0378-y

CrossRef Full Text | Google Scholar

Jiang, A.-M., Ren, M.-D., Liu, N., Gao, H., Wang, J.-J., Zheng, X.-Q., et al. (2021). Tumor Mutation burden, Immune Cell Infiltration, and Construction of Immune-Related Genes Prognostic Model in Head and Neck Cancer. Int. J. Med. Sci. 18 (1), 226–238. doi:10.7150/ijms.51064

CrossRef Full Text | Google Scholar

Jiang, F., Wu, C., Wang, M., Wei, K., Zhou, G., and Wang, J. (2020). Multi-omics Analysis of Tumor Mutation burden Combined with Immune Infiltrates in Melanoma. Clinica Chim. Acta. 511, 306–318. doi:10.1016/j.cca.2020.10.030

CrossRef Full Text | Google Scholar

Jung, Y. H., Lee, D. Y., Cha, W., Kim, B. H., Sung, M.-W., Kim, K. H., et al. (2016). Antitumor Effect of CXCR4 Antagonist AMD3100 on the Tumorigenic Cell Line of BHP10-3 Papillary Thyroid Cancer Cells. Head Neck 38 (10), 1479–1486. doi:10.1002/hed.24461

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

Kim, S., Cho, S. W., Min, H. S., Kim, K. M., Yeom, G. J., Kim, E. Y., et al. (2013). The Expression of Tumor-Associated Macrophages in Papillary Thyroid Carcinoma. Endocrinol. Metab. 28 (3), 192. doi:10.3803/enm.2013.28.3.192

PubMed Abstract | CrossRef Full Text | Google Scholar

La Vecchia, C., Malvezzi, M., Bosetti, C., Garavello, W., Bertuccio, P., Levi, F., et al. (2015). Thyroid Cancer Mortality and Incidence: A Global Overview. Int. J. Cancer 136 (9), 2187–2195. doi:10.1002/ijc.29251

CrossRef Full Text | Google Scholar

Li, B., Cui, Y., Nambiar, D. K., Sunwoo, J. B., and Li, R. (2019). The Immune Subtypes and Landscape of Squamous Cell Carcinoma. Clin. Cancer Res. 25 (12), 3528–3537. doi:10.1158/1078-0432.CCR-18-4085

PubMed Abstract | CrossRef Full Text | Google Scholar

Liotti, F., de Pizzol, M., Allegretti, M., Prevete, N., and Melillo, R. M. (2017). Multiple Anti-tumor Effects of Reparixin on Thyroid Cancer. Oncotarget 8 (22), 35946–35961. doi:10.18632/oncotarget.16412

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., and Cao, X. (2015). Intratumoral Dendritic Cells in the Anti-tumor Immune Response. Cell Mol Immunol. 12 (4), 387–390. doi:10.1038/cmi.2014.130

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, M., Lin, B., Wang, M., Liang, X., Su, L., Okose, O., et al. (2020). Immunotherapy in Anaplastic Thyroid Cancer. Am. J. Transl Res. 12 (3), 974–988.

Google Scholar

Mantovani, A., Dinarello, C. A., Molgora, M., and Garlanda, C. (2019). Interleukin-1 and Related Cytokines in the Regulation of Inflammation and Immunity. Immunity 50 (4), 778–795. doi:10.1016/j.immuni.2019.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Marabelle, A., Fakih, M., Lopez, J., Shah, M., Shapira-Frommer, R., Nakagawa, K., et al. (2020). Association of Tumour Mutational burden with Outcomes in Patients with Advanced Solid Tumours Treated with Pembrolizumab: Prospective Biomarker Analysis of the Multicohort, Open-Label, Phase 2 KEYNOTE-158 Study. Lancet Oncol. 21 (10), 1353–1365. doi:10.1016/S1470-2045(20)30445-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Melaiu, O., Lucarini, V., Cifaldi, L., and Fruci, D. (2020). Influence of the Tumor Microenvironment on NK Cell Function in Solid Tumors. Front. Immunol. 10, 3038. doi:10.3389/fimmu.2019.03038

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller-Hermelink, N., Braumüller, H., Pichler, B., Wieder, T., Mailhammer, R., Schaak, K., et al. (2008). TNFR1 Signaling and IFN-γ Signaling Determine whether T Cells Induce Tumor Dormancy or Promote Multistage Carcinogenesis. Cancer Cell. 13 (6), 507–518. doi:10.1016/j.ccr.2008.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Na, K. J., and Choi, H. (2018). Immune Landscape of Papillary Thyroid Cancer and Immunotherapeutic Implications. Endocr. Relat. Cancer 25 (5), 523–531. doi:10.1530/ERC-17-0532

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagarsheth, N., Wicha, M. S., and Zou, W. (2017). Chemokines in the Cancer Microenvironment and Their Relevance in Cancer Immunotherapy. Nat. Rev. Immunol. 17 (9), 559–572. doi:10.1038/nri.2017.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer Statistics, 2020. CA A. Cancer J. Clin. 70 (1), 7–30. doi:10.3322/caac.21590

CrossRef Full Text | Google Scholar

Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade to Improve Prognosis. J. Natl. Cancer Inst. 98 (4), 262–272. doi:10.1093/jnci/djj052

CrossRef Full Text | Google Scholar

van der Leun, A. M., Thommen, D. S., and Schumacher, T. N. (2020). CD8+ T Cell States in Human Cancer: Insights from Single-Cell Analysis. Nat. Rev. Cancer. 20 (4), 218–232. doi:10.1038/s41568-019-0235-4

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

Veglia, F., and Gabrilovich, D. I. (2017). Dendritic Cells in Cancer: the Role Revisited. Curr. Opin. Immunol. 45, 43–51. doi:10.1016/j.coi.2017.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Z., Li, X., He, Y., Wu, S., Wang, S., Sun, J., et al. (2020). Immune Cell Confrontation in the Papillary Thyroid Carcinoma Microenvironment. Front. Endocrinol. 11, 570604. doi:10.3389/fendo.2020.570604

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, J., Gong, Y., Zhao, W., Han, Z., Guo, S., Liu, H., et al. (2019). Comprehensive Analysis of POLE and POLD1 Gene Variations Identifies Cancer Patients Potentially Benefit from Immunotherapy in Chinese Population. Sci. Rep. 9 (1), 15767. doi:10.1038/s41598-019-52414-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, L., Li, Y., Tang, H., Liu, W., Chen, Y., Dai, T., et al. (2019). CD8+CXCR5+T Cells Infiltrating Hepatocellular Carcinomas Are Activated and Predictive of a Better Prognosis. Aging 11 (20), 8879–8891. doi:10.18632/aging.102308

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, H., Tang, Y., Guo, Y., and Wen, S. (2020). Immune Microenvironment of Thyroid Cancer. J. Cancer 11 (16), 4884–4896. doi:10.7150/jca.44506

CrossRef Full Text | Google Scholar

Zhang, C., Li, Z., Qi, F., Hu, X., and Luo, J. (2019). 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Shen, L., Qi, F., Wang, J., and Luo, J. (2020). Multi‐omics Analysis of Tumor Mutation burden Combined with Immune Infiltrates in Bladder Urothelial Carcinoma. J. Cel Physiol. 235 (4), 3849–3863. doi:10.1002/jcp.29279

CrossRef Full Text | Google Scholar

Zhang, X., Shi, M., Chen, T., and Zhang, B. (2020). Characterization of the Immune Cell Infiltration Landscape in Head and Neck Squamous Cell Carcinoma to Aid Immunotherapy. Mol. Ther. - Nucleic Acids. 22, 298–309. doi:10.1016/j.omtn.2020.08.030

PubMed Abstract | 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