ATCs exhibit lower TDS values and higher infiltration levels of most immune cells than PTCs and normal thyroid tissues.
An overview of our study on thyroid dedifferentiation-associated IRGs is shown in Fig. 1. Loss of differentiation plays an important role in the progression of thyroid carcinoma. The TDS is a 16-metabolic-gene signature and serves as a parameter of thyroid differentiation [3]. We initially compared the TDS levels among ATCs, PTCs and normal thyroid tissues in the GSE33630 cohort. As expected, ATCs had obviously lower TDS values than the control samples for all genes (Fig. 2a-b). It has been reported that the immune landscape was correlated with the TDS in thyroid carcinoma [24]. To evaluate the differences in immune-related signatures, the ESTIMATE algorithm in the R program was applied to analyze transcriptome data and showed that ATCs had higher stromal scores and immune scores than PTCs and normal tissues (Fig. 2c). Consistent with this result, higher infiltration of most immune cells was observed in ATCs than in control groups, as indicated by the R package GSVA (Fig. 2d). Moreover, ATCs were well discriminated from PTCs or normal tissues with the PCA algorithm (Fig. 2e).
The inflammatory response is obviously activated in ATCs.
To further explore the differences in IRGs between ATCs and PTCs or normal thyroid tissues in the GSE33630 cohort, the powerful web tool ImmPort was used; a total of 1240 IRGs were present on our array. The volcano plots of GSEA displayed the top five pathways significantly activated in ATCs compared with normal tissues or PTCs, among which the inflammatory response was the most obvious (ATCs vs normal: normalized enrichment score (NES) = 2.16, p = 0.013; ATCs vs PTCs: NES = 2.15, p = 0.002) (Fig. 3a-b). We then analyzed the expression data of 1240 IRGs and found 207 IRGs were commonly deregulated (126 upregulated and 81 downregulated IRGs; fold change ≥ 2 and p<0.05) in ATCs compared with both normal tissues and PTCs (Fig. 3c-d). GO and KEGG enrichment analysis of consistently deregulated IRGs indicated significant activation of immune-associated signaling pathways in ATCs (Fig. 3e-f). To further screen the deregulated IRGs, another GEO cohort GSE29265 was also included, and a total of 68 genes (36 upregulated IRGs and 32 downregulated IRGs) were screened out and validated with the two cohorts (Fig. 4a-c).
Identification of differentiation-associated and prognostic IRGs.
To further analyze the differentiation-associated IRGs, we initially enrolled 505 PTC samples with complete clinical annotations and transcriptome data from TCGA database, which were then divided into the low-differentiated (low TDS) group and the high-differentiated (high TDS) group by the median TDS value. The volcano plot based on GSEA of 1240 IRGs indicated that the two phenotypes displayed distinct activation of immune-related pathways. Compared with the high-differentiated group, interferon gamma response (NES = 2.15, p = 0.003), inflammatory response (NES = 2.10, p = 0.003) and interferon alpha response (NES = 1.90, p = 0.002) were significantly enriched in the low-differentiated group (Additional file 1: Figure S1). Considering that patients with ATC in GEO database lack clinical information, we thus focused on the low-differentiated group in the TCGA cohort to identify the IRGs relative to prognosis. First, clinical characteristics and 68 IRGs were enrolled in a univariate Cox analysis. Stage III/IV, T3/T4, distant metastasis and the expression of 5 genes, including matrix metalloproteinase-9 (MMP9), fibroblast growth factor receptor 2 (FGFR2), plasminogen activator, urokinase receptor (PLAUR), syndecan-2 (SDC2) and thyroglobulin (TG), were identified as risk factors for the PFI of low-differentiated PTCs in the TCGA cohort (Table 1). In addition, multivariate Cox regression analysis was carried out to correct for confounding factors. As a result, SDC2 (HR = 0.614; 95% CI: 0.391-0.962; p = 0.033) and MMP9 (HR = 1.339; 95% CI: 1.076-1.666; p = 0. 009) were determined as independent prognostic factors of low-differentiated PTCs (Table 2).
We then assessed the expression levels of MMP9 and SDC2 in a combined GEO cohort composed of 52 ATCs, 78 normal tissues and 69 PTCs from the same chip platform and in a TCGA cohort grouped by TDS level. As shown in Fig. 5a and 5b, MMP9 was significantly higher in ATCs and low-differentiated PTCs than in normal tissues and high-differentiated PTCs, whereas SDC2 displayed the opposite expression profile, as its expression was relatively low. qPCR analyses further indicated that ATC cell lines CAL-62, DRO and 8305C possessed lower SDC2 levels; DRO and 8305C had higher MMP9 levels than PTC cell line TPC-1 and normal thyroid cell line Nthy-ori 3-1 (Fig. 5c). In addition, an immune-related risk score was constructed and was obviously higher in ATCs and low-differentiated PTCs than in normal tissues and high-differentiated PTCs (Fig. 5d).
We proceeded to analyze the ability of the risk score to predict dedifferentiation in the combined GEO cohort and TCGA cohort. Correlation analyses (Fig. 5e) showed that the risk score was negatively associated with TDS (r = -0.77, p < 0.001 in the combined GEO cohort; r = -0.46, p < 0.001 in the TCGA cohort). Consistent with this finding, ROC curves confirmed the robust predictive value in both cohorts, and the AUC values were 0.842 (p < 0.001) and 0.707 (p < 0.001) (Fig. 5f). Moreover, we estimated the prognostic value of a single gene and the two-gene risk score. Kaplan-Meier curves (Fig. 5g) indicated that PTC patients with high MMP9 (HR = 2.28; p = 0.005), low SDC2 (HR = 2.47; p = 0.002) or high risk score (HR = 3.05; p = 0.0003) had a significantly shorter PFI than other patients in the TCGA cohort, indicating the potential of the risk score to serve as a prognostic biomarker.
The two-gene risk score is closely correlated with immune-related signatures.
Immunotherapy has revolutionized the treatment of many patients with cancer. However, through exploration of the complicated mechanisms of cancer immunoediting, it was found that some negative costimulatory molecules and immunosuppressive cytokines could promote immune escape of tumor cells by inhibiting immune responses, especially antitumor T cells [25, 26], which stimulated the appearance of immune checkpoint therapy. We wanted to explore the role of the two-gene risk score in immune-related signatures. First, in both the combined GEO cohort and TCGA cohort, a marked positive correlation was found between the risk score and important immune checkpoint molecules, including programmed cell death ligand 1 (PDL1; also known as CD274), cytotoxic T-lymphocyte-associated protein 4 (CTLA4), indoleamine 2,3-dioxygenase 1 (IDO1) and hepatitis A virus cellular receptor 2 (HAVCR2; also known as TIM3), as illustrated in Fig. 6a and 6b. We further analyzed whether the risk score was correlated with the infiltration of immune cells; ssGSEA was applied, and as shown in Fig. 6c, most immune cells displayed a higher level of infiltration in the high-risk score group than in the low-risk score group in the combined GEO cohort. Correlation analysis further confirmed the close relationship between risk score and the infiltration of most immune cell types (Fig. 6d).
With deep research on predictive markers for immunotherapy, emerging studies support the crucial role of the TME in the response to immune checkpoint inhibitors [8, 27]. Considering the positive correlation of the risk score with the expression of common immune checkpoint molecules and the infiltration of immune cells, we further explored whether it helps to predict the efficacy of immunotherapy. At present, increasing researches have demonstrated that some immune checkpoint inhibitors, including the PD-1 inhibitors spartalizumab and pembrolizumab and the PD-L1 inhibitor nivolumab [2, 12-15], are effective in ATCs as a salvage therapy. However, no such cohort with complete prognosis data and expression profiles was obtained for ATCs. We therefore investigated the value of the two-gene risk score for predicting the response to immunotherapy in the IMvigor210 cohort and the GSE63557 cohort, a dataset of urothelial cancer patients who received an anti-PD-L1 agent and a dataset of AB1-HA mesothelioma tumor mice treated with anti-CTLA4 therapy [19, 20]. We found that patients with high risk score tended to have a better response to immunotherapy than those with a low risk score in the IMvigor210 cohort (p = 0.003; Fig. 6e). The number of mice responding to anti-CTLA4 agents was higher in the GSE63557 cohort, but the difference was not statistically significant (p = 0.18; Fig. 6g). In addition, the risk score was higher in patients with complete response (CR) to anti-PD-L1 therapy than in patients with progressive disease (PD) (p < 0.05) (Fig. 6f). Similarly, the risk score was slightly increased in mice that responded to anti-CTLA4 therapy compared with those that did not respond (p < 0.05) (Fig. 6h).
Molecular basis of the immune-related risk score signature.
We further explored the differences between the high-risk score and low-risk score samples. GSEA was performed on the broad hallmark gene sets in the combined GEO cohort. The volcano plots displayed distinct enrichment of signaling pathways. As shown in the GSEA plots with high NES values, compared with the low-risk score group, epithelial-mesenchymal transition, TNFα signaling, and some common immune-related signaling pathways, including the IL-6/JAK/STAT3 pathway, interferon alpha response, interferon gamma response and inflammatory response, were significantly enriched in the high-risk score group (Additional file 2: Figure S2a-b).