CXCL10 is a potential biomarker and associated with immune infiltration in human papillary thyroid cancer

Thyroid cancer is not only the most common malignancy of the human endocrine system, accounting for ∼3% of the total incidence of systemic malignant tumors, but also one of the fastest growing malignant solid tumors. According to the histopathological classification system, thyroid cancer is mainly categorized into papillary thyroid, follicular thyroid, medullary thyroid and undifferentiated thyroid cancer types. Among these, papillary thyroid carcinoma accounts for approximately 80% of all pathological thyroid cancer cases [1]. In 2018, approximately 52,070 new cases of thyroid cancer in the United States were recorded [2]. Due to improvements in diagnostic techniques, such as ultrasound and fine-needle biopsy, the incidence of PTC has been increasingly reported in recent years [3]. The majority of PTC patients have slow disease progression and good prognosis, with a 10-year survival rate of more than 90%. However, patients with papillary thyroid microcarcinoma commonly develop lymph node metastasis at an early stage of disease with an incidence of 30–40%, leading to a 50% decrease in the 10-year survival rate [4–6]. Lymph node metastasis from thyroid cancer cells serves as an independent risk factor for poor prognosis of PTC, with more than 8–10% eventually evolving to distant metastasis or potential cancer-related death [7–10]. Early and effective interventions are critical for diagnosis of PTC and blocking the development of lymph node metastasis in thyroid cancer cells, which would significantly improve patient outcomes. Therefore, identification of biomarkers related to the tumor microenvironment is currently an important consideration in the clinical management of PTC.

Recent studies suggest that immune cells, such as macrophages, mast cells, neutrophils, and lymphocytes, play a tumorigenic role in PTC, supporting the potential efficacy of targeted immunotherapy as a treatment strategy [11]. Immune escape, an important feature of tumor malignancy, plays a critical regulatory role in the evolution of many solid tumor types. The term ‘tumor microenvironment’ is used to describe the milieu surrounding tumor cells composed of inflammatory mediators, extracellular stromal molecules, and immune cells [12]. Essential non-tumor components of the tumor microenvironment are immune and stromal cells, which are critical for transmitting information and stimulating the formation and evolution of tumors. The diagnostic and prognostic value of immune and stromal cells in tumors has been established [13,14]. Tumor immune escape can be modulated by immune cell components in the microenvironment, where by alterations in these components promote the formation of an immunosuppressive state [15]. Clinicopathological typing of thyroid cancer patients is reported to be closely related to immune and inflammatory cells in the tumor microenvironment [16,17]. The ESTIMATE algorithm has been further used to evaluate immune and stromal scores through analysis of expression patterns of specific genes in immune and stromal cells for the purpose of predicting cancer cell infiltration into non-tumor cells [18]. The ESTIMATE algorithm has been increasingly applied to a variety of tumor models, including breast and colon cancer types [19,20]. However, to our knowledge, a few studies have explored the ratio of stromal and immune cells in PTC using the ESTIMATE algorithm.

In the present study, we evaluated immune scores by extracting tumor microenvironment-related genes based on TCGA data of PTC patients and the ESTIMATE algorithm and explored one or several gene modules playing crucial roles in the PTC microenvironment, with a view to providing a theoretical basis for further research.

