High Expression of THBS1 Leads to a Poor Prognosis in Papillary Thyroid Cancer and Suppresses the Anti-Tumor Immune Microenvironment.

Objectives: To study the role of thrombospondin-1 (THBS1) in papillary thyroid cancer (PTC) prognosis and the immune microenvironment.
Methods: A retrospective cohort study was designed, and data from The Cancer Genome Atlas database and PTC tissues from Fudan University Shanghai Cancer Center were used. Weighted gene co-expression network analysis was performed to build a THBS1-immune-related gene prognostic index (T-I index).
Results: High THBS1 expression was correlated with advanced TNM stage, higher recurrence risk, and shorter progression-free interval. High THBS1 expression correlated with MAPK and PD1 pathways indicating a tumor promoting and immunity-inhibiting tendency. The T-I index showed a powerful capacity to predict progression-free survival and immunotherapy benefit.
Conclusion: High expression of THBS1 leads to a poor prognosis in PTCs and suppresses the anti-tumor immune microenvironment.

Introduction

Thyroid cancer is the most common malignant tumor of the endocrine system, and the incidence of thyroid cancer, especially papillary cancer, has been increasing. Most papillary thyroid cancers (PTCs) can be surgically removed and show a good overall prognosis and low fatality rate. However, a considerable number of patients show cervical lymph node metastasis and need to undergo cervical lymph node dissection surgery, which seriously impacts the postoperative quality of life. Advanced PTCs, such as tumors with a large cervical mass and extensive extrathyroidal invasion of tumors, are more difficult to operate on, and postoperative adverse reactions can occur, such as cervical recurrent laryngeal nerve injury. Therefore, the identification of molecular markers of advanced PTCs and the corresponding therapeutic targets is critical improving patient treatment.

The tumor microenvironment (TME) plays an important role in all stages of tumor development and contains multiple components, including tumor cells, tumor stromal cells, and immune cells. Immune cells within the TME have key functions in tumorigenesis, with both tumor-promoting and tumor-suppressive roles. Recent years have seen the development of immunotherapy as cancer treatments that function by activating innate immunity or blocking immune checkpoints. Although a large number of studies have shown that immunotherapy has beneficial effect on a variety of tumors,[7,8] the research on immunotherapy for thyroid cancer, especially for advanced PTCs, is limited.[9,10] Given the limited response to immunotherapy, it is important to identify immunotherapy-related targets and patients who can benefit from immunotherapy.

Thrombospondin-1 (THBS1) is a member of the thrombospondin protein family. It is a large matricellular glycoprotein, with various protein-binding domains. THBS1 has multiple biological functions, such as in wound repair and tissue generation. Recent studies have shown that THBS1 is also involved in cancer development. THBS1 plays 2 major roles in tumors. THBS1 inhibits neovascularization in tumors[13,14] but also promotes tumor invasion and metastasis. This process is regulated by different pathways in different tumors. For example, in breast tumors, THBS1 regulates tumor cell adhesion and invasion by upregulating MMP-9 or TGF-β or activating the urokinase plasminogen system, thus promoting tumor invasion and metastasis, as observed in breast cancer and thyroid cancer. In our imaging histology research on cervical lymph node metastasis ultrasound images of thyroid cancer, high THBS1 expression negatively correlated with lymph node metastasis, which is inconsistent with the previous study. However, our imaging histology research mainly focused on ultrasound images of PTCs and did not explore the function of THBS1 in tumors. We also noted a correlation between the gene module with THBS1 as the hub gene and tumor immunity. In view of studies showing that THBS1 has a role in suppressing the antitumor immune microenvironment in gastric cancer, we believe that the role of THBS1 in PTC and the role it plays in the tumor immune microenvironment is worthy of investigation.

In this study, we performed a preliminary analysis on the expression level and possible role of THBS1 in PTCs using data from a public database and clinical specimens that were retrospectively collected. We further designed a THBS1-immune-related gene prognostic index (T-I index) based on THBS1 and immune-related genes and studied its efficacy as a predictor of thyroid cancer prognosis and immunotherapy response.

Methods and Materials

Patient Data Acquisition

The reporting of this study conforms to STROBE guidelines. Thyroid carcinoma patient datasets, with gene expression profiles and clinical information, were downloaded from the publicly available The Cancer Genome Atlas TCGA-THCA project. This project included 507 cases with 510 tumor samples and 58 paired normal samples. We applied the following exclusion criteria in the study analyses: (1) samples with incomplete clinical data were excluded from the subgroup comparison, and (2) in certain subgroup analyses, normal tissue samples were excluded. After data sorting, RNA sequencing data of 503 tumor samples and 58 paired normal samples were converted to RNA seq data in transcript per million (TPM) format.

For the validation cohort, 53 patients with paired PTC and normal tissues were retrospectively enrolled from Fudan University Shanghai Cancer Center (FUSCC). We obtained frozen tumor tissues and the related pathological data from the Department of Biobank and Pathology at the FUSCC. The inclusion criteria were as follows: (1) tissue specimens were completely preserved in liquid nitrogen or at −80 °C, with sufficient sample size for RNA extraction and sectioning, (2) one or more follow-up visits were recorded, and (3) the medical history and examination data were complete.

