The molecular and gene/miRNA expression profiles of radioiodine resistant papillary thyroid cancer
Journal of Experimental & Clinical Cancer Research volume 39, Article number: 245 (2020)
Papillary thyroid cancer (PTC) is the most frequent endocrine tumor. Radioiodine (RAI) treatment is highly effective in these tumors, but up to 60% of metastatic cases become RAI-refractory. Scanty data are available on either the molecular pattern of radioiodine refractory papillary thyroid cancers (PTC) or the mechanisms responsible for RAI resistance.
We analyzed the molecular profile and gene/miRNA expression in primary PTCs, synchronous and RAI-refractory lymph node metastases (LNMs) in correlation to RAI avidity or refractoriness.
We classified patients as RAI+/D+ (RAI uptake/disease persistence), RAI−/D+ (absent RAI uptake/disease persistence), and RAI+/D- (RAI uptake/disease remission), and analyzed the molecular and gene/miRNA profiles, and the expression of thyroid differentiation (TD) related genes.
A different molecular profile according to the RAI class was observed: BRAFV600E cases were more frequent in RAI−/D+ (P = 0.032), and fusion genes in RAI+/D+ cases. RAI+/D- patients were less frequently pTERT mutations positive, and more frequently wild type for the tested mutations/fusions. Expression profiles clearly distinguished PTC from normal thyroid. On the other hand, in refractory cases (RAI+/D+ and RAI−/D+) no distinctive PTC expression patterns were associated with either tissue type, or RAI uptake, but with the driving lesion and BRAF−/RAS-like subtype. Primary tumors and RAI-refractory LNMs with BRAFV600E mutation display transcriptome similarity suggesting that RAI minimally affects the expression profiles of RAI-refractory metastases. Molecular profiles associated with the expression of TPO, SLC26A4 and TD genes, that were found more downregulated in BRAFV600E than in gene fusions tumors.
The present data indicate a different molecular profile in RAI-avid and RAI-refractory metastatic PTCs. Moreover, BRAFV600E tumors displayed reduced differentiation and intrinsic RAI refractoriness, while PTCs with fusion oncogenes are RAI-avid but persistent, suggesting different oncogene-driven mechanisms leading to RAI refractoriness.
Radioactive iodine (RAI) therapy is a very effective treatment, significantly increasing the life expectancy of differentiated thyroid cancer (DTC) patients . However, about 60% of metastatic tumors (corresponding to less than 5% of all thyroid cancers) will become radioiodine-refractory (RAI-R), with a major impact on patient prognosis. In particular, the 10 years survival rate significantly decreases in patients without any 131I uptake at diagnosis compared to patients with initial 131I uptake but persistent diseases, and even more if compared to RAI-avid (RAI-A) metastatic thyroid tumors .
According to American Thyroid Association (ATA) Guidelines , RAI-refractory DTCs are identified by: a) presence of malignant/metastatic tissue that does not ever concentrate RAI, b) presence of tumor tissue that loses the ability to concentrate RAI after previous evidence of RAI-avid disease, c) RAI uptake in some lesions but not in others, and d) metastatic disease that progresses despite significant concentration of RAI.
It was shown that BRAFV600E mutated PTCs display a greater reduction in the expression of genes involved in iodine uptake and organification, namely sodium/iodide symporter NIS (SLC5A5), thyroperoxidase (TPO) and pendrin (SLC26A4), compared with tumors with other mutations or without known genetic alterations [4,5,6]. In vitro studies showed that the activation of BRAFV600E down-regulates the expression of NIS and reduces RAI uptake in thyroid cells. On the other hand, no difference either in NIS or in other thyroid differentiation genes was found in PTC with or without RET/PTC rearrangements , in accordance with the notion that these tumors rarely progress to aggressive or undifferentiated carcinomas .
Few studies are currently available on the molecular profile of RAI-R thyroid cancers, either well differentiated, poorly differentiated or anaplastic, showing that they are enriched with BRAF mutations, whereas RAS mutations are more represented in RAI-A metastatic DTCs [9, 10]. Based on these in vitro and in vivo data, pilot studies and phase 2 trials demonstrated that two selective BRAF inhibitors, namely Vemurafenib and Dabrafenib, are able to stimulate radioiodine uptake in patients with metastatic BRAFV600E-mutant RAI-R DTCs [11,12,13,14,15,16]. In addition, attempts of reverse RAI refractoriness in TC patients restoring radioiodine uptake, have been conducted with MEK inhibitors, tested as single agent  or, more recently, in combination with BRAF inhibitors [18, 19].
Even though the aberrant expression and/or function of NIS has been identified as factor for the lack of sensitivity to RAI, the mechanisms responsible for RAI resistance are still poorly understood. To get more insights into this topic, we investigated both the molecular profiles of RAI-A and RAI-R PTC patients, and their gene/miRNA expression. Indeed, though some studies have already described miRNA deregulation in PTC suggesting their possible role as diagnostic and prognostic markers [20,21,22,23,24], any report is currently available about miRNA expression in relation to the response to RAI.
Therefore, aim of this study was to find possible molecular tags able to identify the avidity or refractoriness to radioiodine of PTCs, since this information is expected to have a major impact on the selection of the more appropriate treatment and follow up .
We retrospectively analyzed, the clinical records of 1339 PTC patients treated and followed by our tertiary care center during the period 2001–2019 (Fig. 1). Strict selection criteria were applied to identify patients for whom RAI avidity or refractoriness could be definitely determined: first: patients without regional or distant metastases at diagnosis were excluded, due to the impossibility to determine RAI avidity or refractoriness (1339–886 = 453); second: patients with regional or distant metastases at diagnosis but not submitted to RAI treatment for several reasons (patient choice in most cases) were excluded (453–59 = 394); third: patients with RAI uptake limited to the thyroid bed at the time of ablation with following complete remission were excluded (394–324 = 70).
Finally, the selected 70 PTC patients, all with regional or distant metastases at diagnosis, were classified into three different groups according to RAI uptake and response: RAI+/D+: 23 PTCs with initial RAI uptake at the metastatic site, and disease persistence; RAI−/D+: 32 PTCs without RAI uptake at the metastatic site and disease persistence; and RAI+/D-: 15 PTCs with RAI uptake at the metastatic site and disease remission. In RAI+ groups, all metastases, either loco regional or distant, were RAI avid.
Patients were treated, followed up, and defined as in remission or in persistence according to the International and National guidelines [3, 26]. Patients were risk-stratified using AJCC staging system (Stage I, II,III or IV)  and 2015 ATA Guidelines (low, intermediate or high risk of recurrence) .
Molecular profile by PTC-MA assay
Seventy primary PTCs (52 frozen and 18 FFPE tissues) and 22 lymph node metastases (frozen tissues) were submitted to molecular analyses by PTC-Mass Array platform (PTC-MA), previously set up for the simultaneous identification of 13 known hotspot mutations (BRAFV600E; AKT1E17K; EIF1AX c.338-1G > C; NRASQ61R and NRASQ61K; HRASG13C, HRASQ61K, and HRASQ61R; KRASG12V and KRASG13C; TERT c.-124C > T and TERT c.-146C > T; and PIK3CAE542K) and 6 recurrent fusion genes typical of PTC: RET/PTC1 (RET/CCDC6), RET/PTC2 (RET/PRKAR1A), and RET/PTC3 (RET/NCOA4) and TRK (NTRK1/TPM3), TRK-T1 (NTRK-T1/TPR), and TRK-T3 (NTRK1/TFG) [28, 29].
For point mutations the allelic frequencies were recorded and normalized for the cancer cell content, as previously described . For fusion genes, the allelic frequency cannot be evaluated, since their detection involves a selective amplification of the rearranged gene transcript.
Samples selection for gene/miRNA profiling
For 33/70 patients submitted to PTC-MA analysis, frozen samples of the primary tumor and/or the metastasis were still available after the genetic characterization, and were submitted to gene/miRNA microarray analysis. Whenever available, synchronous lymph node metastasis (synLNM) and/or RAI-R LNM diagnosed and collected after RAI treatment (LNM post RAI), and non-neoplastic thyroid (NT) were also studied. In order to increase the number of metastatic samples, 9 additional patients (for a total of 42 patients) were included, for whom only metastatic tumor tissue was available. Overall, 69 frozen tissues derived from those 42 PTC metastatic patients were collected including 23 primary tumors, 8 synLNMs, 20 LNMs post RAI, and 18 matched NTs (Fig. 1).
Gene/miRNA microarray analysis
Nucleic acid extraction was performed on the 69 collected frozen tissues. Genomic DNA and RNA were extracted by Qiagen DNeasy Blood & Tissue Kit and by Qiagen miRNeasy Mini Kit, respectively, and quantified by Qubit 2.0 Fluorometer (Thermo Fisher). DNA and RNA quality and integrity were assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA).
Gene profiles were established by Thermo Fisher Human Clariom S Assay. MiRNA profiles were established by Agilent SurePrint Human miRNA microarrays. Complete microarrays experimental procedures and data processing available in Additional file 2: Supplemental Methods.
Following quality control and data pre-processing, gene expression data were available for 52 thyroid tissues including 39 primary and LN metastatic PTCs and 13 NTs (Fig. 1).
The BRAF−/RAS-like subtype, indicative of MAPK pathway transcriptional activation, was assigned using the 71 gene signature derived from The Cancer Genome Atlas (TCGA) study on PTC  as previously described ; BRAF−/RAS-like subtype were further confirmed using the Nearest Template Prediction (NTP) algorithm  implemented in GenePattern software (https://www.genepattern.org/).
Unsupervised analyses for gene and miRNA arrays were performed using principal component analysis (PCA) and hierarchical clustering with Euclidean distance and Ward linkage method. Heatmaps were visualized using the ComplexHeatmap package . Differential expression analysis was carried out using the limma package . P-values were corrected for multiple testing using the Benjamini-Hochberg false discovery rate method.
CIBERSORT algorithm  was applied to compute tumor purity; p-value outputs were used to stratify samples (p-value < 0.05 were considered significantly infiltrated and assigned to a low tumor purity class, while p-value ≥0.05 were assigned to a high tumor purity class). According to CIBERSORT tumor purity, 21 samples, resulting with low purity, were excluded from the subsequent analysis (Fig. 1).
Gene expression validation in public datasets
Human thyroid cancer gene datasets available on public repositories were investigated for genes expression validation. TCGA study on PTC  and nine additional datasets from GEO were analyzed (GSE27155 , GSE3467 , GSE6004 , GSE33630 [39, 40], GSE35570 , GSE53157 , GSE29265, GSE3678, and GSE60542 ). Details available in Additional file 2: Supplemental Methods. Both normal thyroid and primary PTC samples were tested. Three thyroid specific genes (NIS/SLC5A5, TPO, and SLC24A6) and a thyroid differentiation (TD) score were assessed; only PTCs with reported BRAFV600E mutation or RET and NTRK1 gene fusions were specifically studied. TD score was calculated as mean of log2-transformed and median-centered expression of 16 thyroid function related genes as previously described in TCGA study ; list of the 16 TD genes available in Additional file 1: Suppl. Figure S6E.
Relations between discrete variables were evaluated by means of χ2 test or t-test, as appropriate. Clinicopathological and molecular features were evaluated by univariate analysis.
Statistical significance was defined as P < 0.05. All statistical analyses were performed using MedCalc Analyses using the Version 18.11.3 of the MedCalc Software (B-8400 Ostend, Belgium).
The three groups of patients (RAI+/D+, RAI−/D+, and RAI+/D-) displayed comparable clinicopathological features at diagnosis, with minor differences observed only for the prevalence of the non-conventional histological variants of PTC. There were no differences regarding the most important prognostic criteria of differentiated thyroid cancer: TNM, AJCC stage, ATA risk classes (Fig. 2). The median cumulative dose of radioiodine was similar in the three groups (Table 1).
According to patients’ classification, RAI-R patients (RAI+/D+, RAI−/D+) were in disease persistence at the last follow-up, while RAI-A patients (RAI+/D-) were all in remission after radiometabolic treatment. In addition, the two groups of RAI-R patients had comparable prevalence of both disease progression (RAI+/D+, 3/23 (13%) vs. RAI−/D+, 8/32 (25%)) and stable disease persistence (RAI+/D+, 20/23 (87%) vs. RAI−/D+, 24/32, (75%)) (P = 0.278).
The molecular profile of PTCs correlated to radioiodine resistance
Analysis in 70 primary PTCs
The molecular profile was different in the 3 groups of patients: BRAFV600E was significantly more frequent in RAI−/D+ group (66%) compared with RAI+/D+ (35%) and RAI+/D- (33%) groups (P = 0.032). Consistent with previous findings on the association between RAS mutations and RAI-avidity, RAS mutations were identified in RAI+/D- (7%), and in the RAI+/D+ group (9%), but not in the RAI−/D+ group (P = 0.82, Fig. 2).
The prevalence of the other genetic alterations was not significantly different among the 3 groups, but it is worth to note that, though not statistically significant, probably due to the small number of cases, fusion genes (in particular RET/PTC) were more frequent in RAI+/D+ (39%), compared with RAI−/D+ (16%) and RAI+/D- (13%) patients (P = 0.075). TERT promoter (pTERT) mutations had a higher prevalence, though not statistically significant, in the group without RAI uptake (RAI−/D+, 28%) compared with the RAI+/D+ and RAI+/D- groups (17 and 20%) (P = 0.61).
Finally, the prevalence of PTCs that were wild type for all the analyzed point mutations and fusions was higher in RAI-avid group (RAI+/D-, 40%) than in RAI-refractory PTCs (RAI+/D+ 17% and RAI−/D+ 22%, P = 0.41) (Fig. 2).
In 9/11 patients with available multiple tumor tissues we found concordant genotypes (Fig. 3a), with either BRAFV600E mutations or gene fusions detected both at the primary tumor and at the metastatic sites.
No statistically significant differences were found in the allelic frequencies of point mutations among the 3 groups analyzed: the mean ± SD normalized allelic frequencies were 23.3 ± 10.1 in RAI+/D+, 30 ± 9.2 in RAI−/D+ and 35 ± 9.4 in RAI+/D- for BRAF (P = 0.13); 37 ± 15.1 in RAI+/D+, 41.5 ± 10.2 in RAI−/D+ and 47 ± 32 in RAI+/D- for pTERT (P = 0.37) and 33 ± 1.4 in RAI+/D+ and 32 in RAI+/D- for RAS mutations.
Analysis in the subgroup for gene/miRNA studies
The 39 tumor samples that underwent gene/miRNA analysis were also assessed by PTC-MA profiling (Fig. 3b). They included both primary (n = 17) and lymph node metastases (n = 22, 5 collected before and 17 after RAI treatment). Of note, the 17 primary tumors included in this set displayed a genotype-RAI class distribution, with BRAFV600E mutation more frequent in RAI−/D+ patients, and fusion genes in RAI+/D+ patients, similar to that observed in the 70 primary PTCs (Fig. 2), thus representing a reliable subgroup.
In LNMs collected after RAI (LNM post RAI) BRAFV600E mutation was more frequent in RAI−/D+ (58% vs. 20%, P = 0.14), while fusion genes in RAI+/D+ patients (40% vs. 0%, P = 0.09), consistently with the results found in primary PTCs. A similar enrichment was observed also in synLNMs, even though only few samples were tested. Interestingly, pTERT mutations were detected only in LNMs post RAI, while no RAS mutations were found (Fig. 3b). Overall, driving genetic alterations were identified in 28 (72%) of the 39 tumor tissues.
Gene/miRNA expression in RAI-refractory and RAI-avid PTCs
Gene and miRNA expression was thus investigated in the above described series of 39 samples, including primary PTCs and metastatic LN, and in 13 available matched NTs, included as control (Fig. 1).
First, by using a TCGA derived gene signature  we tested the BRAF−/RAS-related signaling and confirmed that BRAFV600E mutated tumors, independently of the tissue type, were BRAF-like, as well as most of the detected gene fusions, and that both BRAF- and RAS-like subtypes were present in WT samples (Additional file 1: Suppl. Figure S1).
Considering gene profiles, we found, as expected, a significant separation between tumor and NT tissues, both via unsupervised hierarchical clustering and principal component analysis (PCA; Additional file 1: Suppl. Figure S2A and E). Focusing on PTCs, we did not observe a clear separation considering either tumor tissue type (Additional file 1: Suppl. Figure S2A and E), or collection before or after RAI treatment (Additional file 1: Suppl. Figure S2B), or patient RAI class (Additional file 1: Suppl. Figure S2B and S2F), or driving lesion (Additional file 1: Suppl. Figure S2C), but tumor samples were distinctly separated in two major clusters (Tumor cluster 1 and 2, Additional file 1: Suppl. Figure S2A). The two groups had a different tissue composition in term of cancer cells content. Different tumor purity was indeed found by the gene expression based algorithm CIBERSORT (Additional file 1: Suppl. Figure S2D and G) and by tissue histological review (Additional file 1: Suppl. Figure S2H), that confirmed in low purity tissues a reduced cancer cell content and increased presence of microenvironment derived cells. Interestingly, low tissue purity was observed also in a subgroup of NTs (NT 2 cluster; Additional file 1: Suppl. Figure S2D and G); clinical records and histological revision confirmed the presence of an autoimmune thyroiditis in these patients (Additional file 1: Suppl. Figure S2I). Of note, also the primary tumors included in the low purity group (Tumor cluster 2) displayed at histology a more frequently associated thyroiditis (Tumor cluster 2, 63%, vs. Tumor cluster 1, 11%, P = 0.026¸ Additional file 1: Suppl. Figure S2I).
As in gene expression studies reduced tumor cell content can negatively affect the analysis of clinical tumor specimens, since the presence of infiltrating cells can cause the dilution of tumor-specific mRNAs and miRNAs, we decided to focus only on high-purity tissues (based on CIBERSORT class; Additional file 1: Suppl. Figure S2D and G). In addition, as these latter included only two PTCs from RAI-avid patients and could be only partially representative of this group, we excluded these two samples from the subsequent analyses and focused on RAI-refractory cases, as below described.
Gene/miRNA expression in RAI-refractory primary and metastatic PTCs
Gene/miRNA expression was thus specifically reanalyzed in RAI-refractory patients (RAI+/D+ and RAI−/D+), comprising a total of 29 tissues (10 primary and 10 metastatic LN, and 9 NTs; Fig. 1).
Considering gene expression, we confirmed also in this set a significant separation between tumor and NT (Fig. 4a and e). No clear separation was found between the different tumor tissues (Fig. 4a), even though the majority of LNMs post RAI (5/7, 72%) localized in the same cluster (Fig. 4e). No significant separation was observed according to patient RAI-refractoriness class (Fig. 4b and f). On the other hand, a significant stratification was observed according to the BRAF−/RAS-subtype (Fig. 4c) and to the driver genetic alteration (Fig. 4d). In particular, RAI-refractory BRAFV600E tumors appeared as a quite homogenous group (Fig. 4d) with the majority of samples clustering together (8/11, 73%; Fig. 4h) regardless of either the tissue type (i.e. primary tumors and LFN post RAI; Fig. 4e) or the collection before or after RAI treatment (Fig. 4g).
Similar results were obtained by miRNA profiles. Significant miRNA deregulation was observed in PTCs compared with NTs, confirming the upregulation of broadly reported oncogenic miRNAs (miR-146b-5p, miR-222-3p, miR-221-3p, miR-21-5p/−3p) and the downregulation of miR-7-2-3p/−5p and miR-486-5p (Fig. 4l). However, no significant differences were found in their expression across tumor samples either based on tissue type, or patient RAI class or collection before or after RAI treatment (Fig. 4e-g). Along with these, we confirmed the up- and down-modulation of additional miRNAs (Additional file 1: Suppl. Figure S3), whose deregulation has been already reported in PTC (revised in ). Considering these PTC related miRNAs (Additional file 1: Suppl. Figure S3C), we observed samples stratification according to BRAFV600E mutation (Additional file 1: Suppl. Figure S3D) and BRAF−/RAS-like subtype (Additional file 1: Suppl. Figure S3E). Interestingly, we found that while miRNA upregulation was consistent across the various driving lesions and BRAF−/RAS-subtypes, miRNA downregulation was more evident in BRAFV600E and in BRAF-like tumors (Additional file 1: Suppl. Figure S3H-I).
In our series of RAI refractory PTCs along with miRNA, we found consistent gene deregulation. Many genes were differentially expressed in PTCs compared with NTs (Additional file 1: Suppl. Figure S4A), and their deregulation was confirmed in the large PTC series from TCGA (Additional file 1: Suppl. Figure S4B). Upregulated genes were associated with cell junction/adhesion, integrin and extracellular matrix remodeling, and inflammation (Additional file 1: Suppl. Figure S5A); they included also MAPK pathway regulators (Additional file 1: Suppl. Figure S5B) and genes of BRAF−/RAS-like subtype signature. Downregulated genes were involved in metabolic processes, and, more interestingly, included pendrin (SLC26A4), thyroid peroxidase (TPO) and several other genes related to thyroid function and differentiation (Additional file 1: Suppl. Figure S5C-D).
As the expression of NIS (SLC5A5), TPO and SLC26A4 has been already proposed to affect RAI uptake and response, we tested whether their expression could be different in our PTC samples stratified for RAI class of refractoriness (Additional file 1: Suppl. Figure S6A). Interestingly, a significant difference was observed for SLC26A4 gene. Moreover, we tested the possible association between gene expression and molecular profile of our RAI-refractory cases. We found that TPO and SLC26A4 genes were consistently different between BRAFV600E and gene fusions positive tumors (Fig. 5a), while, by contrast, NIS expression appeared to be similarly reduced. This was confirmed in other PTC series derived from public repositories (GEO datasets, Fig. 5c and e; and TCGA dataset, Fig. 5g), where significant differences were observed for TPO and SLC26A4 genes, whereas NIS expression was significantly reduced only in TCGA dataset. Along with these, other genes related to thyroid metabolism and function exist and TCGA proposed a Thyroid Differentiation (TD) score  that summarizes the expression of 16 thyroid related genes (also including the above mentioned SLC5A5, TPO, and SLC26A4). We applied the TCGA established TD score and not only confirmed its expression in PTC tissues (Additional file 1: Suppl. Figure S6C), but also its high level in normal thyroid, which was unexplored for TD score in TCGA study (Additional file 1: Suppl. Figure S6D). The same TD score was applied to our samples series (Additional file 1: Suppl. Figure S6E). While no significant differences were observed according to RAI refractoriness class (Additional file 1: Suppl. Figure S6B), it was significantly different in BRAFV600E compared to gene fusions, with BRAFV600E PTCs displaying reduced thyroid differentiation levels (Fig. 5b). The same results were obtained in public datasets (Fig. 5d, f and h). In addition, TD score was significantly associated with sample clustering observed in both gene and miRNA expression profiles (Fig. 4m and Additional file 1: Suppl. Figure 3G).
In this study we investigated the molecular profiles and gene/miRNA expression in metastatic PTC patients stratified for RAI response in accordance with their RAI-refractoriness/avidity.
By the PTC-MA platform we found a different molecular profile in RAI-avid and RAI-refractory PTCs. BRAFV600E was significantly more frequent in PTCs unable to uptake RAI at the metastatic site already at the time of the diagnosis. An association between BRAF mutation and RAI-R disease has been previously reported [4,5,6, 9], though a whole molecular profile was not obtained and cases were not classified according to the “type” of RAI refractoriness, evident at the time of initial treatment or developed during time. Our data indicate that BRAFV600E associates with an intrinsic tumor RAI refractoriness, being significantly more prevalent in tumors RAI−/D+. On the contrary, we report for the first time a higher prevalence of fusion genes in PTCs with initial RAI uptake, but without therapeutic efficacy (RAI+/D+), suggesting that also these molecular alterations interfere with the iodine organification process, though not impairing iodine uptake. Accordingly, the partial response recently reported in an advanced DTC harboring an EML4-NTRK3 fusion and treated with Larotrectinib, a NTRK inhibitor, is consistent with the restoration of radioiodine organification . pTERT mutations were more frequent in the PTC group without RAI-uptake, in accordance with data that report for this genetic alteration a negative prognostic role in thyroid cancer . On the other hand, cases wild type for all the tested mutations and gene fusions were more prevalent in patients with RAI-avid metastatic disease in remission, indicating that gene alterations not analyzed by our assay are more frequently associated with a more differentiated and less aggressive disease.
In patients with available multiple tumor specimens, we found in most cases concordant genetic patterns between primary PTC and lymph node metastases, not only for BRAFV600E mutation, in agreement with previous reports , but also for RET/PTC fusions. This was observed not only in synchronous LNMs but also in RAI-refractory LNM collected after RAI treatment.
The mutational profiles of primary and synchronous metastases in thyroid cancers, including PTC, have been already investigated [46,47,48,49,50], describing high and significant concordance between the primary tumor and the corresponding metastasis, especially the LNM. On the other hand, only few studies assessed the global gene/miRNA expression in these tissue types [43, 51, 52], confirming transcriptome similarity between primary TC and matched synchronous LNMs. To our knowledge, this is the first study investigating gene/miRNA expression in PTC tissues including RAI-refractory LNMs collected following RAI treatment.
We initially investigated a series including primary and LN metastatic PTCs, from either RAI-refractory or -avid patients. The gene profiles identified a clear separation between neoplastic and normal tissues, and tumor samples were also found to stratify into 2 major clusters, according to the tissues composition, and in particular to presence or not of an associated autoimmune thyroiditis or other types of immune cells infiltrating the tissue for primary tumors and NTs, and, of residual lymphoid tissue not replaced by neoplastic cells for LNMs. Indeed, PTC tissues including very high or low amount of immune cells displayed a distinguishable transcriptomic profile, which should be further evaluated for possible clinical impact in larger series.
On the other hand, the presence of these tumor infiltrating cells, of various origin, might have affected our analyses impairing a clear identification of PTC-specific mRNAs and miRNAs, as already reported [31, 43, 53,54,55].
With the aim of identifying tumor related features, that could be not perturbed by the presence of other microenvironment derived cellular components, we thus focused on the subset of samples displaying high tissue purity from RAI-refractory patients.
In this sample set gene profiles confirmed not only the separation between tumors and NTs, but also and more interestingly, the transcriptome similarity among BRAFV600E mutated samples, either primary or LNMs post RAI, thus suggesting that radioiodine treatment minimally affects the expression profiles of RAI-refractory lymph metastases with respect to primary tumors.
In this samples group we also confirmed the up- and down- regulation of previously reported PTC-related miRNAs (revised in ). In addition, we found a correlation with genotype and a more pronounced miRNA downregulation in BRAFV600E and BRAF-like tumors, suggestive of more dedifferentiated tissues.
To the best of our knowledge this is the first study investigating gene/miRNA expression in tumor samples from RAI refractory patients, and including not only primary PTCs, but also synchronous and RAI-refractory LNMs. We found gene and miRNA deregulations consistent with those described by previous reports, such as the TCGA study , which analyzed primary PTCs, regardless of the RAI classification. In addition, according with TCGA , that describes how PTC subclassification is influenced by BRAF−/RAS-subtype and by the expression of thyroid related genes (TD score), along with the histological variant and driving lesion, we found that in our series from RAI-refractory patients (RAI+/D+ and RAI−/D+) the major classifying features are still the driving lesion (and in particular BRAFV600E), the BRAF−/RAS-subtype and the TD score.
Among TD genes, NIS was found to be significantly dowregulated in RAI-refractory tumors, and particularly in RAI−/D+ cases, consistent with the clinical behavior. As an original finding, also SLC26A4 was significantly more downregulated in RAI−/D+ than in RAI+/D+ cases, suggesting that these cases are even more dedifferentiated and that if a small amount of RAI is uptaken, it is immediately washed out due to the lack of the apical transporter . Finally, TPO was downregulated in both RAI refractory classes, as expected for tissues not able to organify and thus retain iodine into the colloid. By evaluating these crucial genes in relation to the genetic pattern, we found that NIS appeared to be similarly reduced in BRAFV600E and fusions positive tumors, whereas TPO, SLC26A4 and TD score were significantly more downregulated in BRAFV600E tumors than in fusion positive tumors. Interestingly, these data were further confirmed by analyzing other series from public repositories and are consistent with the already reported dedifferentiation  and RAI refractoriness associated with BRAFV600E mutation [9, 10].
Further studies to get more insights into the functional relationships among BRAF and RET/PTC oncogenes, RAI uptake and expression of thyroid differentiation genes such as NIS and SLC26A4, would be useful. In this context it was demonstrated that in mouse thyroid cancers with conditional BRAFV600E expression the lack of radioiodine incorporation depend on BRAF activity . In the clinical setting, promising data on redifferentiation have been obtained with new compounds targeting BRAF or MEK, such as Vemurafenib, Dabrafenib, Trametinib, and Selumetinib, though the anti-tumor effect is variable, ranging from a partial to a null response, possibly due to the histopathological and genetic heterogeneity of mutated tumors [11,12,13,14,15,16,17,18,19]. Finally, based on the original data obtained in the present work, future studies will be particularly focused on the possible effects of RET/PTC oncogenes on RAI uptake, since to the best of our knowledge no data are available on this topic.
We studied the molecular profile and gene/miRNA expression in primary PTCs, synchronous and RAI-refractory LNMs in correlation to RAI avidity or refractoriness.
BRAFV600E tumors were found more frequently unable to uptake RAI at the metastatic site already at the time of the diagnosis, suggesting an intrinsic RAI refractoriness. Accordingly, the genes related to RAI refractoriness, e.g. those responsible for iodine apical transport or organification, were highly downregulated, regardless of either the different tissue type (i.e. primary tumor, or synchronous or RAI refractory lymph nodal metastases) or the collection before or after RAI treatment, indicating that the dedifferentiation related to this mutation is present at the very early stages of tumor growth. Moreover, a transcriptome similarity was found among BRAFV600E mutated samples, either primary or lymph nodal metastases post RAI, thus suggesting that radioiodine treatment minimally affects the expression profiles of RAI-refractory lymph metastases with respect to primary tumors.
On the other hand, genes related to RAI refractoriness were less downregulated in fusions positive tumors which were, consistently, more frequently able to uptake RAI, but still with persistent disease, suggesting different mechanisms leading to RAI-refractoriness.
Availability of data and materials
The datasets generated and/or analysed during the current study are available in the NCBI Gene Expression Omnibus (GEO) repository (www.ncbi.nlm.nih.gov/geo/) with the accession number GSE151179. Gene microarray data include 52 thyroid tissues (39 primary and LN metastatic PTCs and 13 NTs); the corresponding miRNA data were available for 47/52 samples.
Differentiated thyroid cancer
Gene Expression Omnibus
Lymph node metastases
- LNM post RAI:
Radioiodine-refractory LNM diagnosed and collected after RAI treatment
Papillary thyroid cancer
Absent RAI uptake/disease persistence
RAI uptake/disease remission
RAI uptake/disease persistence
Synchronous lymph node metastasis
The Cancer Genome Atlas
Verburg FA, Hanscheid H, Luster M. Radioactive iodine (RAI) therapy for metastatic differentiated thyroid cancer. Best Pract Res Clin Endocrinol Metab. 2017;31:279–90.
Durante C, Haddy N, Baudin E, Leboulleux S, Hartl D, Travagli JP, et al. Long-term outcome of 444 patients with distant metastases from papillary and follicular thyroid carcinoma: benefits and limits of radioiodine therapy. J Clin Endocrinol Metab. 2006;91:2892–9.
Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, et al. 2015 American Thyroid Association management guidelines for adult patients with thyroid nodules and differentiated thyroid Cancer: the American Thyroid Association guidelines task force on thyroid nodules and differentiated thyroid Cancer. Thyroid. 2016;26:1–133.
Durante C, Puxeddu E, Ferretti E, Morisi R, Moretti S, Bruno R, et al. BRAF mutations in papillary thyroid carcinomas inhibit genes involved in iodine metabolism. J Clin Endocrinol Metab. 2007;92:2840–3.
Mian C, Barollo S, Pennelli G, Pavan N, Rugge M, Pelizzo MR, et al. Molecular characteristics in papillary thyroid cancers (PTCs) with no 131I uptake. Clin Endocrinol. 2008;68:108–16.
Bastos AU, Oler G, Nozima BH, Moyses RA, Cerutti JM. BRAF V600E and decreased NIS and TPO expression are associated with aggressiveness of a subgroup of papillary thyroid microcarcinoma. Eur J Endocrinol. 2015;173:525–40.
Romei C, Ciampi R, Faviana P, Agate L, Molinaro E, Bottici V, et al. BRAFV600E mutation, but not RET/PTC rearrangements, is correlated with a lower expression of both thyroperoxidase and sodium iodide symporter genes in papillary thyroid cancer. Endocr Relat Cancer. 2008;15:511–20.
Nikiforov YE. RET/PTC rearrangement in thyroid tumors. Endocr Pathol. 2002;13:3–16.
Ricarte-Filho JC, Ryder M, Chitale DA, Rivera M, Heguy A, Ladanyi M, et al. Mutational profile of advanced primary and metastatic radioactive iodine-refractory thyroid cancers reveals distinct pathogenetic roles for BRAF, PIK3CA, and AKT1. Cancer Res. 2009;69:4885–93.
Sabra MM, Dominguez JM, Grewal RK, Larson SM, Ghossein RA, Tuttle RM, et al. Clinical outcomes and molecular profile of differentiated thyroid cancers with radioiodine-avid distant metastases. J Clin Endocrinol Metab. 2013;98:E829–36.
Brose MS, Cabanillas ME, Cohen EE, Wirth LJ, Riehl T, Yue H, et al. Vemurafenib in patients with BRAF(V600E)-positive metastatic or unresectable papillary thyroid cancer refractory to radioactive iodine: a non-randomised, multicentre, open-label, phase 2 trial. Lancet Oncol. 2016;17:1272–82.
Dunn LA, Sherman EJ, Baxi SS, Tchekmedyian V, Grewal RK, Larson SM, et al. Vemurafenib Redifferentiation of BRAF mutant, RAI-refractory thyroid cancers. J Clin Endocrinol Metab. 2019;104:1417–28.
Dadu R, Shah K, Busaidy NL, Waguespack SG, Habra MA, Ying AK, et al. Efficacy and tolerability of vemurafenib in patients with BRAF(V600E) -positive papillary thyroid cancer: M.D. Anderson Cancer center off label experience. J Clin Endocrinol Metab. 2015;100:E77–81.
Aashiq M, Silverman DA, Na'ara S, Takahashi H, Amit M. Radioiodine-refractory thyroid Cancer: molecular basis of Redifferentiation therapies, management, and novel therapies. Cancers (Basel). 2019;11:1382.
Rothenberg SM, McFadden DG, Palmer EL, Daniels GH, Wirth LJ. Redifferentiation of iodine-refractory BRAF V600E-mutant metastatic papillary thyroid cancer with dabrafenib. Clin Cancer Res. 2015;21:1028–35.
Falchook GS, Millward M, Hong D, Naing A, Piha-Paul S, Waguespack SG, et al. BRAF inhibitor dabrafenib in patients with metastatic BRAF-mutant thyroid cancer. Thyroid. 2015;25:71–7.
Ho AL, Grewal RK, Leboeuf R, Sherman EJ, Pfister DG, Deandreis D, et al. Selumetinib-enhanced radioiodine uptake in advanced thyroid cancer. N Engl J Med. 2013;368:623–32.
Iravani A, Solomon B, Pattison DA, Jackson P, Ravi KA, Kong G, et al. Mitogen-activated protein kinase pathway inhibition for Redifferentiation of radioiodine refractory differentiated thyroid Cancer: an evolving protocol. Thyroid. 2019;29:1634–45.
Leboulleux S, Dupuy C, Lacroix L, Attard M, Grimaldi S, Corre R, et al. Redifferentiation of a BRAF(K601E)-mutated poorly differentiated thyroid Cancer patient with Dabrafenib and Trametinib treatment. Thyroid. 2019;29:735–42.
Boufraqech M, Klubo-Gwiezdzinska J, Kebebew E. MicroRNAs in the thyroid. Best Pract Res Clin Endocrinol Metab. 2016;30:603–19.
Fuziwara CS, Kimura ET. MicroRNAs in thyroid development, function and tumorigenesis. Mol Cell Endocrinol. 2017;456:44–50.
Wojcicka A, Kolanowska M, Jazdzewski K. Mechanisms in endocrinology: MicroRNA in diagnostics and therapy of thyroid cancer. Eur J Endocrinol. 2015;174:R89–98.
Minna E, Romeo P, De CL, Dugo M, Cassinelli G, Pilotti S, et al. miR-199a-3p displays tumor suppressor functions in papillary thyroid carcinoma. Oncotarget. 2014;5:2513–28.
Minna E, Romeo P, Dugo M, De CL, Todoerti K, Pilotti S, et al. miR-451a is underexpressed and targets AKT/mTOR pathway in papillary thyroid carcinoma. Oncotarget. 2016;7:12731–47.
Fugazzola L, Elisei R, Fuhrer D, Jarzab B, Leboulleux S, Newbold K, et al. 2019 European thyroid association guidelines for the treatment and follow-up of advanced radioiodine-refractory thyroid Cancer. Eur Thyroid J. 2019;8:227–45.
Pacini F, Basolo F, Bellantone R, Boni G, Cannizzaro MA, De PM, et al. Italian consensus on diagnosis and treatment of differentiated thyroid cancer: joint statements of six Italian societies. J Endocrinol Investig. 2018;41:849–76.
Edge SB, Byrd DR, Compton CC, Fritz AG, Greene FL, Trotti A. American Joint Committee on Cancer (AJCC) cancer staging manual. 7th ed. Chicago: Springer, Inc; 2010.
Pesenti C, Muzza M, Colombo C, Proverbio MC, Fare C, Ferrero S, et al. MassARRAY-based simultaneous detection of hotspot somatic mutations and recurrent fusion genes in papillary thyroid carcinoma: the PTC-MA assay. Endocrine. 2018;61:36–41.
Colombo C, Muzza M, Proverbio MC, Tosi D, Soranna D, Pesenti C, et al. Impact of mutation density and heterogeneity on papillary thyroid Cancer clinical features and remission probability. Thyroid. 2019;29:237–51.
The Cancer Gene Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159:676–90.
Minna E, Brich S, Todoerti K, Pilotti S, Collini P, Bonaldi E, et al. Cancer associated fibroblasts and senescent thyroid cells in the invasive front of thyroid carcinoma. Cancers (Basel). 2020;12:112.
Hoshida Y. Nearest template prediction: a single-sample-based flexible class prediction with confidence assessment. PLoS One. 2010;5:e15543.
Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32:2847–9.
Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK. ROBUST Hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Ann Appl Stat. 2016;10:946–63.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7.
Giordano TJ, Kuick R, Thomas DG, Misek DE, Vinco M, Sanders D, et al. Molecular classification of papillary thyroid carcinoma: distinct BRAF, RAS, and RET/PTC mutation-specific gene expression profiles discovered by DNA microarray analysis. Oncogene. 2005;24:6646–56.
He H, Jazdzewski K, Li W, Liyanarachchi S, Nagy R, Volinia S, et al. The role of microRNA genes in papillary thyroid carcinoma. Proc Natl Acad Sci U S A. 2005;102:19075–80.
Vasko V, Espinosa AV, Scouten W, He H, Auer H, Liyanarachchi S, et al. Gene expression and functional evidence of epithelial-to-mesenchymal transition in papillary thyroid carcinoma invasion. Proc Natl Acad Sci U S A. 2007;104:2803–8.
Dom G, Tarabichi M, Unger K, Thomas G, Oczko-Wojciechowska M, Bogdanova T, et al. A gene expression signature distinguishes normal tissues of sporadic and radiation-induced papillary thyroid carcinomas. Br J Cancer. 2012;107:994–1000.
Tomas G, Tarabichi M, Gacquer D, Hebrant A, Dom G, Dumont JE, et al. A general method to derive robust organ-specific gene expression-based differentiation indices: application to thyroid cancer diagnostic. Oncogene. 2012;31:4490–8.
Handkiewicz-Junak D, Swierniak M, Rusinek D, Oczko-Wojciechowska M, Dom G, Maenhaut C, et al. Gene signature of the post-Chernobyl papillary thyroid cancer. Eur J Nucl Med Mol Imaging. 2016;43:1267–77.
Pita JM, Banito A, Cavaco BM, Leite V. Gene expression profiling associated with the progression to poorly differentiated thyroid carcinomas. Br J Cancer. 2009;101:1782–91.
Tarabichi M, Saiselet M, Tresallet C, Hoang C, Larsimont D, Andry G, et al. Revisiting the transcriptional analysis of primary tumours and associated nodal metastases with enhanced biological and statistical controls: application to thyroid cancer. Br J Cancer. 2015;112:1665–74.
Groussin L, Clerc J, Huillard O. Larotrectinib-enhanced radioactive iodine uptake in advanced thyroid Cancer. N Engl J Med. 2020;383:1686–7.
Yang X, Li J, Li X, Liang Z, Gao W, Liang J, et al. TERT promoter mutation predicts radioiodine-refractory character in distant metastatic differentiated thyroid Cancer. J Nucl Med. 2017;58:258–65.
Melo M, Gaspar da RA, Batista R, Vinagre J, Martins MJ, Costa G, et al. TERT, BRAF, and NRAS in Primary Thyroid Cancer and Metastatic Disease. J Clin Endocrinol Metab. 2017;102:1898–907.
Sohn SY, Park WY, Shin HT, Bae JS, Ki CS, Oh YL, et al. Highly concordant key genetic alterations in primary tumors and matched distant metastases in differentiated thyroid Cancer. Thyroid. 2016;26:672–82.
Jeon MJ, Chun SM, Lee JY, Choi KW, Kim D, Kim TY, et al. Mutational profile of papillary thyroid microcarcinoma with extensive lymph node metastasis. Endocrine. 2019;64:130–8.
Song E, Song DE, Ahn J, Kim TY, Kim WB, Shong YK, et al. Genetic profile of advanced thyroid cancers in relation to distant metastasis. Endocr Relat Cancer. 2020;27:285–93.
Fakhruddin N, Jabbour M, Novy M, Tamim H, Bahmad H, Farhat F, et al. BRAF and NRAS mutations in papillary thyroid carcinoma and concordance in BRAF mutations between primary and corresponding lymph node metastases. Sci Rep. 2017;7:4666.
Saiselet M, Gacquer D, Spinette A, Craciun L, Caussin-Petrucci M, Andry G, et al. New global analysis of the microRNA transcriptome of primary tumors and lymph node metastases of papillary thyroid cancer. BMC Genomics. 2015;16:828.
Romeo P, Colombo C, Granata R, Calareso G, Gualeni AV, Dugo M, et al. Circulating miR-375 as a novel prognostic marker for metastatic medullary thyroid cancer patients. Endocr Relat Cancer. 2018;25:217–31.
Yoo SK, Lee S, Kim SJ, Jee HG, Kim BA, Cho H, et al. Comprehensive analysis of the transcriptional and mutational landscape of follicular and papillary thyroid cancers. PLoS Genet. 2016;12:e1006239.
Na KJ, Choi H. Immune landscape of papillary thyroid cancer and immunotherapeutic implications. Endocr Relat Cancer. 2018;25:523–31.
Kim K, Jeon S, Kim TM, Jung CK. Immune gene signature delineates a subclass of papillary thyroid Cancer with unfavorable clinical outcomes. Cancers (Basel). 2018;10:494.
Fugazzola L, Cerutti N, Mannavola D, Vannucchi G, Beck-Peccoz P. The role of pendrin in iodide regulation. Exp Clin Endocrinol Diabetes. 2001;109:18–22.
Chakravarty D, Santos E, Ryder M, Knauf JA, Liao XH, West BL, et al. Small-molecule MAPK inhibitors restore radioiodine incorporation in mouse thyroid cancers with conditional BRAF activation. J Clin Invest. 2011;121:4700–11.
The authors thank from the Platform of Integrated Biology, DRAST, Dr. Marialuisa Sensi for her valuable advices and for critical reading of the manuscript and Mr. Edoardo Marchesi for performing gene/miRNA expression wet analysis.
This work was supported by Ministero della Salute (Grant RF-2013-02354985 to M.G. Borrello) and by Fondazione Umberto Veronesi (Post-Doctoral Fellowship granted to E. Minna).
Ethics approval and consent to participate
This was a retrospective cohort study approved by the Ethical Committee of IRCCS Istituto Auxologico Italiano (Institutional Review Board approval #2018_09_25_04). Patient informed consent was obtained for the use of thyroid tumor tissues and collection of clinicopathologic information.
Consent for publication
All authors give consent for the publication of the manuscript.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Tumor samples driving lesion by PTC-MA and BRAF-/RAS-like subtype. Suppl. Figure S2. Gene expression in RAI-refractory and RAI-avid primary and LN metastatic PTCs. Suppl. Figure S3. Expression of deregulated miRNAs in our series of RAI-refractory PTCs. Suppl. Figure S4. Differentially expressed genes in our series of RAI-refractory PTCs. Suppl. Figure S5. Expression of selected genes according to the most deregulated pathways identified in our series of RAI-refractory PTCs. Suppl. Figure S6. Thyroid differentiation (TD) genes expression.
About this article
Cite this article
Colombo, C., Minna, E., Gargiuli, C. et al. The molecular and gene/miRNA expression profiles of radioiodine resistant papillary thyroid cancer. J Exp Clin Cancer Res 39, 245 (2020). https://doi.org/10.1186/s13046-020-01757-x