Molecular and Proteomic Profiles of Radioiodine Refractory Papillary Thyroid Cancer

2.1 Patient cohort

Based on the inclusion and exclusion criteria in the Patients and Methods section, forty-eight patients with locoregional or distant metastatic PTCs were recruited for our study (Table 1). The average age of all patients was 38.4 ± 12.69 years, and the patients in the RR-PTC group were slightly older than those in the control group (P = 0.421). Interestingly, all five patients ≥ 55 years old were in the RR-PTC group. Although females were more common in the control group (72.7% vs. 54.1%, P = 0.319), the difference was not significant. In addition, no differences were found in the tumor diagnosis, smoking status or alcohol consumption, body mass index or the incidence of concomitant thyroid diseases between the two groups.

The histopathological subtypes of PTC were identically distributed between the two groups. No significant differences were found in terms of tumor size, multifocality or bilateral lesions. The incidence of extrathyroidal extension (ETE) was greater in the RR-PTC group with marginal statistical significance (P = 0.095). Owing to differences in age and ETE, T2-4 categories were more common in patients with RR-PTC (91.9% vs. 63.6%, P = 0.039). The presence of lymph node metastatic lesion (LNM) was similar in the two groups, but distant metastasis was much more common in the patients with RR-PTC (56.8% vs. 9.1%, P = 0.006). Consequently, higher AJCC stages were observed in RR-PTC patients (P = 0.001). Compared with patients in the control group, patients in the RR-PTC group received larger doses of I-131 after the initial resection of primary tumors. Notably, all six patients who died were in the refractory group.

2.2 Molecular profile of PTCs correlated with I-131 resistance

A total of 23 genes were sequenced (Sup table S1). Four genes were found to harbor driver mutations (Table 2). The BRAF^V600E mutation was identified in 25 samples (56%), with no significant difference between the RR-PTC group and the control group (54.3% vs. 60.0%, P = 1.000). All six samples with TERTp mutations (17.1% vs. 0%, P = 0.312) were in the RR-PTC group. In addition, RET and TP53 mutations were found in one and two patients, respectively. RAS mutations, including HRAS, KRAS and NRAS, were not detected in any samples. In addition, four kinds of gene fusions were found, namely, eight samples with NCOA4-RET, two with CCDC6-RET, three with ETV6-NTRK3 and one with STRN-ALK. Fusions were more frequent in the control group (28.6% vs. 40.0%, P = 0.700), whereas RET-related fusions were almost equally distributed (22.9% vs. 20.0%, P = 1.000). However, no significant differences in the expression of these genes were found between the two groups.

The forty-five samples were then divided into four subgroups. The BRAF^V600E mutation was less common in the continuous RAI uptake (ID) group (11.1%, P < 0.05), whereas the TERTp mutation was more frequent in the loss of uptake at the first I-131 treatment (iDF) group (33.3%, P < 0.05). Notably, NCOA4-RET fusion was more common in the ID group (44.4%, P < 0.05).

The present study further evaluated the associations between gene variants and I-131 uptake. The samples in the ID and RAI-sensitive PTC (Id) groups were categorized as positive for I-131 uptake. BRAF^V600E mutation was more frequent in the negative I-131 uptake group with statistical significance (69.2% vs. 37.8%, P = 0.039). The incidence of TERTp mutations was also greater in this group (19.2% vs. 5.3%, P = 0.222). In contrast, fusions were more common in the positive uptake group (23.1% vs. 42.1%, P = 0.206).

The distribution of gene variants with unknown clinical significance was also assessed (Sup table S2). Four variants were subjected to nonparametric tests after those with a frequency < 3 were filtered out. The four mutations were located on ATM (c.1236-2A>T), PTEN (c.402G>A) and RET (c.2071G>A and c.2671T>G). Three variants were more common in the RR-PTC group, except for PTEN c.402G>A, although the difference was not statistically significant. In subgroup analysis, the incidence of PTEN c.402G>A was greater in the ID group (22.2%, P < 0.05). In addition, PTEN c.402G>A was associated with positive I-131 uptake with marginal significance (0% vs. 15.8%, P = 0.068).

2.3 Overview of the proteomic profiling