Gene expression data for thyroid cancer patients were retrieved from the TCGA database (https://tcga-data.nci.nih.gov/tcga/). The RNA expression profiles of PTC patients were further analyzed and annotation data of the Affymetrix HT human genome U133 array plate set extracted. Clinical information, including gender, age, histologic grade, pathologic stage, survival and outcome, were additionally extracted from TCGA. After downloading the data, the ESTIMATE algorithm was implemented in the form of the ‘ESTIMATE’ package in R for calculation of the immune and stromal scores [18].

To improve screening for differentially expressed genes, the limma software package was employed [21]. Differentially expressed genes (DEGs) were identified based on cutoffs of fold change > 2 and adj. P<0.05.

To construct cluster analysis and immune stromal heatmaps, we used pheatmap R software package [22].

To analyze protein–protein interactions (PPI), the STRING database was employed [23]. To clarify the interactions between genes, protein–protein interaction networks of genes that were up- or down-regulated were constructed based on immune and stromal scores and reconstructed via Cytoscape software. The active interaction source settings were as follows: text mining, experiment, database, co-expression, neighborhood, gene fusion and co-occurrence. The minimum interaction score requirement was 0.4. Only a single network consisting of more than 10 nodes was shown and allowed to proceed to the next level of analysis. The MCODE plug in was subsequently used to analyze closely related gene modules.

To evaluate the prognostic value of DEGs in thyroid cancer, Kaplan–Meier survival analysis was used and overall survival curves obtained using Kaplan–Meier plotter.

To establish the potential biological functions of intersecting DEGs, KEGG pathway enrichment [24] and GO [25] analyses were performed using the Profiler package. GO assessed biological processes (BP), molecular functions (MF) and cellular components (CC). The results from both GO and KEGG pathway analyses were processed using ‘clusterProfiler’, ‘enrichplot’ and ‘ggplot2’ packages. The results were visualized with the ggplot2 package. The P value cut-off was set at 0.05.

The LinkedOmics (http://www.linkedomics.org/) database is a unique platform for accessing, analyzing and comparing multi-omics data within and between tumor types, providing a convenient resource for biologists and clinicians [26]. To explore the kinase target of the PTC hub gene, CXCL10, the LinkInterpreter module of LinkedOmics was applied in the present study.

The TIMER (https://cistrome.shinyapps.io/timer/) database is a practical platform for systematic evaluation of the interactions between tumor and immune cells [27]. In the current study, we aimed to explore the associations among immune infiltrates, immunomodulatory factors and CXCL10 gene expression. All analyses were based on the TCGA TCHA dataset (N= 501).

Open targets provide a target-centric workflow to aid in identifying diseases potentially related to specific targets [28]. Here, we used Open targets to explore diseases related to CXCL10.

GeneMANIA (www.genemania.org) is a portal to establish protein–protein interaction (PPI) networks and provide insights into the functions of submitted genes [29]. In the present study, GeneMANIA was used to visualize gene networks of CXCL10 and its related genes.

GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) is a portal website providing a platform for analysis of gene sets in cancer. Genomic aberrations can affect clinical response to treatment and serve as potential biomarkers for drug screening [30]. We integrated 265 small molecules from Cancer Drug Sensitivity Genomics (GDSC) and analyzed the expression of each gene in the genome in relation to small molecule/drug sensitivity (IC₅₀) via Spearman correlation analysis. P-values less than 0.05 were considered statistically significant.

GSEA is a calculation method that performs genome-wide expression profiling of two sample types, with the aim of identifying significant and consistent differences. In this study, according to CXCL10 expression patterns, samples of PTC were divided into high and low expression groups. GSEA was subsequently applied to evaluate the differences between the high and low CXCL10 expression groups. Data were considered significant at P<0.05 and the false discovery rate (FDR) was <25%. According to the normalized enrichment score (NES), an enrichment approach related to the biological process of papillary thyroid carcinoma was selected.

For IHC, PTC tissue microarrays were firstly deparaffinized with xylene and rehydrated with a gradient of ethanol to distilled water. Antigen retrieval was performed by heating the tissue microarrays in 10 mM sodium citrate buffer for 30 min. Next, endogenous peroxidase was blocked with 3% H₂O₂ and the microarray treated with normal goat serum working solution for 30 min to reduce non-specific binding. After overnight incubation with the rabbit monoclonal anti-CXCL10 antibody (1:100 dilution, BOSTER, Wu Han, China, BA4723) at 4°C, 120-point PTC tissue microarrays were washed with PBS and incubated with broad-spectrum secondary antibodies. Slides were further washed with PBS, followed by incubation with horseradish enzyme-labeled streptavidin working solution at 37°C for 30 minutes. After a further three washes with PBS, signal detection was performed with the DAB staining kit. According to Remmele's semi-quantitative immune response score (IRS) scale, the overall IHC score (1–5) was evaluated [31].

In March 2020, the gene expression profiles and corresponding clinicopathological data from 518 PTC patients, including 369 (70.89%) female patients and 149 (29.11%) male patients, were extracted from the TCGA database. The minimum age of initial diagnosis was 15 years and average age was 46 years. PTC stages were distributed as follows: I (55.98%, n=290), II (11%, n=57), III (23.16%, n=120) and IV (9.84%, n=51). Using the ESTIMATE method, immune scores were in the range of -1304.72 to 3233.39, stromal scores were -1714.40 to 1608.17and ESTIMATE scores were -2470.47 to 4199.50. With increase in the tumor staging level, the corresponding average immune (P=0.044), stromal and ESTIMATE (P=0.046) scores were increased (Figure 1). Next, Kaplan–Meier survival analysis was conducted to establish the associations of immune, stromal and ESTIMATE scores with prognosis. Patients with PTC were divided into high and low score groups. However, none of the scores were associated with survival outcomes in PTC (Figure 1).

To obtain the gene expression profiles of patients with PTC in association with immune and stromal scores, Affymetrix microarray data analysis was conducted on 518 PTC samples from patients in the TCGA database. The population with a higher immune score, presented in heatmap2A, had a total of 1269 DEGs, including 954 up-regulated and 315 down-regulated genes. Based on the higher stromal score (Figure 2B), 986 DEGs were identified, including 914 up-regulated and 72 down-regulated genes. Venn diagram analysis disclosed that compared with the low group, the high group had a total of 854 DEGs, including 789 up-regulated and 65 down-regulated genes (Figure 2C,D).

To establish the functions of DEGs, enrichment analysis using GO and KEGG functions was conducted. GO analysis demonstrated that the main ‘Biological Process’ (BP) ontology was enriched in leukocyte migration, leukocyte chemotaxis, cell chemotaxis, chemokine-mediated signaling pathway, response to chemokine, cellular response to chemokine, myeloid leukocyte migration, neutrophil chemotaxis, neutrophil migration and granulocyte chemotaxis functions. For the ‘Cellular Component’ (CC) ontology, enriched terms included the external side of the plasma membrane, plasma membrane receptor complex, receptor complex, membrane raft, membrane microdomain, membrane region, T-cell receptor complex, immunological synapse, protein complex involved in cell adhesion and plasma membrane raft. With regard to ‘Molecular Function’ (MF), cytokine receptor binding, chemokine receptor binding, cytokine activity, G protein-coupled receptor binding, receptor ligand activity, chemokine activity, CCR chemokine receptor binding, CXCR chemokine receptor binding, C-C chemokine receptor activity and C-G chemokine binding terms were enriched (Figure 3A). KEGG pathway analysis additionally disclosed involvement of DEGs in cytokine–cytokine receptor interactions and the chemokine signaling pathway (Figure 3B). PPI network analysis was an important analytical tool to identify the critical genes associated with PTC interactions from a systems perspective. In the present study, the top three PPI networks of CXCL10, KRT, and ZAP70 modules were selected using Cytoscape software (Figure 3C).Greater connectivity was signified by larger gene nodes. The top five connected genes were CXCL10, CCL20, CCL5, CD4, and CXCL1 (Figure 3D).

To establish the relationships between the 789 up-regulated and 65 down-regulated genes and prognostic survival of PTC patients, Kaplan–Meier survival analysis was conducted. All TC tumor samples were divided into high and low expression groups based on specific genes. Overall, 15 genes were associated with prognosis and survival. Groups with high expression of GATA5 (HR = 2.95, 95% CI: 1.11–7.89, P=0.023), LRRN4CL (HR = 10.13, 95% CI: 1.34–76.75, P=0.0055), OGDHL (HR = 4.2, 95% CI: 1.55–11.93, P=0.0022), PSAT1 (HR = 3.57, 95% CI: 1.32–9.65, P=0.0073), SLC25A47 (HR = 3.45, 95% CI: 1.28–9.32, P=0.0093) and TWIST2 (HR = 3.23, 95% CI: 1.04–10.03 and P=0.032) showed poorer overall survival rates while higher survival rates were associated with high expression of CXCL10 (HR = 0.24, 95% CI: 0.07–0.86, P=0.017), ADGRG5 (HR = 0.17, 95% CI: 0.05–0.59, P=0.0016), ASB2 (HR = 0.37, 95% CI: 0.14–0.97, P=0.036), DMBT1 (HR = 0.15, 95% CI: 0.03–0.66, P=0.0038), GPR34 (HR = 0.17, 95% CI: 0.05–0.61, P=0.0019), GZMK (HR = 0.3, 95% CI: 0.11–0.8, P=0.011), GZMM (HR = 0.2, 95% CI: 0.07–0.58, P=0.00094), HTRA1 (HR = 0.25, 95% CI: 0.09–0.72, P=0.0056) and TBX21 (HR = 0.28, 95% CI: 0.09–0.88, and P=0.02; Figure 4).Accordingly, we speculate that expression levels of these 15 genes are critical for prognosis of PTC patients. Based on the highest degree, CXCL10 was selected as the core gene in the PTC immune microenvironment for further exploration.

The genes related to CXCL10 and differentially expressed in PTC were investigated via LinkedOmics to establish the specific mechanisms of action of CXCL10 in thyroid cancer. As a result, 19,928 genes related to CXCL10 were identified. Among these, 10,157 genes (dark red dots) were positively correlated and 9769 genes (dark green dots) were negatively correlated with CXCL10 (Figure 5A; false discovery rate <0.05). The top 50 genes that were significantly positively and negatively correlated with CXCL10 are presented in Figure 5B,C, respectively. As shown in Figure 5D, CXCL10 and neighboring genes were mainly implicated in adaptive immune response, cellular defense response, T-cell activation, response to interferon-gamma, immune response, regulation of signaling pathways, antigen processing and presentation, and regulation of the innate immune response. The results of KEGG pathway analysis showed that the functions of CXCL10 and its neighboring genes were mainly enriched in the NF-kappaB signaling pathway, allograft rejection, primary immunodeficiency and leishmaniasis (Figure 5E). To further establish the diseases caused by aberrant expression of the CXCL10 gene, an analysis using the Open Targets website was conducted. As expected, CXCL10 was significantly implicated in diseases of the immune and endocrine systems (Figure 6). Our findings demonstrate that the CXCL10 gene plays an important regulatory role in immune pathways of the PTC microenvironment, supporting its utility as a novel immunotherapeutic target for PTC. The potential mechanisms of action of CXCL10 were further explored.

To determine the CXCL10-associated signaling pathways involved in PTC, we conducted gene set enrichment analysis (GSEA) to determine differences in CXCL10 gene expression. The most significantly enriched signal pathway was identified based on the normalized enrichment scores (NES). High expression of CXCL10 was markedly associated with activation of oxidative phosphorylation, glycolysis gluconeogenesis, glycerolipid metabolism, PPAR signaling, fatty acid metabolism and calcium signaling pathways and inhibition of p53 signaling, cell cycle, pathways in cancer, and colorectal cancer signal pathways (Figure 7).

The potential kinase targets of CXCL10 in thyroid carcinoma were analyzed using the LinkedOmics database. The kinase network corresponding to CXCL10 included LCK, LYN, HCK and SYK (Table 1). To clarify the potential regulatory mechanisms between these four kinases and CXCL10, PPIs were constructed using GeneMANIA to establish their specific functions. The gene sets were mainly implicated in regulation of lymphocyte activation, regulation of leukocyte activation, regulation of T-cell activation, regulation of cell activation, antigen receptor-mediated signaling and T-cell receptor signaling pathways (Figure 8). Next, we performed tissue microarray immunohistochemical staining to evaluate expression of CXCL10 in PTC and matched adjacent tissues. Compared with adjacent tissues, expression of CXCL10 in PTC tissues was clearly decreased (Figure 9).

We further analyzed the correlation between CXCL10 and immune cells infiltrating the tumor microenvironment using the TIMER database. CXCL10 expression was positively associated with B cell infiltration (Cor = 0.394, P = 2.16e-19), CD8 + T cells (Cor = 0.535, P = 1.98e-37), CD4 + T cells (Cor = 0.177, P = 8.55-05), macrophages (Cor = 0.338, P = 1.53e-14), neutrophils (Cor = 0.624, P = 5.84e-54) and dendritic cells (Cor = 0.749, P = 1.90e-88; Figure 10). The relationship between CXCL10 expression and a number of immune markers was further evaluated. As shown in Table 2, significant correlations (correlation coefficient> 0.6) were observed between CXCL10 and expression of PD-1 (PDCD1), CD2, CD19, CD79A, CD68, IL10, VSIG4, HLA-DPB1, CD11c (ITGAX), T-bet (TBX21), STAT1, IFN-g (IFNG), FOXP3 and CTLA4.

For drug sensitivity analysis, we assessed the correlation between CXCL10 expression and IC₅₀ values of molecules from the Cancer Drug Sensitivity Genomics (GDSC) database. Our data showed that four drugs or small molecules were effective in association with reduced expression of CXCL10 (Figure 11). Specifically, CXCL10 was negatively regulated by docetaxel, BRD-K30748066, BRD-K01737880, and BRD-A86708339, supporting its utility as a potential therapeutic drug target for PTC.

PTC is a common malignant tumor of the endocrine system with increasing annual incidence. The tumor immune microenvironment is known to play a critical role in the evolution of malignancy and related to various biological behaviors of cancer, affecting tumor growth and spread. In the present study, we extracted data from the TCGA database to explore genes in the tumor microenvironment affecting overall survival, tumorigenesis and development of PTC. Furthermore, biomarkers related to prognosis and treatment in the microenvironment of PTC were identified. The relationships between clinicopathological variables and immune and stromal scores in PTC were initially evaluated. As shown in Figure 1, immune and stromal scores were related to tumor stage. With increase in the tumor staging level, the corresponding average immune and stromal scores as well as ESTIMATE score were increased, supporting the potential value of immune and stromal scores in the treatment and prediction of PTC. After grouping of all samples based on high and low immunity/stromal scores, 854 differentially expressed genes (DEGs) were identified. To clarify the roles of these DEGs in immune and stromal systems, we explored their functions in PTC. GO analysis disclosed the involvement of DEGs in leukocyte migration, leukocyte chemotaxis, cell chemotaxis, chemokine-mediated signaling pathways, response to chemokines, cellular response to chemokines and myeloid leukocyte migration. KEGG pathway analysis showed that DEGs are involved in cytokine–cytokine receptor interactions and the chemokine signaling pathway. Both GO and KEGG data suggest that the functions of these DEGs are closely related to specific pathways and molecular functions in the tumor microenvironment in PTC. Our results are in agreement with previous studies showing that aberrant expression of specific genes contributes significantly to development of the PTC tumor microenvironment through regulatory effects on the activities of various immune cells together with extracellular matrix molecules [32,33].

To further clarify the interactions between DEGs in PTC, we used STRING and Cytoscape software to construct protein–protein interaction networks and identify central genes. As shown in the results, the top five connected degree genes were CXCL10, CCL20, CCL5, CD4 and CXCL1. The top three PPI networks of CXCL10, KRTs, and ZAP70 modules were further selected as core modules. Immune evasion is a signature of tumorigenesis [34]. An earlier study by Marina et al. [35] reported that the gene encoding the chemokine CCL20 is overexpressed in PTC compared with normal thyroid tissue and independent of oncogene (RET/PTC or BRAF) status, suggesting the involvement of this chemokine in tumor-associated inflammation and the microenvironment. Several researchers have additionally demonstrated that four chemokines, CXCL1, CCL5, CCL20 and CCL21, biologically affect thyroid cancer cells by modulating inflammation and the immune microenvironment [36–40]. Abnormal secretion of CXCL10 in PTC and its involvement in regulation of malignant biological behavior have been previously reported [41,42]. Here, we further explored the association between DEGs and prognosis of PTC and identified 15 genes related to the tumor microenvironment. High expression of GATA5, LRRN4CL, OGDHL, PSAT1, SLC25A47, and TWIST2 was associated with poorer overall survival rates and high expression of CXCL10, ADGRG5, ASB2, DMBT1, GPR34, GZMK, GZMM, HTRA1 and TBX21with longer overall survival rates in thyroid carcinoma patients. Accumulating evidence suggests that aberrant expression of CXCL family peptides leads to abnormal recruitment of immune cells in tumors and various types of immune cells can increase CXCL peptide secretion by modulating the corresponding signaling pathways (particularly NF-kappaB), which suppress the immune surveillance status of the organism, ultimately resulting in a vicious cycle and reduced survival of cancer patients [43,44]. In view of the finding that the degree of connectivity of CXCL10 gene was up to 23 and associated with PTC survival, we selected CXCL10 as a core gene for analysis. Our data suggest that the CXCL10 gene plays an important regulatory role in the immune microenvironment of PTC, supporting its utility as a novel immunotherapeutic target. Next, we focused on the mechanisms underlying the function of CXCL10 in PTC. Immunochemical staining experiments revealed that CXCL10 was expressed at lower levels in PTC than paracancerous tissues. GSEA pathway enrichment analysis further showed low expression of CXCL10 in multiple cancer pathways (p53 signaling pathway, Cell cycle, Pathways in cancer, Colorectal cancer signal pathway).The LinkedOmics database was used to explore the top 50 positively and negatively correlated genes and kinase targets of CXCL10.As shown in Figure 5B,C, the top five genes negatively associated with CXCL10 were CXCL11, CXCL9, GBP5, UBD, and IL12RB1 and the top five genes positively associated with CXCL10 were CMTM4, TOM1L2, LCLAT1, CRY2, and PCBD2 with CXCL10. LCK, LYN, HCK, and SYK were the main kinase targets of CXCL10. All four kinases, along with CXCL10 were significantly involved in T-cell activation, regulation of lymphocyte activation, regulation of leukocyte activation, regulation of T-cell activation, regulation of cell activation, antigen receptor-mediated signaling and T-cell receptor signaling pathways. The kinase LCK belongs to the Src family and is critical for T-cell proliferation and activation in physiological conditions. Under pathological conditions, aberrantly expressed LCKs can directly trigger an inflammatory immune response [45]. In addition, LCK has been shown to be overexpressed in many cancer types and LCK inhibitors have been clinically used to treat a number of solid cancers [46–48]. Another study by Silvia and co-workers demonstrated that LCK kinases are involved in regulating the progression of thyroid cancer [49]. The collective results suggest that kinases play regulatory roles in the evolution of cancer. In many solid tumors, the degree of immune cell infiltration is strongly associated with patient prognosis [50,51]. In our experiments, Significant correlations were detected between CXCL10 and expression of immune markers, including PD-1 (PDCD1), CD2, CD19, CD79A, CD68, IL10, VSIG4, HLA-DPB1, CD11c (ITGAX), T-bet (TBX21), STAT1, IFN-g (IFNG), FOXP3 and CTLA4, suggesting that the CXCL10 gene exerts its effects in the thyroid cancer microenvironment by influencing immune cells and related molecules. Our results showed that CXCL10 is related to the overall survival of PTC patients and positively correlated with the infiltration level of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils and dendritic cells. While the CXCL10 gene could be used to provide additional information on the correlation between immune cell infiltration and clinical outcomes of PTC patients, the specific mechanisms of action of CXCL10 in PTC remain to be elucidated. The tumor microenvironment is composed of immune and stromal cells that participate in the progression of cancer and interact with PTC, thereby affecting occurrence, development, migration, metastasis and prognosis of cancer. The present study focused on the correlation between immune/stromal scores and PTC clinicopathological variables and explored the functions and prognostic value of genes related to the tumor microenvironment. We selected CXCL10 as a potential core gene of the PTC microenvironment and showed its significant correlation with prognosis and immune cells in PTC. CXCL10 expression was additionally regulated by a number of small-molecule drugs, including docetaxel, BRD-K30748066, BRD-K01737880, and BRD-A86708339. The collective evidence supports the utility of CXCL10 as a potential therapeutic drug target and novel potential biomarker for improving management and prognosis of PTC.

The data used to support the findings of this study are available from the corresponding author upon request.

The authors declare that there are no competing interests associated with the manuscript.

This study was funded by the Hebei North University 2019 Basic Research Business Expenses Project [grant number JYT2019015]; Hebei Provincial Department of Finance Specialist Capacity Building and Specialist Leadership Program, China [grant number 361009]; and Zhangjiakou City's 2019 Key R & D Plan Projects [grant number 1911019D].

All authors approved the final draft of the manuscript. Xu Lin designed the study. Qin Xiao-Jing, Xu Lin and Gang Xue contributed equally to the writing of the manuscript and should be considered co-first authors. Hui-Li Fan and Hao-Yu Wang directed the experiments.

DEG

differentially expressed gene

IHC

immunohistochemical

PPI

protein–protein interaction

PTC

papillary thyroid cancer

TC

thyroid cancer