All research protocols were approved by the ethical committee of FUSCC (approval 2101-ZZK-41) and all enrolled patients signed informed consent forms.

Survival and Clinical Correlation Analysis

To maintain consistency of grouping criteria within the same study, TNM stage and extrathyroidal extension (ETE) were defined according to the 7th edition of American Joint Committee on Cancer guidelines. In this research, ETE includes both cases of tumor breaking through the perithelium and invasion by extra-thyroidal peripheral adipose tissue; this type of tumor is classified as T3. Since thyroid carcinoma patients show a good overall survival, we selected progression-free survival (PFS) as the prognostic end point. The gene expression values were grouped according to the best cut-off values in Kaplan–Meier analysis with R package “survival.” Cox regression was used to analyze the association between PFS and clinicopathologic characteristics in TCGA and the FUSCC cohort. Risk stratification of thyroid cancer recurrence was classified into low, moderate, and high grades according to the American Thyroid Association Guidelines and used as a prognostic indicator of FUSCC cohort. Receiver operating characteristic (ROC) curve was used to analyze the ability of THBS1 to predict high recurrence risk.

Quantitative Real-Time PCR (qPCR) and Immunohistochemistry

qPCR was performed on the 53 PTC tissues in the FUSCC cohort, and the 2−ΔΔCT method was used to determine the relative expression of THBS1. For immunohistochemistry (IHC), paraffin sections were prepared from frozen tissue. The Thrombospondin-1 antibody (A6.1, catalogue: NB100-2059, 1:100 dilution) was purchased from Novus Biologicals (Bio-Techne China Co., China). Specimens were scored according to the intensity of staining and the number of positive cells in 5 high magnification views selected at random, and an immunohistochemical score was calculated for each specimen by adding the scores. See Supplemental Appendix 1 for detailed experimental procedures.

Gene Set Enrichment Analysis (GSEA)

Differential gene screening was performed according to THBS1 expression level using R package “DESeq2” (|log2(FC)| >1, P.adjust <.05). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for differentially expressed genes, and GSEA was performed using the R package “clusterProfile.” Significance was defined as a false discovery rate (FDR) <0.25 and P.adjust <.05.

Establishment of the T-I index

The immune-related gene lists were obtained from ImmPort, and after intersection with former identified THBS1-related differentially expressed genes, THBS1-immune-related genes were obtained.

Weighted gene co-expression network analysis (WGCNA) was used to identify hub genes with the R package “WGCNA.” First, RNA seq data of 503 TCGA tumor samples were used to calculate the Pearson correlation coefficient between 2 genes, and the similarity matrix was constructed. Next, the similarity matrix was transformed into an adjacency matrix with a network type of signed and a soft threshold of β = 3. The adjacency matrix was converted into a topological overlap matrix to reduce noise and false correlation, and the new distance matrix was obtained. After building the dynamic pruning tree to identify the modules, and setting the module membership (MM) cut-off criteria as |MM| > 0.8, hub genes were obtained. Genes significantly affecting PFS were identified by Kaplan–Meier analysis and multivariate Cox regression analysis along with THBS1. The T-I index score of each sample was calculated by multiplying the expression values of genes and adding the scores together. The coefficient of each gene was determined by its weight in the Cox model. Chord diagram was used to show the relationship between THBS1 and other genes in T-I index by R package “Circhize.”

Value of T-I index in Prognosis and Immunotherapy Benefit Prediction