To profile the proteomic characteristics of RR-PTC, we collected 168 tissue slides from 73 patients. A total of 9769 proteins were identified (Sup figure 1A), which indicated the high quality of our data. Since the number of proteins identified in all the samples exceeded the threefold interquartile range (≤ 1614), no sample was excluded from subsequent analysis. The principal component analysis (PCA) plot annotated by batch groups indicates no observable batch effect during data acquisition (Sup figure 1B). Although samples were well resolved by their site, tumor tissues mixed with little difference but negative lymph nodes were clustered (Sup figure 1C). The R values of correlation analyses exceeded 0.97 in both intra-batch and inter-batch tests, indicating high technical stability across the data acquisition (Sup figure 1D&E). Overall, no observable batch effect was detected (Sup figure 1F).

Next, the raw files were converted to an expression matrix. We applied a filtering criterion to select proteins with missing values less than 50%. These 7406 filtered proteins were then subjected to data normalization and outlier substitution (Sup figure 2A-C). After ID conversion, an expression profile encompassing 7394 proteins across 168 samples was compiled (Sup figure 2D). As indicated by the ESTIMATE algorithm, the median tumor purity was 76.2% (Sup figure 2E), indicating the reliability of our dataset for further analysis.

The global proteomic profiles of the original PTC lesions were subjected to dimensionality reduction to evaluate the clustering of the samples. The area covered by the samples represented similarity within a group. The RAI-sensitive samples were more tightly clustered than the RR-PTC samples (Figure 1A). The overlap of the two groups might be explained by the similarity between tumor tissues and the intrinsic heterogeneity of RR-PTC samples. To minimize potential heterogeneity, we furtherly divided RR-PTCs into three subgroups (iDF, iDG and ID) as described in the patients and tissue selection section. The diversity of the ID group was the lowest among the three subgroups, while the diversity of the RAI uptake lost gradually (iDG) group was the highest (Figure 1B). Interestingly, when labeled with the BRAF^V600E-mutation status, samples in the two groups were separately clustered (Figure 1C).

• Download figure
• Open in new tab

Figure 1.

Proteomic overview of RAI-refractory and RAI-sensitive PTC samples. (A-C) Dimension reduction visualization of the analyzed PTC samples. The samples were labeled according to their response to I-131 treatment (refractory or sensitive), subgroup (iDF, iDG, ID or Id) and BRAF mutation status. Three samples whose quality control for sequencing failed were not included in Figure C. (D) Differential expression analysis of samples in the RR-PTC and control groups. Unadjusted P values < 0.05 and absolute FC > 1.414 were defined as thresholds for DEPs. (E) Overview of the proteomic profiles. Z scores of the DEPs were calculated. (F) Correlation of scores and groups. Error bar denoted 95% CI. The TDS was lower in the RR-PTC group (−0.17 ± 0.562 vs. 0.03 ± 0.512, P = 0.288). Conversely, the ERK score was greater in the RR-PTC group (0.03 ± 0.317 vs. −0.10 ± 0.297, P = 0.215). The means and SDs are presented with boxes and error bars.

Abbreviations: age at D, age at diagnosis; DM, distant metastasis; ETE, extrathyroidal extension; LNM, lymph node metastasis.

Subsequently, 107 differentially expressed proteins (DEPs) with P values < 0.05 and absolute fold change (FC) > 1.414 were identified via the multiple linear regression algorithm (Figure 1D). Eighty-nine proteins were upregulated in the RR-PTCs, whereas the remaining eighteen proteins were downregulated. The top upregulated protein in the RR-PTCs was KRT71 (keratin, type II cytoskeletal 71, log₂FC = 1.422). Interestingly, the top downregulated protein also belonged to the keratin family (keratin, type I cytoskeletal 16, KRT16, log₂FC = −0.955), which suggested the potential involvement of the cytoskeleton in the progression of RAIR. Samples were clustered unsupervisely based on the 107 DEPs and then labeled with clinicopathological features (Figure 1E). The borders of the clinicopathological labels were generally indistinct. However, the thyroid differentiation score (TDS) and extracellular signal-regulated kinase (ERK) score were correlated with the groups to some extent (Figure 1F). The TDS was lower in the RR-PTC group, which could indicate poorer differentiation (P = 0.288). Conversely, the ERK score was greater in the RR-PTC group (P = 0.215). Unfortunately, none of these scores reached statistical significance.

