Identification of Key Genes in Thyroid Cancer Microenvironment.

BACKGROUND Tumor microenvironment (TME) plays important roles in the development of cancer. However, the roles of TME in thyroid cancer are not well studied. In our study, we aimed to identify genes related to thyroid cancer microenvironment. MATERIAL AND METHODS We combined The Cancer Genome Atlas (TCGA) and Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) datasets to identify differentially expressed genes in thyroid cancer microenvironment. Then, using these differentially expressed genes, we constructed protein-protein interaction (PPI) network and conducted functional enrichment analysis. Genes with degree beyond 12 in the PPI network were regarded as hub genes. Finally, we conducted Kaplan-Meier curve and log-rank test and functional enrichment analysis on these hub genes. RESULTS There were 793 differentially expressed genes identified to be associated with immune score and stromal score in thyroid cancer microenvironment. We screened out 30 hub genes by construction of PPI network. The functions of these hub genes were enriched in immune cell activity, cytokine and chemokine activity, cell adhesion molecules, and extracellular matrix, which provided further insight into the roles of these genes in the tumor microenvironment. CXCL10, with the highest degrees in the PPI network, were positively related to overall survival of thyroid cancer patients (P=0.02467). CONCLUSIONS We identified 30 tumor microenvironment related genes in thyroid cancer. Among these hub genes, CXCL10 can be regarded as a prognostic biomarker in thyroid cancer.

Background

Thyroid cancer is common all over the world and the incidence continues to rise. According to the statistics of International Agency for Research on Cancer (IARC), 567 233 new cases of thyroid cancer occurred worldwide in 2018. About 76.9% of new cases were women, accounting for 5.1% of all cancer types in women [1]. Thyroid cancer owes a broad range of prognosis. Commonly, thyroid cancer has a good outcome following standard treatments; however, thyroid cancer can also exhibit very aggressive behavior with high mortality [2]. So, it is of great importance to find out how novel biomarkers contribute to clinical decision-making and prognosis prediction in thyroid cancer.

Tumor microenvironment (TME), one of the hallmarks of cancer, is defined as the cellular and physical environment surrounding the primary tumor. There are many kinds of cells and molecules in TME including inflammatory and immune cells, extracellular matrix proteins, soluble factors and so on [3]. Tumor cells act in close interaction with molecules and cells in the TME, which profoundly influences the biological behavior of cancer cells and the prognosis of cancer patients [4,5]. A number of therapeutic strategies that target the TME are being test in clinical trials and showed promising results [6,7]. Therefore, we aim to find hub genes associated with thyroid cancer microenvironment. The Cancer Genome Atlas (TCGA) database includes gene expression quantification data and clinical information of thyroid cancer. Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) website provides an easy access to predicting infiltration of immune cell and stromal cells in TME [8]. So, we combine these two datasets to screen out key genes in the thyroid cancer microenvironment.

Material and Methods

Data collection

Gene expression quantification data and associated clinical characteristics of thyroid cancer patients were downloaded from TCGA thyroid carcinoma (THCA) project (). The immune score and stromal score of thyroid cancer were retrieved from ESTIMATE (). Only samples with immune score and stromal score were included in our study.

Identification of differentially expressed genes

Based on median immune score as the cutoff criteria, we assigned thyroid cancer samples into a high immune score group and low immune score group. Using limma package [9] in R software (version 3.5.2), we identified the differentially expressed genes between high immune score group and low immune score group. Likewise, the differentially expressed genes between high stromal score group and low stromal score group were analyzed in the same method. The threshold of differentially expressed genes were: |log2fold change (log2FC)| >1 and false discovery rate (FDR) <0.05. Then, we identified the consensus differentially expressed genes in both groups. The results were exhibited using Venn plot.

Construction of protein-protein interaction (PPI) network

Using STRING database [10] (, version 11.0), we construct the PPI network of consensus gene. Gene-gene interactions with integrated scores bigger than 0.95 were selected and visualized using Cytoscape (version 3.7.1) [11]. The degree of genes indicated the number of edges (interactions) that linked to a given node (gene). We calculated the degree of each gene and selected the hub genes whose degree is beyond 12.

Functional enrichment analysis

Using clusterProfiler package [12], we conducted gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. GO terms of biological process (BP), cellular components (CC) and molecular functions (MF), and KEGG pathways with FDR <0.05 were considered to be significant.

Results

Stromal score and immune score of thyroid cancer samples and baseline data of thyroid cancer patients

