Integrative proteogenomic characterization reveals therapeutic targets in poorly differentiated and anaplastic thyroid cancers

Abstract

Poorly differentiated thyroid cancer (PDTC) and anaplastic thyroid cancer (ATC) present major challenges in treatment owing to extreme aggressiveness and high heterogeneity. In this study, deep-scale analyses spanning genomic, proteomic, and phosphoproteomic data are performed on 348 thyroid-cancer and 119 tumor-adjacent samples. TP53 (48%), TERT promoter (36.5%), and BRAF (23%) are most frequently mutated in PDTC and ATC. Ribosome biogenesis is identified as a common hallmark of ATC, and RRP9 silencing dramatically inhibits tumor growth. Proteomic clustering identified three ATC/PDTC subtypes. Pro-I subtype is characterized with aberrant insulin signaling and low immune cell infiltration, and Pro-II is featured with DNA repair signaling, while Pro-III harbors high frequency of TP53 and BRAF mutation and intensive C5AR1⁺ myeloid infiltration. Targeting C5AR1 synergistically improves antitumor effect of PD-1 blockade against ATC cell-derived tumors. These findings provide systematic insights into tumor biology and opportunities for drug discovery, accelerating precision therapy for virulent thyroid cancers.

Fig. 1. Multi-omics landscape of thyroid cancer.

A Workflow and sample types for molecular profiling of thyroid cancer. In this study, three types of thyroid cancer tissue samples (ATC, PTC, and PDTC) and two corresponding adjacent non-cancerous tissue samples (ATC NAT and PTC NAT) were collected. Whole-exome sequencing (WES) and TMT-labeled LC-MS/MS were respectively used for gene mutation and protein quantification profiling. 467 FFPE samples (124 ATC, 47 PDTC, 177 PTC (including 3 biological replicates), and 119 NAT samples) were detected in proteomic profiling. 217 FFPE samples (63 ATC, 25 PDTC, 81 PTC, and 48 NAT samples) were collected for phosphoproteomic profiling, and 97 FFPE samples for WES profiling (51 ATC, 28 PDTC, and 18 NAT samples). B An overview of detected samples in three molecular profiles (WES, proteomics, and phosphoproteomics). C The mutational landscape and associated clinical information of 51 ATC and 28 PDTC patients. Top 20 significantly mutant genes based on MutSig analysis (q value < 0.05). D Relationship between pathological features and gene mutations in 51 ATC and 28 PDTC patients. In boxplot, two-sided T-test is used, the central line represents median, the bounds of boxes represent the first and third quantiles and the upper and lower whiskers extend to the highest or the smallest values with 1.5 interquartile range. In contingency tables, the value represents the number of samples. Statistical significance was examined by Fisher’s exact test. E Kaplan–Meier curves showing the relationship between gene mutation and patient survival (n = 51 samples). Source data are provided as a Source Data file.

Fig. 2. Proteomics and phosphoproteomics analyses in tumor tissues compared with NATs.

A Volcano plot showing the differential protein abundance between 113 ATC and 20 NAT samples (FDR < 0.05, T/N > 1.5 or <0.67). Statistical significance was examined by two-sided Wilcoxon rank-sum test. P values were adjusted by Benjamini–Hochberg method. B Heatmap showing the abundance of significantly different proteins between 113 ATC and 20 NAT samples. Pathways significantly enriched with upregulated and downregulated proteins are labeled on the right side of the heatmap. C Scatter plot showing the common and specific protein changes in ATC (n = 113) and PTC (n = 177) compared to NAT (n = 119) samples. Statistical significance was examined by two-sided Wilcoxon rank-sum test. D Kaplan–Meier curves (Log-rank test) showing the relationship between the protein abundance of ARF4 and ARHGAP9 and patient overall survival (OS) (n = 101 samples). E Scatter plot showing the relationship between protein abundance changes and phosphorylation abundance changes in ATC compared to NAT samples. F Bubble plot by Hypergeometric test showing significantly enriched pathways for proteins with fold changes greater than 1.5 in both phosphorylation and protein levels. G Bar plot showing significantly changed kinase activity in ATC compared to NAT samples. Kinase activity is estimated depending on the substrate phosphorylation abundance and calculated using the KSEA method. H The changes of potential substrate phosphorylation abundance in patients with BRAF^V600E or non-BRAF^V600E mutation. I Network diagram showing the relationship between potential substrates influenced by BRAF^V600E mutation and the pathways they participate in. J Kaplan–Meier curves showing the OS of patient with ATC based on BRAF^V600E mutation and phosphorylation abundance of ERBB2, respectively. K Multivariate analyses (Cox proportional hazards models) of gene mutation and phosphorylated ERBB2 in patients with ATC. Whiskers represent 95% confidence interval. In results of (E–G), 113 ATCs and 20 NATs in proteomics and 63 ATCs and 20 NATs in phosphoproteomics were analyzed. In results of (H–K), n = 51 samples. *P < 0.05, **P < 0.01. Source data are provided as a Source Data file.