Gene set enrichment analysis (GSEA) based on the Reactome database revealed 64 pathways with statistical significance (Sup figure 3A). Twenty-nine pathways with positive normalized enrichment scores (ESs) were unregulated in RR-PTC samples. The pathway with the greatest increase was gamma carboxylation hypusine formation and arylsulfatase activation (nES = 1.687). In addition, Akt phosphorylates targets in the cytosol was also upregulated (nES = 1.552). The remaining 39 pathways were downregulated. Regulation of mRNA stability by proteins that bind Au-rich elements was mostly suppressed for RR-PTCs (nES = −1.960). Notably, RAS processing was also downregulated (nES = - 1.728). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was also used for the annotation of pathways closely related to oncogenesis and dedifferentiation (Sup figure 3B). P53 signaling was significantly enriched (nES = 1.5440, P = 0.0136). Interestingly, the MAPK signaling pathway showed a slight positive correlation with RR-PTCs (nES = 1.2623, P = 0.1642). The gene set variation analysis (GSVA) results supported the above results (Sup figure 3C).

A previously published study assessed the transcript profiles of RR-PTC via a gene microarray ²⁸. The DEPs identified in the present study were assessed via this profile to determine whether they were robust in another cohort. Ninety-nine of the 107 DEPs were identified in this cohort, and only three of them were dysregulated (Sup figure 3D). BET1 and C11orf96 were downregulated in the RR-PTC group (BET1 median: 6.534 vs. 6.105, P = 0.032; C11orf96 median: 7.338 vs. 7.028, P = 0.032). The expression of SLC4A9 was greater in the RR-PTC group (median: 5.216 vs. 5.392, P = 0.045). Notably, the regulation direction of BET1 and SLC4A9 was consistent with the results of the present study, whereas C11orf96 was oppositely expressed in the two analyses. GSVA based on the Reactome and KEGG databases revealed two differentially downregulated pathways (N-glycan antennae elongation and GAB1 signalosome), none of which were enriched in our study (Sup figure 3E&F).

To assess the spectral change in RR-PTC, the Id, ID and iDF subgroups were defined as having low, intermediate and high severity, respectively. A total of 403 proteins with significance difference in analysis of variance were divided into six clusters (Sup figure 4). In two clusters, 160 proteins were monotonically regulated (Sup figure 5A). Collagen biosynthesis and modifying enzymes were significantly enriched in Metascape (Sup figure 5B). In addition, regulation of protein autophosphorylation and regulation of extrinsic apoptotic signaling pathway were also associated with the severity of RR-PTC. The results were supported by the ingenuine pathway analysis (IPA) database (Sup figure 5C).

2.4 Proteomic biomarkers and pathways associated with I-131 uptake

For analysis of proteomic biomarkers and the underlying mechanisms for the loss of I-131 uptake in PTC, two subgroups, ID and iDF, were included in the subsequent analysis of this section. The ability of these two groups to take up I-131 was completely different at the first radiotherapy. The demographic and clinicopathologic characteristics of the two groups were roughly comparable (Sup table S3). No significant differences were observed in histopathological features, tumor staging or risk stratification. Interestingly, four patients in the iDF group drank alcohol and/or smoked, while no individual in the ID group reported these behaviors (P = 0.033). In addition, the patients in the iDF group had greater weights than their counterparts (P = 0.010). Since the diagnosis of RR-PTC in the ID group required multiple rounds of I-131 treatment, the cumulative dose of I-131 was inevitably greater than that in the iDF group (P = 0.001).

Although the ID and iDF groups overlapped with each other to some extent, a total of 145 DEPs were identified with P values < 0.05 and absolute FC > 1.414 (Figure 2A&B). Ninety-five proteins were upregulated in the ID group, whereas the remaining 50 proteins were highly expressed in the iDF group. The 145 DEPs could barely separate the two groups. The TDS was greater in the ID group, whereas the ERK score was lower. However, none of these scores reached statistical significance (Figure 2C&D).

• Download figure
• Open in new tab

Figure 2.

Differential expression analysis of the ID and iDF groups. (A) PCA plot of the two groups. (B) Differential expression analysis of samples in the two groups. Unadjusted P values < 0.05 and absolute FC > 1.414 were defined as thresholds for DEPs. (C) Correlations of scores and groups. Error bar denoted 95% CI. The TDS was greater in the ID group (0.03 ± 0.599 vs. - 0.20 ± 0.454, P = 0.356), while the ERK score was greater in the iDF group (−0.07 ± 0.345 vs. 0.07 ± 0.300, P = 0.386). (D) Expression profiles of 145 DEPs. (E&F) GSEA based on the KEGG and Reactome databases. Sixteen pathways are presented in the plot.