The prognosis of different T-I subgroups was evaluated by a nomogram and Kaplan–Meier curves in TCGA cohorts. To validate the independent prognostic value of T-I index, univariate and multivariate Cox regression analyses were performed. The Tumor Immune Dysfunction and Exclusion (TIDE) score of each sample in TCGA cohort was calculated online (http://tide.dfci.harvard.edu/) to predict the likelihood that patients will benefit from immune checkpoint inhibitors (ICI) therapy.[28,29] ROC curve was used to compare T-I index and single THBS1 expression level to predict the ability of ICI therapy responder.

Molecular and Immune Characteristics Analysis in Different Expression and T-I Subgroups

To get a more complete insight of the index, we investigated the genetic mutations and immune landscape associated with the index. To analyze the genes mutated in different T-I index subgroups, information on genetic alterations was obtained from the cBioPortal database (http://www.cbioportal.org/). The R package “estimate” was used to calculate the Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) score of each sample. Single sample GSEA (ssGSEA) was used to analyze the relationship between THBS1 and 24 classic tumor immune cell subtypes (R package “GSVA”). In addition, the CIBERSORTX (https://cibersortx.stanford.edu/) website and the LM22 signature were used to calculate the tumor immune cell infiltration score of each case from TCGA cohort. We then used TIMER (Tumor Immune Estimation Resource) website (http://timer.cistrome.org/)[30-32] and the clinical data from TCGA cohort to screen immune cell subtypes that have a significant impact on the prognosis of thyroid cancer. Spearman was used to analyze the correlation between THBS1 expression and these immune cells.

Immunofluorescence (IF) staining was performed to determine the infiltration of immune cells in the PTC microenvironment. The CD4 (catalogue: GB13064-1, 1:100 dilution) and FoxP3 (catalogue: GB11093, 1:100 dilution) antibodies were purchased from Servicebio (Servicebio Co., China).

Statistical Analysis

Statistical data acquired from TCGA were merged and converted by R-3.6.3. A P value <.05 was set as the cut-off criterion. The R package “pheatmap” was used to draw a tumor-infiltrating immune cell heatmap. The R package “ggplot2” was used for data visualization and image rendering. The Mann–Whitney U-test, Kruskal–Wallis H–test, and Dunn's test were used for nonparametric tests of independent samples. Wilcoxon signed rank test was used for nonparametric test of paired samples.

Results

THBS1 Expression Level is Correlated with Invasion, Metastasis, and Poor Prognosis of Thyroid Cancer

We obtained clinical and gene expression data of 503 thyroid carcinoma cases, including 503 tumor samples and 58 paired normal samples, from TCGA-THCA project. Figure 1 shows a summary of the overall analysis performed in this research. The baseline data of all cases are listed in Supplemental Appendix 2. The THBS1 expression level in the tumor group was significantly higher than that in the normal group in unpaired samples (U = 0.631, P = .005, Figure 2A and B). However, no significant difference in THBS1 expression level was detected in the 58 paired samples (Supplemental Appendix 3). There were 357 classical PTC (C-PTC), 101 follicular variant PTC (FV-PTC), 36 tall cell variant PTC (TCV-PTC), and 9 other histological type cases. THBS1 expression in FV-PTC was lowest, and its expression in TCV-PTC was highest (P < .001, Figure 2C). THBS1 expression in the T3&T4 group was significantly higher than that in the T1&T2 group (U = 0.524, P < .001, Figure 2D). In addition, THBS1 expression was higher in patients with lymph node metastasis than in patients without lymph node metastasis (U = 0.669, P < .001, Figure 2E). The expression of THBS1 in the pathological stage III&IV group was significantly higher than that in the stage I&II group (U = 0.45, P = .001, Figure 2F). Furthermore, THBS1 expression was significantly higher in patients with ETE than in those patients without ETE (U = 0.727, P < .001, Figure 2G).

Figure 1.

Graphic abstract of this study. Abbreviations: TCGA, The Cancer Genome Atlas; FUSCC, Fudan University Shanghai Cancer Center; T-I index, THBS1-immune-related gene prognostic index.

Figure 2.

THBS1 was highly expressed in advanced thyroid cancer. Data in Figure (A) to (H) were from TCGA cohort. (A) Expression of THBS1 in normal and tumor tissues in different cancer types. (B) The expression level of THBS1 in the tumor group was higher than that in the normal group. (C) The expression on THBS1 in FV-PTC was lowest, and it in TCV-PTC was highest. (D) The expression level of THBS1 in the T1&T2 group was lower than that in the T3&T4 group. (E) THBS1 expression was higher in patients with lymph node metastasis than in patients without lymph node metastasis. (F) The expression level of THBS1 in the TNM stage I&II group was lower than that in the stage III&IV group. (G) THBS1 expression was higher in patients with extrathyroidal extension* than in those who without. (H) Kaplan–Meier survival analysis of THBS1 expression groups. Data in Figure (I) to (P) were from FUSCC cohort. (I) Tumors with high THBS1 expression were larger than with low THBS1 expression. (J) Patients with high THBS1 expression got more metastatic lymph nodes in the neck. (K) The expression level of THBS1 in the TNM stage I&II group was lower than that in the III&IV group. (L) THBS1 expression was higher in patients with extrathyroidal extension* than in those who without. (M) ROC curve for THBS1 expression level to predict high risk of recurrence. (N to O) Representative images of immunohistochemical staining of THBS1 protein in normal thyroid and papillary thyroid carcinoma tissue, scale bar: 100 μm. (P) Comparison of immunohistochemical scores of normal-tumor paired samples. Abbreviations: TCGA, The Cancer Genome Atlas; FUSCC, Fudan University Shanghai Cancer Center; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LGG, brain lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; THCA, thyroid carcinoma; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma; C-PTC, classical papillary thyroid carcinoma; FV-PTC, follicular variant papillary thyroid carcinoma; TCV-PTC, tall cell variant papillary thyroid carcinoma. *, extrathyroidal extension was defined according to 7th American Joint Committee on Cancer guidelines. *, P < .05; **, P < .01; ***, P < .001; ns: not significant.

The cut-off value of THBS1 expression with the smallest P value was selected and patients were categorized into high and low expression groups using this value. Kaplan–Meier curve analysis revealed that the prognosis of the high THBS1 expression group (n = 323) was significantly worse than that of the low THBS1 expression group (P = .046, Figure 2H). However, multivariate Cox regression analysis showed that THBS1 expression level was not an independent risk factor for the prognosis of thyroid cancer (P = 0.125, Table 1).

Table 1.

Univariate and Multivariate Cox Regression Analysis of THBS1 Expression and Thyroid Papillary Cancer Progression-free Survival.

Characteristics	Total (N)	Univariate analysis		Multivariate analysis
		Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
Age at initial pathologic diagnosis	503	1.019 (1.003-1.037)	.023	1.018 (1.000-1.037)	.055
Gender	503
Female	368	Reference
Male	135	1.535 (0.866-2.719)	.142
Tumor volume (mm3)	403	1.011 (1.001-1.022)	.039	1.005 (0.991-1.019)	.476
Primary neoplasm focus type	493
Unifocal	266	Reference
Multifocal	227	1.028 (0.591-1.788)	.923
Primary thyroid gland neoplasm location	497
Right lobe	213	Reference
Left lobe	176	1.057 (0.573-1.951)	.859
Bilateral	86	1.145 (0.523-2.505)	.735
Isthmus	22	0.409 (0.055-3.044)	.383
Histological type	503
C-PTC	357	Reference
FV-PTC	101	0.601 (0.254-1.422)	.246
TCV-PTC	36	2.114 (0.942-4.741)	.069
Other	9	1.060 (0.145-7.727)	.954
Pathologic stage (T stage)	503
T1&2	308	Reference
T3&4	193	2.520 (1.441-4.407)	.001	1.528 (0.585-3.992)	.387
Tx	2	0.000 (0.000-Inf)	.996	0.000 (0.000-Inf)	.998
Pathologic stage (N stage)	503
N1	224	Reference
N0	229	0.637 (0.357-1.137)	.127
Nx	50	0.585 (0.206-1.665)	.315
Pathologic stage (M stage)	502
M0	282	Reference
M1	9	7.767 (2.956-20.412)	<.001	3.731 (1.010-13.790)	.048
Mx	211	1.300 (0.734-2.305)	.368	1.587 (0.847-2.974)	.149
With ETE	485
No	332	Reference
Yes	153	2.004 (1.157-3.473)	.013	0.964 (0.385-2.416)	.938
THBS1 expression group	503
Low	180	Reference
High	323	1.833 (1.06-3.494)	.046	1.735 (0.859-3.506)	.125

C-PTC, classical thyroid papillary cancer; FV-PTC, follicular variant thyroid papillary cancer; TCV-PTC, tall cell variant thyroid papillary cancer; ETE, extrathyroidal extension ETE. Variables with P values less than 0.1 in the univariate regression analysis were included in the multivariate regression equation. P values less than .05 are indicated in bold.

We next examined the expression of THBS1 in the 53 PTC tissue samples in the validation cohort from FUSCC (Supplemental Appendix 4). Tumors with high THBS1 expression showed a larger tumor size (U = 7, P = .011, Figure 2I) and more numbers of lymph node metastases than the low THBS1 expression group (U = 4 P = .002, Figure 2J). We also observed that with the increase of TNM grade, THBS1 expression level also increased (U = 3.928, P < .001, Figure 2K). In addition, THBS1 expression in the samples with ETE was higher than that in the other group (U = 2.285, P = .004, Figure 2L). There were 32 patients at low or moderate recurrence risk, while 21 patients at high recurrence risk. The sensitivity of THBS1 expression to predict high recurrence risk was 1, and the specificity was 0.719 (area under ROC curve = 0.821, Figure 2M). Immunohistochemical staining also showed that THBS1 protein expression was higher in tumor tissue with a large number of lymph node metastases and ETE (Figure 2N to O). The IHC score of tumor samples was higher than paired normal samples (Figure 2P, Z = 1, P = .007).

MAPK, Tumor Adhesion, and Immune-Related Pathways Enriched in the THBS1 High Expression Phenotype

There were 1826 genes that were differentially expressed according to THBS1 expression level, including 1320 genes highly expressed in high THBS1-expression group and 506 genes expressed at low levels in the low THBS1-expression group (Figure 3A). GSEA of data from thyroid cancers with low THBS1 and high THBS1 expression was used to identify THBS1-related signaling pathways. There were 432 datasets with an FDR (q value) <0.25 and P.adj <.05. Because of space limitations, we selected several pathways associated with high and low expression to display in Figure 3B. In addition to the classical thyroid cancer-related MAPK, NF-kB pathway and various tumor-related pathways, THBS1 expression levels were also associated with many immune-related pathways, such as REACTOME_IMMUNOREGULATORY_INTERACTIONS_BETWEEN_A_LYMPHOID_AND_A_NON_LYMPHOID_CELL and REACTOME_PD_1_SIGNALING pathway. KEGG and GO enrichment analyses were also performed, and several pathways are shown in Figure 3C (the complete GO results are listed in Supplemental Appendix 5). Based on the conditions of P.adj <.05 and q value <0.2, there were 14 pathways in KEGG, 179 pathways in biological process (BP), 14 pathways in cellular component (CC), and 38 pathways in molecular function (MF) subgroups. Similar to the results of GSEA, KEGG and GO analyses also revealed a close relationship between high THBS1 expression and tumor immunity. The Wnt signaling pathway was found in the enrichment analysis of KEGG, GO and MF. Furthermore, cell adhesion, integrin binding, and other pathways were related to high THBS1 expression.

Figure 3.

THBS1-related signaling pathways based on GSEA and GO enrichment analysis. (A) Differentially expressed genes associated with THBS1 expression. 1826 genes met |log2 (FC)| >1 and P.adj < .05 threshold value including 1320 high expression genes (log FC positive, red points) and 506 low expression genes (log FC negative, blue points). (B) Four noteworthy signaling pathways THBS1-related signaling pathways based on GSEA. (C) 10 noteworthy signaling pathways THBS1-related signaling pathways based on KEGG, MF, and BP.

High Expression of THBS1 Suppresses the Anti-Tumor Immune Microenvironment

The enrichment analysis results suggested that THBS1 expression level may be related to immune molecules and immune pathways. Therefore, we examined whether THBS1 expression was associated with immune infiltration in thyroid cancer. We first calculated the ESTIMATE score of each sample to reflect the tumor immunity score and tumor purity. As shown in Figure 4A, THBS1 expression was positively correlated with both stromal score and immune score, indicating that patients with higher THBS1 expression level had more stromal cells and immune cells infiltrated in TME (P < .001). We then used ssGSEA to analyze the relationship between THBS1 expression and 24 classic tumor immune cell subtypes (Figure 4B). We found that THBS1 expression was positively correlated with the level of infiltration of most of the immune cells. We also used CIBERSORTx to analyze the level of immune cell infiltration of each sample in TCGA cohort with the LM22 signature for classification. The proportions of 22 types of immune cells in different THBS1 expression groups are shown in Figure 4C. Based on the results of these 2 analyses, we concluded that THBS1 expression level was correlated with T cell infiltration level of different subtypes. Therefore, we used TIMER website and clinical data from TCGA cohort to screen immune cells that could affect the prognosis of thyroid cancer. We found that high fractions of follicular helper T cells (TFH) and regulatory T cells (Tregs) were significantly associated with shorter progression-free interval (PFI) in thyroid cancer (TFH: HR = 3.02, P = 0.011, Tregs: HR = 2.88, P = .009, Figures 4D and E). Furthermore, Figures 4F and G show the positive relationship between THBS1 expression level and TFH and Tregs infiltration, indicating that high expression of THBS1 might lead to a poorer outcome in thyroid cancer via suppressing the anti-tumor immune microenvironment.

Figure 4.

Relationship between THBS1 expression and tumor-infiltrating immune cells. (A) Tumor immunity score according to THBS1 expression group using ESTIMATE algorithm. Correlation analysis between THBS1 expression and ESTIMATE score. Tumors with high THBS1 expression got higher stromal, immune, and ESTIMATE score (P < .001). (B) Relationship between THBS1 expression and 24 subtypes of tumor-infiltrating immune cells according to ssGSEA method. (C) The infiltrating level of immune cells in different THBS1 expression groups according to CIBERSORTx and LM22 signature. (D to E). Kaplan–Meier survival analysis of TFH and Tregs infiltrating level groups (HR = 3.02, P = 0.011, HR = 2.88, P = .009). (F to G). Correlation analysis results between THBS1 expression and immune cell infiltration fraction according to ssGSEA: TFH (r = 0.27, P < .001); Tregs (r = 0.380, P < .001). ssGSEA, single sample gene set enrichment analysis; ESTIMATE, estimation of stromal and immune cells in malignant tumors using expression data; TFH, follicular helper T cells, Tregs, regulatory T cells. ns, not significant. *, P < .05; **, P < .01; ***, P < .001.

We performed IF staining of tissue sections from patients in the FUCSS cohort to show Tregs infiltration in papillary thyroid carcinoma. As shown in Figure 5, in samples with high THBS1 expression, there was a high infiltration of CD4-positive/FoxP3-positive Tregs in the tumor tissue.

Figure 5.

Immunofluorescence staining for Tregs infiltration in high THBS1-expressing samples. (A) Representative image of immunohistochemical staining of THBS1 protein in a high THBS1-expression papillary thyroid carcinoma tissue, scale bar: 100 μm. (B to F) Immunofluorescence staining of Tregs from the same sample (scale bar in B to E: 20 μm, scale bar in F: 10 μm). Arrows point to CD4-positive/FoxP3-positive co-localized Tregs.

THBS1-Immune-Related Hub Genes and a Prognostic index

Since high THBS1 expression level is not an independent risk factor for PFS in thyroid cancer, we then designed a THBS1-immune-related gene prognostic index. We examined differentially expressed genes between 323 high THBS1 expression cases and 180 low THBS1 expression cases, and a total of 1826 differentially expressed genes were obtained. After intersecting these genes with the list of immune-related genes acquired from ImmPort, a total of 292 THBS1-immune-related genes were selected for WGCNA (Figure 6A). The optimal soft-thresholding power was 3 based on the scale-free network (Supplemental Appendix 6A-B). A total of 75 hub genes and 3 modules were allocated (30 genes in module blue, 5 genes in module brown, 40 genes in module turquoise, 1 gene in module grey, Supplemental Appendix 6C-D). Because of the small number of hub genes screened, we did not conduct phenotypic correlation analysis on these 3 modules. However, we did perform Kaplan–Meier analysis to identify the genes that correlated with the PFS of thyroid cancer. Twelve genes were significantly associated with PFS, and their expression levels were grouped according to the optimal cut-off value. Multivariate Cox regression analysis was performed, and 4 genes along with THBS1 were identified for index construction with their coefficient in the Cox model. The THBS1-immune-related gene prognostic index was calculated by the formula = expression level of IGHV3-49 *0.18 + expression of CXCL13 *0.006 + TRAV8-3 *0.3 + expression of TRBV30 *0.22 + expression of THBS1 *0.45. Figure 6B shows the co-expression heatmap of THBS1 and these 4 genes, and Figure 6C shows the relationship between them. A flow chart showing the establishment of the THBS1-immune-related gene prognostic index is depicted in Supplemental Appendix 6E, and Kaplan–Meier curves of 4 selected genes are shown in Figure 6D-G.

Figure 6.

Design of THBS1-immune-related gene prognostic index and K-M plots of 4 selected genes. (A) Venn diagram of 1826 THBS1-related genes and 1793 immune-related genes with 292 genes in intersection region. (B) Co-expression heatmap of THBS1 and other 4 T-I index genes obtained from WGCNA (Spearman analysis, P < .001). (C) Chod diagram of THBS1 and other 4 T-I index genes. The red circle arc represented 2 genes were positively correlated. (D to G) Kaplan–Meier survival analysis of 4 selected genes significant in the univariate Cox analysis (P ≤ .05). T-I index, THBS1-immune-related gene prognostic index; WGCNA, weighted gene co-expression network analysis.

Univariate and multivariate Cox analyses indicated that high T-I index score was an independent risk factor for shorter PFI in thyroid cancer (hazard ratio = 2.738, P = .025, Table 2). A nomogram of clinical characteristics and T-I index group is shown Figure 7A. Patients in the high T-I subgroup had a significantly shorter PFI than those in the low T-I subgroup (P = .035, Figure 7B). TIDE website provides an online ICI therapy benefit prediction algorithm including dysfunction and exclusion score to reflect the dysfunction and exclusion of T cells. Higher TIDE score represented a higher potential for immune evasion, which suggested that the patients were less likely to benefit from ICI therapy. According to the Wilcoxon test, both dysfunction and exclusion scores in the high T-I index subgroup were higher than those in the low subgroup (Figure 7C), indicating that patients with a higher T-I index score may be more likely to benefit from ICI therapy. According to TIDE results, there were 426 potential responders to ICI therapy (Supplemental Appendix 7). The distribution of different histological subtypes of PTC in T-I index groups and TIDE responder groups is shown in the Appendix. Compared with THBS1 expression alone, the T-I index exhibited a better ability to screen these responders (Figure 7D).

Table 2.

Univariate and Multivariate Cox Regression Analysis of T-I Index and Thyroid Papillary Cancer Progression-free Survival.

Characteristics	Total(N)	Univariate analysis		Multivariate analysis
		Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
Age at initial pathologic diagnosis	503	1.019 (1.003-1.037)	.023	1.019 (1.001-1.038)	.043
Gender	503
Male	135	Reference
Female	368	0.651 (0.368-1.154)	.142
Tumor volume (mm3)	403	1.011 (1.001-1.022)	.039	1.006 (0.992-1.019)	.410
Primary neoplasm focus type	493
Unifocal	266	Reference
Multifocal	227	1.028 (0.591-1.788)	.923
Primary thyroid gland neoplasm location	497
Right lobe	213	Reference
Left lobe	176	1.057 (0.573-1.951)	.859
Bilateral	86	1.145 (0.523-2.505)	.735
Isthmus	22	0.409 (0.055-3.044)	.383
Histological type	503
C-PTC	357	Reference
FV-PTC	101	0.601 (0.254-1.422)	.246
TCV-PTC	36	2.114 (0.942-4.741)	.069
Other	9	1.060 (0.145-7.727)	.954
Pathologic stage (t stage)	503
T1&2	308	Reference
T3&4	193	2.520 (1.441-4.407)	.001	1.550 (0.590-4.074)	.374
Tx	2	0.000 (0.000-Inf)	.996	0.000 (0.000-Inf)	.998
Pathologic stage (n stage)	503
N0	229	Reference
N1	224	1.569 (0.880-2.798)	.127
Nx	50	0.918 (0.312-2.700)	.877
Pathologic stage (m stage)	502
M0	282	Reference
M1	9	7.767 (2.956-20.412)	<.001	4.264 (1.162-15.652)	.029
Mx	211	1.300 (0.734-2.305)	.368	1.648 (0.882-3.082)	.118
With ETE	486
No	332	Reference
Yes	153	2.004 (1.157-3.473)	.013	0.946 (0.377-2.372)	.905
T-I index	503
Low	134	Reference
High	369	2.353 (1.061-5.217)	.035	2.738 (1.133-6.616)	.025

C-PTC, classical thyroid papillary cancer; FV-PTC, follicular variant thyroid papillary cancer; TCV-PTC, tall cell variant thyroid papillary cancer; ETE, extrathyroidal extension ETE. Variables with P values less than .1 in the univariate regression analysis were included in the multivariate regression equation. P values less than .05 are indicated in bold.

Figure 7.

Prognosis analysis of different T-I subgroups and the prognostic value of T-I index in patients with immune therapy. (A) Nomogram of clinical characteristics and T-I index group. (B) Kaplan–Meier survival analysis of T-I index groups. (C) The TIDE, dysfunction and exclusion score in different T-I index groups. (D) ROC curve of T-I index groups and single THBS1 expression groups to predict the responder of immune checkpoint suppression therapy according to the TIDE results. TIDE, the tumor immune dysfunction and exclusion score. ***, P < .001.

Molecular and Immune Characteristics of Different T-I Subgroups

To gain a complete insight into T-I index, we first analyzed gene mutations in different T-I subgroups using data from cBioPortal database. Missense mutation was the most common mutation type in both high and low T-I subgroups. The 15 genes with the highest mutation rates in all cases are listed in Figure 8A. BRAF was the most common mutated gene and showed the biggest difference in expression between high and low T-I subgroups (P < .001). While NRAS mutation occurred more in the low T-I group (P < .001).

Figure 8.

Molecular characteristics and TME landscape in different T-I subgroups. (A) Significantly mutated genes in the mutated TCGA-THCA cohort samples of different T-I subgroups. (B) Heatmap of tumor immune cell infiltrating score from 503 thyroid cancer samples in TCGA. Cases were divided into the HIGH and LOW groups according to the expression level of THBS1 and T-I index score. Cases are displayed on the X-axis; the Y-axis is clustered according to the immune cell subpopulation and infiltration fraction. (C) The proportions of 22 subpopulations of immune cells in different T-I index subgroups. The scattered dots represent the immune score of the 2 subgroups, and the score between them were compared through the Wilcoxon test. ns, not significant; *, P < .05; **, P < .01; ***, P < .001.

A heatmap of tumor immune cell infiltrating score from 503 thyroid cancer samples in TCGA cohort according to CIBERSORTx is shown in Figure 8B, with THBS1 expression level and T-I index grouping as patient annotations. Similar to the above results, we found that TFH and Tregs were more abundant in the high T-I subgroup (Figure 8C). Therefore, we speculated that the prognostic value of T-I index might result from both worse immune control and more aggressive cancer growth.

Discussion

THBS1 exhibits various roles in different tumors. Previous studies showed that THBS1 not only promotes tumor invasion and metastasis in PTC, but also inhibits angiogenesis in early PTC. THBS1 also regulates resistance to anaplastic thyroid carcinoma-targeted therapy. Therefore, the various roles of THBS1 in thyroid cancer deserve further elucidation.

In our research, analyses of TCGA cases showed that THBS1 expression level was higher in the thyroid tumor group compared with the normal group. Given the high positive correlation between THBS1 expression level and higher tumor pathological stage, lymph node metastasis, and ETE, we speculate that THBS1 plays an important role in tumor invasion and metastasis in PTCs. For the validation cohort from FUSCC, 53 paired frozen specimens were acquired from the tissue bank department. High THBS1 expression was positively correlated with a larger tumor size, higher TNM stage, more lymph node metastasis, and ETE. Although patients with higher THBS1 expression level showed a shorter PFI in TCGA cohort, THBS1 expression was not an independent risk factor in multivariate Cox regression analysis. Because of the short follow-up time of the FUSCC cohort, we analyzed the ability of THBS1 expression to predict high cancer recurrence risk. These results indicate that high expression of THBS1 in advanced thyroid cancer is associated with tumor invasion and metastasis, but its value as an independent prognostic biomarker is limited.

From the gene enrichment results, we identified a correlation between THBS1 expression and MAPK pathway. The activation of MAPK pathway is closely related to the occurrence and poor outcome of thyroid cancer. Previous studies have shown that inhibition of THBS1 expression can reduce the phosphorylation levels of ERK and MEK, thus inhibiting the invasion and metastasis of PTC. Therefore, we speculated that the high expression of THBS1 may promote cancer by promoting the activation of MAPK pathway. PD1 is a classic immune checkpoint molecule, and closely related to the role of T cells in tumor immunity. The enrichment of PD1 pathway in the high THBS1 expression group also suggested that the high expression of THBS1 may be closely related to tumor immunosuppression. Furthermore, KEGG and GO enrichment analysis results showed that the Wnt signaling pathway was enriched in the THBS1 high expression group. This signaling pathway is related to ETE, epithelial–mesenchymal transformation, tumor immunosuppression, and drug resistance in thyroid cancer.[37-39] In addition, MF enrichment analysis showed that high THBS1 expression was associated with tumor adhesion, which was consistent with previous findings in breast cancer. We also found that the integrin-binding pathway was enriched in the group with high THBS1 expression. Therefore, we performed immunohistochemical analysis of clinical specimens for co-expression of the neovascular marker CD31 and THBS1. However, we did not find clear evidence of THBS1 inhibiting neovascularization in thyroid cancer (data not shown). This result was not in line with findings from a previous study.

The prognosis of thyroid cancer, especially PTCs, is good and most of them can be cured by surgical resection. However, a few patients with advanced PTCs may be lost to surgery because of reasons such as huge neck masses, and immunotherapy may provide a new treatment option for these patients. Several studies have demonstrated higher levels of immune cell infiltration in PTC than in normal thyroid tissue, suggesting that PTC patients may benefit from immunotherapy.[40,41] Studies on how key genes in PTCs regulate its immune microenvironment have also been reported in the literature. For example, inhibition of BRAF sensitizes thyroid carcinoma to immunotherapy by enhancing tsMHCII-mediated immune recognition. In addition, in a trial of pediatric PTCs, investigators found that gene-fusion-driven PTCs were less differentiated and associated with more overrepresentation of mutations in tumor-immune crosstalk pathways. Therefore, we suggest that some genes that affect prognosis in PTCs may act simultaneously by regulating the immune microenvironment. It is essential to analyze such genes or molecular markers in conjunction with their relationship with tumor immunity.

Previous studies have examined the role of THBS1 in tumor immunity,[23,45] and our latest imaging histology study also found THBS1 is related to the tumor immunity-related pathways. However, the role for THBS1 in tumor immunity in thyroid cancer has not been studied. We first calculated the stromal score and immune score in different THBS1 expression groups using ESTIMATE. The higher ESTIMATE score in high THBS1 expression group indicated more immune cells in the TME. We then used ssGSEA and CIBERSORTx to analyze the tumor-infiltrating immune cells related to THBS1 expression in thyroid cancer. Both ssGSEA and CIBERSORTx results showed that the expression level of THBS1 was correlated with multiple immune cell subtypes, especially T cells. However, not all types of immune cells played a key role in the prognosis of thyroid cancer. Therefore, we further examined the association of immune cells with PFS using the TIMER website and clinical data from TCGA cohort. We found that TFH and Tregs infiltration were positively correlated with poor outcome in thyroid cancer. Given that high expression of THBS1 was positively correlated with the infiltration level of TFH and Tregs, we speculate that THBS1 may play a role in tumor immunosuppression via regulating these 2 types of immune cells.

Since the ability of THBS1 expression alone as a prognostic biomarker was limited, we next constructed a multigene prediction index with THBS1 and immune-related genes to improve the efficacy of THBS1 in prognosis prediction. After WGCNA, 4 THBS1-immune-related genes that significantly affected prognosis in thyroid cancer were selected to build the T-I index. Patients with a higher T-I index score showed a shorter PFI, and T-I index was proven to be an independent risk factor in multivariate Cox analysis. Furthermore, patients with a higher T-I index score showed a higher TIDE score including both Dysfunction score and Exclusion score, indicating that they are less likely to be the responder in ICI therapy. Therefore, the T-I index can also be used to predict the possibility of immunotherapy benefit.

The T-I index consists of 5 genes: IGHV3-49, CXCL13, TRAV8-3, TRBV30, and THBS1. C-X-C motif chemokine ligand 13 (CXCL13) is a well-known cancer-related gene. CXCL13 and its receptor CXCR5 (C-X-C motif chemokine receptor 5) are highly expressed in various tumors.[46,47] Overabundance of CD4( + ) CXCR5( + ) follicular helper T cells in thyroid cancer promoted tumor metastasis. This conclusion is consistent with the results that high THBS1 expression level positively correlated with TFH infiltration and led to a shorter PFI as described above. IGHV3-49 (immunoglobulin heavy variable 3-49) is a member of the IGHV family and mainly associated with B cell immunity and chronic lymphocytic leukemia. TRAV8-3 and TRBV30 belong to T cell receptor alpha and beta subgroups and directly affect T cell immunity. In a study of T cell receptors in patients with goiter cancer, the proportion of TRBV in tumor tissue, peripheral blood, and lymphocytes varied, supporting further study of immunity mechanisms against PTC. In the calculation formula of T-I index, the coefficient of these genes and THBS1 were positive numbers, indicating a positive correlation between T-I index and these genes. These findings indicate that T-I index is a biomarker associated with tumor promotion and tumor immunity suppression.

To obtain a complete insight of the T-I index, we then analyzed the gene mutation situation in different T-I subgroups with cBioPortal dataset. BRAF was the top mutated gene (77% mutation count) and most occurred in the high T-I group (P < .001). Notably, BRAF is the most common gene mutated in thyroid cancer, especially PTCs, and has been proven to be an important prognostic biomarker. The high BRAF mutation rate in the high T-I group may be due to the higher prevalence of C-PTC in this subgroup and be responsible for the poorer prognosis. In contrast, the higher mutation rate of NRAS in the low T-I group may have originated from the fact that FV-PTC was mainly concentrated in this group. This subtype is characterized by both papillary and follicular carcinoma molecular landscapes and therefore exhibits more NRAS mutations characteristic of follicular carcinoma. We also used the results of CIBERSORTx to analyze the immune characteristics in different T-I subgroups. Similar to result in the THBS1 expression level group, there were more TFH and Tregs infiltrations in the high T-I index group. As shown in the heatmap in Figure 6, some cases with low THBS1 expression were categorized into the high T-I subgroup. This is equivalent to a decrease in “false negative cases.” Thus, the T-I index showed effectiveness at predicting the benefit of immunotherapy.

The main shortcoming of this study is the short follow-up time in the validation cohort. In addition, we did not conduct gene sequencing on the validation cohort, so we could not obtain the absolute expression value of each gene or verify the infiltration of immune cells. In addition, to maintain consistency in the analysis of clinical information of patients in this study, we used the ETE and TNM grading criteria of the 7th edition of the AJCC guidelines, which is different from how the latest 8th edition of the guidelines grade ETE and tumors, which may lead to insufficient evidence of the correlation between THBS1 and ETE.

Conclusion

THBS1 is highly expressed in thyroid cancer, especially high-grade thyroid cancer. High THBS1 expression is associated with a shorter PFI and tumor immunosuppression. The T-I index is a valid prognostic biomarker for outcome and immunotherapy benefit in thyroid cancer.

Click here for additional data file.

Supplemental material, sj-docx-1-tct-10.1177_15330338221085360 for High Expression of THBS1 Leads to a Poor Prognosis in Papillary Thyroid Cancer and Suppresses the Anti-Tumor Immune Microenvironment by Anqi Jin, Jin Zhou, Pengcheng Yu, Shichong Zhou and Cai Chang in Technology in Cancer Research & Treatment