Fig. 3. Molecular signature of different thyroid cancer subtypes.

A Heatmap showing the abundance changes of proteins significantly associated with thyroid cancer progression. Pathways significantly enriched with proteins specifically highly expressed at different stages of the disease are labeled on the right side of the heatmap. B Boxplot showing proteins gradually increasing or decreasing in abundance as thyroid cancer progresses. The central line represents median, the bounds of boxes represent the first and third quantiles and the upper and lower whiskers extend to the highest or the smallest values with 1.5 interquartile range. N indicates number of increasing or decreasing proteins. C Displaying protein abundance changes in pathways significantly enriched with downregulated proteins as cancer progresses from normal thyroid towards PTC (differentiated phenotype), PDTC (poorly differentiated phenotype), and ATC (dedifferentiated phenotype). D Displaying protein abundance changes in pathways significantly enriched with upregulated proteins as thyroid cancer progresses. E Kaplan–Meier curves (Log-rank test) showing the relationship between the estimated activity changes of the ribosome biogenesis pathway based on ssGSEA and patient overall survival (n = 101 samples). F Cell viability after knockdown of RRP9 in 8505C (n = 3) and CAL62 (n = 3), respectively. The experiments were repeated for three times. G Colony formation assay after silence of RRP9 in 8505C and CAL62, respectively. Data are a representative of three separate experiments. H–J The tumor growth curve, tumor size, and tumor weight of 8505C-derived ATC xenograft in nude mice (n = 8). In results of (A–D), 113 ATCs, 177 PTCs, 46 PDTCs, and 119 NATs were analyzed. Results (F, H, J) show the mean of biological replicates. Error bars indicate SD for in vitro experiments and SEM for in vivo experiments. Statistical significance was determined by two-sided T-test. Source data are provided as a Source Data file.

Fig. 4. Proteomic Clustering and molecular characteristics in ATC/PDTC tumors.

A Three subtypes (Pro-I, Pro-II, and Pro-II) of ATC and PDTC patients were divided based on protein expression profiles. The bottom heatmap shows the abundance of featured proteins in each subtype. The top heatmap shows the pathological characteristics of patients in three subtypes. B Kaplan–Meier curves (Log-rank test) showing significant differences in patient overall survival among different protein subtypes (Pro-I: 38 samples; Pro-II: 18 samples; Pro-III: 45 samples). C Bubble plot by Hypergeometric test showing pathways significantly enriched with differentially expressed proteins of three subtypes. D Heatmap showed the abundance of characteristic proteins under key regulatory pathways in three protein subtypes. E Radar plot showed differences in kinase family activity among three protein subtypes of patients. F Bar plot showed the activity of key regulatory kinases in three protein subtypes. G In the three protein subtypes, downstream regulatory substrates of significantly activated kinases involved in different signaling pathways. Source data are provided as a Source Data file.

Fig. 5. Immune subtyping and regulatory factors in ATC/PDTC.