Pathway enrichment analysis was then conducted. Only seven pathways in the KEGG database were significantly enriched. Among the seven pathways, glutathione metabolism and Notch signaling were enriched in the ID group whereas TGF-beta signaling and ECM receptor interaction were enriched in the iDF group (Figure 2E). Fifty-one pathways in the Reactome database were significantly enriched. The top enriched pathway was “activated TAK1 mediates P38 MAPK activation” with an nES of 1.9182 (Figure 2F). In addition, several MAPK-related pathways, including ERK/MAPK targets and prolonged ERK activation events, were correlated with positive I-131 uptake. In addition, TNF signaling and glutathione conjugation were associated with positive I-131 uptake.

The results of pathway enrichment were then assessed via GSVA. Sixty-three pathways were significantly enriched and 18 of these pathways overlapped with the GSEA results (Figure 3A). After the associations between pathways and tumor progression were evaluated, sixteen pathways were found to play roles in the loss of I-131 uptake in PTC, eleven of which overlapped with the GSEA results (Figure 3B). Six pathways were closely related to MAPK signaling. Three pathways are associated with transforming growth factor-β activated kinase 1 (TAK1) signaling. In addition, one pathway was correlated with both MAPK and TAK1 signaling. The other six pathways involved glutathione metabolism, TNF signaling and MET proto-oncogene (MET) signaling. The expression of proteins involved in these thirteen pathways was compared between the RAI-positive and RAI-negative groups. Proteins with an absolute FC > 1.2 were displayed (Figure 3C).

• Download figure
• Open in new tab

Figure 3.

Expression of selected genes according to the enriched pathways identified in the ID vs. iDF and ID vs. Id comparisons. (A&D) Dysregulated pathways identified in the ID vs. iDF and ID vs. Id comparisons. The P value threshold was 0.05. No threshold was set for fold change. (B) The ten most significantly enriched pathways and the thirteen selected pathways. (C) The selected proteins from the thirteen selected pathways. The threshold of absolute fold change was > 1.2. (E) The eight significantly enriched pathways and thirteen selected pathways. (F) The selected proteins from the thirteen selected pathways. The TNF signaling pathway and the TAK1-related pathways were combined in a single plot since the proteins strongly overlapped.

* P < 0.05, ** P < 0.01.

^a The pathway “MET activates PTPN1” was also selected for subsequent identification of proteins.

^{b, c} Involved in the selected pathways “glutathione synthesis and recycling” and “MET activates PTPN11”, respectively.

^d Both CDK1 and MAP2K1 participate in the pathways “MAPK3 (ERK1) activation” and “RAF-independent MAPK1/3 activation”, and MAP2K1 is also a participant in “prolonged ERK activation events”.

^e TAB2 and TAB1 are involved in both MAPK and TAK1-associated pathways, while CD14 participates in TAK1-associated pathways.

^f Molecules in these figures might participate in several pathways. The associations of molecules and pathways are provided in Sup table S5.

The results of our proteomic profile were validated with a transcriptomic profile of RR-PTCs (GSE151179). The samples in GSE151179 were categorized as ID or iD based on their I-131 uptake. The iD group was not subdivided into iDF or iDG since the information was not provided. A linear regression algorithm identified three common differentially expressed genes (DEGs) (Sup figure 6A). However, the regulatory direction was completely different. In the proteomic profile, MYH7 was up-regulated, whereas in the transcriptomic profile, its expression was greater in the iD group. S100B and UXT also exhibited inverse changes in the two profiles. Only MYH7 significantly differed according to the nonparametric test (P = 0.035). GSVA revealed five common pathways (Sup figure 6B), and the direction of the change in expression was consistent between the two profiles. Notably, three TAK1-associated pathways were significantly enriched. All genes involved in the 16 selected pathways and these five common pathways were tested. Only three genes were differentially expressed between the two groups. PI4K2B and PIK3R4 participate in the synthesis of PIPs at the early endosome membrane, and CEACAM8 is involved in fibronectin matrix formation (Sup figure 6C).

Immune analysis based on the ssGSEA matrix revealed an abundancy of memory B cell and a deficiency of macrophage and immature T cell in the ID samples (Sup figure 7A). Notably, eosinophils were more abundant in RAI-sensitive samples compared to RR-PTCs. However, immunohistochemistry results did not show any statistically significant differences (Sup figure 7B).

2.5 Proteomic biomarkers and pathways for different responses to positive I-131 uptake