A total of 501 patients with stromal score and immune score were included in our study. The stromal score and immune score of each sample were displayed in Supplementary Table 1. The median age was 46 years old (range: 15 to 89 years old). The majority of patients were female, and more than half of the patients were in stage I. Only a small number of patients developed distant metastases. The median follow-up time was 1.99 years (range: 0.09–14.86 years) (Supplementary Table 2). The detailed clinicopathologic characteristics of thyroid cancer patients were presented in Table 1.

Table 1

Clinicopathologic characteristics of thyroid cancer patients.

Characteristics	Number of case (%)
Age (years)
<45	228 (45.5)
≥45	273 (54.5)
Gender
Female	366 (73.1)
Male	135 (26.9)
Tumor stage
I	283 (56.5)
II	51 (10.2)
III	110 (22.0)
IV	55 (10.0)
NR	2 (0.4)
T stage
T1	142 (28.3)
T2	165 (32.9)
T3	169 (33.7)
T4	23 (4.6)
Tx	2 (0.4)
N stage
N0	226 (45.1)
N1	225 (44.9)
Nx	50 (10.0)
M stage
M0	279 (55.7)
M1	9 (1.8)
Mx	213 (42.5)

NR – not reported.

To pick out the differentially expressed genes with immune score and stromal score, we conducted differentially expression analysis. For differentially expressed genes with immune score, we identified 941 upregulated genes and 308 downregulated genes. 889 upregulated genes and 46 downregulated genes were associated with stromal score. The heatmap showed samples between high and low immune score group/stromal score group could be clearly distinguished by the differentially expressed genes (Figure 1A, 1B). The Venn plot showed the consensus gene with immune score and stromal score (Figure 1C, 1D). There were 793 consensus genes with immune score and stromal score (Supplementary Table 3).

Figure 1

Differentially expressed genes with immune score and stromal score. (A) Heatmap of the differentially expressed genes associated with immune score based on the median immune score. (B) Heatmap of the differentially expressed genes associated with stromal score based on the median stromal score. (C) Venn plot of downregulated consensus genes with immune score and stromal score. (D) Venn plot of upregulated consensus genes with immune score and stromal score.

Functional enrichment analysis of differentially expressed genes

To further investigate the biological functions of these consensus genes, we conducted functional enrichment analysis. The top 10 BP terms, CC terms and MF terms were exhibited in Figure 2A. T-cell activation, cell adhesion, cytokine activity, and extracellular matrix were the main enriched GO terms. Top 30 KEGG pathways were displayed in Figure 2B. The main enriched pathways included Th1 and Th2 cell differentiation, chemokine signaling pathway, cell adhesion molecules (CAMs), T cell receptor signaling pathway.

Figure 2

Functional enrichment analysis of differentially expressed genes. (A) Top 10 biological process (BP) terms, cellular components (CC) terms, molecular functions (MF) terms. (B) Top 30 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Protein-protein interaction (PPI) network and identification of hub genes

The PPI network of differentially expressed genes contained 296 nodes and 670 edges. We displayed the sub-network with most nodes and edges in Figure 3A. According to the PPI network, 30 genes (degree ≥12) were selected as hub genes (Figure 3B). Among these hub genes, 12 genes belonged to two families of chemokines named CXC and CC, 5 genes belonged to chemokine receptors, and two genes (IL10 and TNF) belonged to cytokines. Moreover, CXCL10 owed the highest degrees, suggesting that it may play crucial role in thyroid cancer microenvironment. Kaplan-Meier curve showed that low expression of CXCL10 was associated with poor overall survival of thyroid cancer patients (P-value of log-rank test=0.02467, Figure 3C).

Figure 3

Protein-protein interaction (PPI) network of differentially expressed genes and identification of hub gens. (A) PPI network of differentially expressed genes with integrated scores bigger than 0.95. (B) Hub genes with degree beyond 12. (C) Kaplan-Meier curve of CXCL10. High expression of CXCL10 was associated with good overall survival of thyroid cancer patients.

Functional enrichment analysis of hub genes

The PPI network of these 30 hub genes was displayed in Figure 4A. Biological process related to immune cells including Leukocyte migration and chemotaxis, myeloid leukocyte migration, and granulocyte chemotaxia were the mainly enriched BP terms (Figure 4B). The mainly enriched MF terms were cytokine activity and chemokine activity. As Figure 4C showed, these hub genes were mainly enriched in KEGG pathways that associated with immune response.