A The immune subtypes of ATC and PDTC patients were divided into three based on immune cell scores calculated using xCell. The bottom heatmap shows the abundance of immune cell markers and the activity of pathways (estimated by the GSEA method based on protein abundance) in patients with three immune subtypes. B Heatmap of immune score and stromal score in three immune subtypes. C Kaplan–Meier curves (Log-rank test) showing differences in patient overall survival among three immune subtypes (Cluster1: n = 46 samples, Cluster2: n = 9 samples, Cluster3: n = 46 samples). D Sanky plot showing the composition of ATC and PDTC samples in immune subtypes and proteomic subtype. E Boxplot of the immune score in three immune subtypes (Cluster1: n = 74 samples, Cluster2: n = 11 samples, Cluster3: n = 74 samples). F The survival risk scores (Cox proportional hazards models) of 343 transcription factors (TFs) collected from the ENCODE Project Consortium database based on protein abundance. G Bubble plot of transcription factors in terms of patient survival risk scores and differentially changed in Cluster 3 compared to Cluster 1. Two-sided Kruskal–Wallis test is used to estimate statistical significance. H Boxplot showed the differences in protein abundance of transcription factor STAT1 in three immune subtypes (Cluster1: n = 74 samples, Cluster2: n = 11 samples, Cluster3: n = 74 samples). In boxplot (E, H), two-sided Kruskal–Wallis test is used, the central line represents median, the bounds of boxes represent the first and third quantiles and the upper and lower whiskers extend to the highest or the smallest values with 1.5 interquartile range. Source data are provided as a Source Data file.

Fig. 6. Infiltration of neutrophil resulted in poor prognosis in ATC.

A Bar plot showing significantly altered (q value < 0.1) chromosomal arms based on copy number alterations (CNA) in ATC and PDTC samples. B Heatmap showing the correlation between copy number alterations and protein abundance in ATC and PDTC samples. C Scatter plot showing the Spearman correlation between CNA of genes in the 1q chromosomal arm region and their protein abundance. D Venn diagram showing the overlap between cis-effect gene in the 1q region and significantly upregulated protein in ATC. E Expression and survival risk scores of top three proteins from the overlapping section. Heatmap showing protein expression in 20 NT and 113 ATC samples. Whiskers represent 95% confidence interval. F Immunofluorescence staining showing representative micrographs of MPO and citrullinated histone H3 (cit-H3) in normal tissues (NT), PTC, and ATC. n = 3 samples for each group. G Scatter plot showing the Spearman correlation (two-sided) between the protein abundance of NCF2 and neutrophil score calculated based on xCell (n = 159 samples). H, I Spearman correlation (two-sided) of NCF2, MPO, and PAD4 abundance in ATC/PDTC samples (n = 159). J, K Immunofluorescence staining of C5AR1 in NT (n = 30), PTC (n = 35), and ATC (n = 19). ATCs comprised of 9 ATCs that included in the proteomic analysis and belonged to 3 proteomic subtypes (3 samples for each subtype), and 10 ATCs that were not included in the proteomic analysis. Data are presented as mean ± SD, and statistical analysis is performed via one-way ANOVA. L The relevance of C5AR1 and overall survival (Log-rank test) in patients with ATC or PDTC. M, N Spearman correlation (two-sided) of C5AR1, MPO, and PAD4 abundance in ATC/PDTC samples (n = 159). O Heatmap illustrates the expression and activity of STAT1 and regulatory effect on its downstream substrate proteins correlated with neutrophil infiltration. P Flowchart illustrating the potential biological mechanisms associated with poor prognosis in ATC patients. Source data are provided as a Source Data file.

Fig. 7. Potential targets for three molecular subtypes of ATC/PDTC.

A Proposal of potential targets depending on proteome-based ATC/PDTC subtyping. B, C Immunohistochemistry staining of BCL2 in three proteomic subtypes (n = 4 samples in each group). Data are presented as mean ± SD, and statistical analysis was performed via one-way ANOVA, two-sided. D, E Immunofluorescence staining of C5AR1 and myeloid marker CD11b in three proteomic subtypes (n = 3 samples in each group). Data are presented as mean ± SD, and statistical analysis was performed via one-way ANOVA, two-sided. F, G The tumor growth curve and tumor size in an immunocompetent model of murine ATC after treating with avacopan, and/or anti-PD-1 as single agent and combination therapy (n = 6). Results show the mean of biological replicates. Error bars indicate SEM. H Immunofluorescence staining showing representative micrographs of MPO and cit-H3 in mATC-derived tumors after avacopan treatment. Experiment was repeated independently using 3 samples for each group with similar results. Source data are provided as a Source Data file.