The data of twenty-one patients in the ID group and the Id group (or control group) were compared to determine the potential mechanisms underlying the different responses to positive I-131 uptake in PTC patients. The baseline characteristics of the two groups were compared (Sup table S4). No significant differences in demographic or tumor features were found, except for a greater incidence of ETE in the ID group (P = 0.035). More distant metastases were found in the ID group (80% vs. 9.1%, P = 0.002), and consequently, the disease stages were greater (P < 0.001). According to the criteria of the ID group, the cumulative I-131 dose was inevitably greater (P < 0.001).

Similar to the previous comparison, the two groups could be distinctly separated from each other (Figure 4A). Twenty-one samples were roughly divided into two groups with 125 DEPs identified via the multiple linear regression algorithm (Figure 4B&D). No significant differences were found in the TDS, ERK score, T cell infiltration score (TIS), immune infiltration score (IIS) or immune cytolytic activity score (CYT) (Figure 4C).

• Download figure
• Open in new tab

Figure 4.

Differential expression analysis of the ID and Id groups. (A) PCA plot of the two groups. (B) Differential expression analysis of samples in the two groups. Unadjusted P values < 0.05 and absolute FC > 1.414 were defined as thresholds for DEPs. (C) Correlations of scores and groups. Error bar denoted 95% CI. No significant difference was found in the TDS (0.03 ± 0.599 vs. 0.03 ± 0.512, P = 0.985) or ERK score (−0.07 ± 0.345 vs. −0.10 ± 0.297, P = 0.796). (D) Expression profiles of 125 DEPs. (E&F) GSEA based on the KEGG and Reactome databases. Sixteen pathways are presented in the plot.

GSEA revealed twelve enriched pathways based on KEGG annotations. The p53 signaling pathway (P = 0.0498) and mTOR signaling pathway (P = 0.0346) were both promoted in the RAI-resistant samples (Figure 4E). Interestingly, inositol phosphate (IP) metabolism was also activated in the ID group (P = 0.0190). No significant difference was revealed in the MAPK signaling pathway. Forty-five pathways were significantly enriched according to the Reactome annotation file. TP53-related pathways, including regulation of TP53 through association with co-factors (P = 0.0260) and methylation (P = 0.0463), were increased in ID samples (Figure 4F). TNF signaling (P = 0.0188) and its participant TNFR1-induced NF-kappa-B signaling pathway (P = 0.0021) were both associated with RAI uptake and resistance. In addition, the pathway “activated TAK1 mediates p38 MAPK activation” was closely related to RAI-resistance (P < 0.0001). An increase in IP metabolism was also detected (P = 0.0081).

GSVA based on the Reactome annotation file was conducted to verify the previous results. A total of 49 pathways were enriched with significant differences, eleven of which overlapped with previous results (Figure 3D). After the search and selection, four clusters comprising thirteen pathways were pooled for subsequent analysis (Figure 3E). Eight of the thirteen pathways were also dysregulated according to GSEA. A total of 39 proteins involved in TNF signaling, TAKI-related pathways, plasma lipoprotein remodeling and IP-related pathways were identified with an absolute FC > 1.2 (Figure 3F).

The expression levels in the transcriptomic profiles of the ID and Id groups were then compared for validation. Only PGS1 was downregulated in the ID group (median: 6.080 vs. 6.456, P = 0.012; Sup figure 6D), which was consistent with the proteomic data. According to the pathway analysis, only GAB1 signalosome was significantly decreased in the ID group (P = 0.012, Sup figure 6E). However, none of the genes in the GAB1 signalosome were differentially expressed between the two groups. Among those involved in the above thirteen pathways, three genes (MAPKAPK2, NFKB1 and PLKHA6) were downregulated in the ID group, in contrast to the expression pattern observed in our proteomic profile (Sup figure 6F).

2.6 Association of gene variance and protein expression

The 7394 proteins were divided into 21 modules via weighted correlation network analysis (WGCNA) (Figure 5A&B). BRAF mutation was significantly correlated with seven gene modules, among which the blue module had the highest absolute correlation coefficient (CC) (CC = 0.61, P = 8.4e-6). The most prominent module associated with TERTp mutations was pink (CC = −0.45, P = 1.9e-3). The green-yellow module was closely related to the gene fusions (CC = 0.42, P = 3.7e-3) and the NCOA4-RET fusion (CC = 0.32, P = 0.03). For other gene variants of unknown significance, ATM c.1236-2A>T was related to the green module (CC = 0.41, P = 5.0e-3). In addition, RET c.2671T>G was closely associated with the yellow module (CC = 0.50, P = 4.5e-4).

• Download figure
• Open in new tab

Figure 5.