Figure 4

Protein-protein interaction (PPI) network of hub genes and functional enrichment analysis. (A) PPI network of hub genes. (B) Top 10 biological process (BP) terms, cellular components (CC) terms, molecular functions (MF) terms. (C) Top enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Discussion

Interactions between tumor cells and various components of TME are significant and contribute to all the hallmarks of cancer [13]. TME can affect how a tumor grows and spreads. Therefore, identification of key genes in thyroid cancer microenvironment aids appropriate management and treatment of thyroid cancer. In our study, we screened out 793 differentially expressed genes associated with immune score and stromal score. Then, these differentially expressed genes were demonstrated to relating to tumor microenvironment. Finally, we identified 30 hub genes by construction of PPI network. Of these hub genes, CXCL10 was positively associated with overall survival of thyroid cancer patients.

The functions of these differentially expressed genes were involved in Th1 and Th2 cell differentiation, T cell and lymphocyte cell activation, and T cell receptor signaling pathway, which is consistent with previous studies that immune cells are important components in TME and play a critical role in thyroid cancer invasion and metastasis [14]. By construction of a PPI network, we identified 30 hub genes in a thyroid cancer microenvironment. The functions of these hub genes are related to cytokines, chemokines and their receptors, which is in accordant with previous studies that cytokines, chemokines and their receptors are pivotal in building TME [4,15,16]. Chemokines and their receptors are involved in recruitment and maintenance of immune cells within the TME [4,17]. Through this process, they exert influence on the biological behavior of cancer. As reported, CXCR3 and CXCR4, 2 chemokine receptors, are highly expressed in thyroid cancer [18]. Moreover, CXCR3 is associated metastatic papillary thyroid cancer [19]. CXCR4, which is the most studied chemokine receptor, is involved in the diagnosis and treatment of thyroid cancer [20]. Among 12 chemokines, 4 chemokines named CXCL1, CCL5, CCL20, and CCL21 have been reported to influence several biological effects of thyroid cancer [17,21-24]. Cytokines such as interleukin (IL)-10 and tumor necrosis factor (TNF) have been reported to participate in many important events of carcinogenesis including epithelial mesenchymal transition, angiogenesis, invasion, and metastasis [16]. In addition, cytokines and chemokines are able to interact with each other. In papillary thyroid cancer, interferon (IFN)γ and TNF-α induce COCLÉ and CXCL11 secretion [25], which leads to the inhibition of papillary thyroid cancer proliferation and migration.

As CXCL10 is the highest interconnected node in the PPI network and positively associated with overall survival of thyroid cancer patients, it deserves more attention. CXCL10, CXC motif ligand 10 belongs to CXC family. Binding to CXCR3, it can lead to various effects, including stimulation of monocytes, natural killer and T-cell migration, and modulation of molecule expression. There is growing evidence for its antitumor effects in many kinds of cancer [26,27]. In vitro and in vivo studies have indicated that retroviral CXCL10 gene transduction to tumor cells led to reduction in tumor microvessel density and tumor growth in melanoma [28]. In lung cancer and sarcoma, in vitro and in vivo experiments showed that CXCL10 inhibit tumor angiogenesis, mitotic activity, and tumor growth [29]. CXCL10 exhibit similar antitumor effects in myeloma and glioma [30,31]. Moreover, CXCL10 can promote efficacies of other cytotoxic and targeted treatments, such as dacarbazine, temozolomide plus cisplatin for melanoma [32], all-trans retinoic acid (ARTA) for melanoma [33], and lapatinib with doxorubicin for breast cancer [34]. Importantly, CXCL10 can interact with other immune pathways and enhance them. For example, anti-programmed cell death-1 (PD-1) therapy could not only enhance T cell mediated tumor regression but also increase CXCL10 expression in breast cancer, colon cancer and melanoma [35,36]. COX-inhibitors increase CXCL10 expression and enhance the anti-tumor effect of CXCL10 in breast cancer and ovarian cancer [37,38]. However, the role of CXCL10 in thyroid cancer remains to be elucidated. Data from our studies indicated that CXCL10 may serve as protect gene in thyroid cancer, which is accordant with the tumor suppressing role of CXCL10 in other types of cancer.

Conclusions

Admittedly, our study has some limitations because it was based on the bioinformatics analysis only. The roles of these hub genes deserve further in vitro and in vivo studies because of their core roles in the PPI network. In conclusion, our study identified 30 hub genes in thyroid cancer microenvironment.