Gene variants and weighted gene co-expression network analysis. (A) Gene variance landscape. (B) Correlations of gene modules and gene phenotypes. The cells were divided into upper left (correlation coefficient) and lower right (P value) parts. (C) The intersection of LASSO-, RF- and SVM-RFE-identified biomarkers for BRAF, and TERTp mutations, gene fusions and NCOA4-RET fusion. (D) Validation of six genes in the TCGA-THCA cohort and merged GEO cohort. PTCD1 was not detected in the merged GEO dataset.

^a The low number of TERTp-positive samples (n = 6) prevented model construction and error analysis for TERTp mutations via the SVM algorithm. The first ten proteins with the top average rank were processed for biomarker identification.

Three machine learning algorithms (LASSO, RF and SVM-RFE) were used to select biomarker proteins from the most correlated modules of each gene mutation and fusion. Since other gene mutations were generally rare (< 10%) in the present study, samples with BRAF mutation, TERTp mutation, gene fusion and NCOA4-RET fusion were retained for subsequent analysis. A total of 28, 32 and 39 proteins were used for model construction of BRAF mutations via the least absolute shrinkage and selection operator (LASSO) regression, random forest (RF), and support vector machine-recursive feature elimination (SVM-RFE) algorithms, respectively, and four of these proteins were common biomarkers (Figure 5C). In addition, only PTCD1 and AHCYL2 were revealed to be common biomarkers for gene fusions. According to the receiver operating characteristic (ROC) curves of The Cancer Genome Atlas Thyroid Cancer (TCGA-THCA) dataset, the areas under the curve (AUCs) of three genes exceeded 0.9 (Figure 5D). PTPRE had the highest accuracy (AUC = 0.93), followed by TMEM43 (AUC = 0.92) and FN1 (AUC = 0.91). According to the merged Gene Expression Omnibus (GEO) profile, the most accurate gene was TMEM43 (AUC = 0.80), followed by FN1 (AUC = 0.73) and SUSD2 (AUC = 0.69). Interestingly, PTPRE did not significantly differ. The AUCs of AHCYL2 for gene fusions in TCGA-THCA and merged GEO cohorts were 0.65 and 0.68, respectively. Disease free survival analysis revealed that dysregulation of FN1 (HR = 2, P = 0.026) and AHCYL2 (HR = 0.53, P = 0.038) was associated with poor prognosis in patients with thyroid cancer (Sup figure 8).

To explore the potential biomarkers and pathways in case of BRAF mutation and wild-type TERTp, only 19 BRAF-mutated and TERTp-unmutated samples were analyzed. A total of 184 proteins were significantly different between the RAIR-refractory and RAIR-sensitive groups (Sup figure 9A&B).

Metascape revealed that several pathways related to the immune response, including T cell activation, regulation of lymphocyte activation and phagocytosis, were enriched with statistical significance (Sup figure 9C). IPA results also supported these findings (Sup figure 9D). Granzyme A signaling was downregulated in BRAF-wild-type samples, whereas neutrophil degranulation, B cell development and phagosome formation were upregulated in BRAF-mutated TERTp-unmutated samples.

4.1 Study design and ethics approval

This retrospective, case-control study was conducted at Renmin Hospital of Wuhan University. The Institutional Ethical Committee of the hospital reviewed and approved the study design (No. WDRY2021-K032). The requirement for obtaining informed consent from the involved patients was waived due to the retrospective nature of the study. The study was conducted in accordance with the Declaration of Helsinki ⁶². The study was performed in accordance with the STROBE checklist for case-control study (version 4).

4.2 Patient and tissue selection

The medical records of patients who were diagnosed with thyroid cancer at our tertiary center from Jan. 1^st, 2016 to Dec. 30^th, 2022 were reviewed (Figure 6A). Patient demographic and clinicopathological characteristics, laboratory test results and radiological results were collected by two researchers and independently evaluated. During the initial screening, patients who a) had incomplete medical records and b) did not undergo I-131 WBS were excluded. In the second round of assessment, histopathological data and radiological images, including chest computed tomography (CT), cervical magnetic resonance imaging (MRI), whole-body bone scanning and positron emission tomography (PET), were used to evaluate local and distant metastasis at the first admission for surgery. Patients without cervical or distant metastasis at diagnosis were excluded from the subsequent evaluation. Patients with RAI uptake limited to thyroid bed and disease remission were also excluded. In the final assessment, RAI uptake was determined by I-131 WBSs or I-131 single photon emission computerized tomography (SPECT). A total of 73 patients were identified to have locoregional or distant metastasis, which was confirmed by radiological, cytopathological or histopathological solid evidence. Patients with primary lesions not resected in our center or other histopathological types other than PTC were also excluded from further statistical analysis. Finally, forty-eight PTC patients with locoregional and/or distant metastasis were included in this study and their primary tumor samples (one sample from each patient) were used for subsequent analysis.

• Download figure
• Open in new tab

Figure 6.

(A) Flow chart of patient and sample selection. Strict criteria were followed to identify patients with RAI-refractory or RAI-sensitive thyroid cancers. Individuals with insufficient radiological evidence were also excluded from the second round of selection. A total of 73 patients (in the dashed box) were revealed to have regional or distant metastatic lesions with solid evidence. The FFPE tissues were collected for subsequent proteomic analysis and/or targeted deep sequencing. Notably, 25 patients were excluded from this study due to a lack of primary lesions or not having PTC. (B) Grouping of included patients. Three samples are shown here for the Id, ID and iD groups. The first patient belonged to the Id group. Her or his lung metastases, which were confirmed by I-131 WBS and chest CT, were eliminated by two I-131 treatments. The second patient was included in the ID group. She or he received total thyroidectomy but her or his cervical regional metastatic lesion progressed during the four rounds of I-131 treatment. The first and second I-131 WBSs of the third patient revealed no positive focus. However, PET/CT during the same period revealed several suspicious pulmonary nodules (red arrow), which were confirmed to be metastases by subsequent surgery and histopathology. (C) Information on the collected FFPE samples. The protein profiles of 168 samples were evaluated in 10 batches. The sample types included primary lesions, synchronous and post I-131 cervical LNM, negative cervical lymph node and regional or distant metastatic lesions. Notably, the regional metastatic lesions were tumors near the primary site, such as the esophagus, trachea and cervical skin. LNM was not included in regional metastatic lesions. In addition, forty-eight primary tumors were examined via targeted deep sequencing.

^a The primary lesion of one patient was neither resected in our center nor examined to be PTC.

The included patients were then divided into two groups based on their response to RAI treatment: RAI-sensitive (RAI uptake positive and disease remission, Id) and RAI-refractory (disease persistence, Figure 6B). RAI refractoriness was identified in accordance with the 2015 American Thyroid Association (ATA) guideline ⁸. Briefly, patients with metastatic lesions that a) did not ever concentrate RAI, b) lost the ability to concentrate RAI after previous evidence of RAI avidity, and c) concentrated in some lesions but not in others were considered to have negative RAI uptake. Individuals with d) metastatic disease despite a significant concentration of RAI were defined as having positive uptake. Owing to the potential heterogeneity of RR-PTC, patients in the RAI-refractory group were further categorized into three subgroups: a) negative RAI uptake at the first RAI treatment with disease persistence (iDF), b) RAI uptake lost gradually after previous RAI treatments with disease persistence (iDG), and c) positive RAI uptake but with disease persistence (ID). Tumor stage and the risk of recurrence were stratified using the American Joint Committee on Cancer (AJCC) staging system and the 2015 ATA guideline ^8,63.

Resected primary and metastatic lesions were preserved in formalin-fixed paraffin-embedded (FFPE) tissue blocks at room temperature. A total of 168 FFPE samples were collected from the above 73 patients (Figure 6C). The tissue types included primary tumors, synchronous cervical LNM, LNM after I-131 treatment, negative lymph nodes, and regional and distant metastatic lesions. Hematoxylin and eosin-stained slides from all samples were reviewed by two independent pathologists with expertise in thyroid pathology. The tissue blocks with the highest tumor purity in each patient were selected for subsequent analysis. The borderlines of the tumor and adjacent tissues were carefully marked by expert pathologists.

4.3 Proteomic data acquisition and preprocessing

The FFPE samples were prepared for subsequent proteomic analysis as described previously ^26,64. The samples were allocated into ten batches. The tissues were subjected to a series of manipulations, including dewaxing, rehydration and lysis, for peptide extraction and digestion via pressure cycling technology (PCT). Peptides were then quantified via a liquid chromatography (LC) system coupled with a trapped ion mobility spectrometry mass spectrometer (MS). Data-independent acquisition (DIA) files were acquired and then analyzed against a thyroid tissue specific spectral library ²⁶ using DIA-NN (v1.8.1) ⁶⁵. Correlation analysis was conducted to assess the inter-batch and intra-batch stability of data acquisition. The profile was then processed for protein filtration, imputation and ID conversion. A detailed description of proteomic data acquisition, quality control and preprocessing can be found in the supplementary methods.

4.4 Targeted next-generation sequencing (TNGS)

The primary lesions from 48 patients with PTC were subjected to molecular profiling via TNGS. The protocol was performed as described previously ⁶⁶. Briefly, DNA was extracted from FFPE tissue sections using the QIAamp DNA FFPE Tissue Kit (Qiagen, Germany) following the manufacturer’s instructions. After amplification of targeted DNA fragments and removal of primers, the products were purified using an Ion AmpliSeq Library Kit (Thermo Fisher Scientific, USA). The concentrations of samples were quantitated by a NanoDrop system (Thermo Fisher Scientific, USA). The quality of the purified DNA was evaluated by 1% agarose gel electrophoresis. Severe degradation was detected in three samples, which were excluded from further analysis.

The remaining 45 library products were sequenced via 150 bp paired-end runs on the NextSeq 500 platform (Illumina, Inc., USA). The medians of sequencing depth and coverage were 5136× and 98.7%, respectively. Sequencing data were aligned to a reference human genome dataset (hg19/GRCh37). Subsequently, read mapping, quality control, variant calling and genotyping were performed following the protocols of the OncoAim® Thyroid Cancer Multigene Assay Kit (Singlera Genomics, Inc., China). Mutations and fusions were evaluated for 23 genes (Sup table S1). The minimum confidence threshold for variant allele frequency was 5%. The ENSEMBL Variant Effect Predictor (v90) was used for variant functional annotation. In addition, the gene variants were searched against the ClinVar database (v2020006) ⁶⁷. Somatic gene variants were categorized for clinical significance based on the Catalogue Of Somatic Mutations In Cancer database (COSMIC, v97) ⁶⁸.

4.5 Bioinformatic analysis

DEPs were identified using a multiple linear regression algorithm and analysis of variance (ANOVA). Differential expression analysis was conducted using the limma R package. Proteins with an absolute FC > 1.414 and a P value < 0.05 were defined as DEPs. The dimension of the profile was reduced and visualized using the PCA algorithm. ANOVA was used to identify proteins correlated with the severity of RR-PTC. Proteins with statistical significance were clustered using the mfuzz R package. Proteins in monotonically regulated clusters were regarded as DEPs.

Pathway enrichment analysis and gene ontology (GO) were conducted using four tools. GSEA was conducted using GSEA software to evaluate potential pathways and molecular mechanisms. For GSVA, ESs were calculated with the GSEA R package. Predefined gene sets were downloaded from the Reactome and KEGG databases. The top enriched GO processes were identified via the Metascape web-based platform. Another network tool, IPA, identifies most significantly relevant pathways with the overall activation or inhibition states based on DEPs.

Immune infiltration analysis was performed via the CIBERSORTx, ESTIMATE and single sample GSEA (ssGSEA) algorithms. Several scores, including the TDS, ERK score and CYT, were calculated based on the normalized protein profile. TIS and IIS were calculated using the Z score-standardized ssGSEA matrix. BRAF-RAS score (BRS) was not calculated since no RAS-mutated samples were identified.

WGCNA was performed to cluster genes with high correlation and assess the correlation between protein modules and gene alterations. Proteins in the modules that were highly associated with gene mutations or fusions were selected. Least absolute shrinkage and selection operator (LASSO) regression, random forest (RF), and support vector machine-recursive feature elimination (SVM-RFE) were used to identify biomarker proteins based on selected proteins.

Several external datasets were used for the validation of our results (Sup Table S7). GSE151179 from the GEO database included 52 samples derived from radioiodine-refractory and radioiodine-avid PTC patients. This dataset was used to verify the results related to RAI refractoriness. The associations between gene alterations and protein expression were assessed with the THCA (thyroid cancer) program from TCGA and four additional datasets from the GEO database. The predictive performance of the genes in the external datasets was evaluated with ROC curves.

A detailed description of the bioinformatic analysis and references to methodological articles can be found in the supplementary methods.

4.6 Statistical analysis

Quantitative variables were displayed as the means ± standard deviations or medians ± quartiles, whereas qualitative variables were presented as numbers and ratios. The significant differences in the quantitative variables were determined via two-tailed independent t test or Mann-Whitney U test, as appropriate. The chi-square test was used to evaluate the differences in the distributions of qualitative variables. When multiple groups were present, the z test with a Bonferroni correction was used to assess the intergroup difference in every group. Statistical analysis was conducted using SPSS software (IBM, US, v26).