Purpose:

The determinants of response or resistance to radioiodine (RAI) are unknown. We aimed to identify genomic and transcriptomic factors associated with structural responses to RAI treatment of metastatic thyroid cancer, which occur infrequently, and to test whether high MAPK pathway output was associated with RAI refractoriness.

Experimental Design:

Exceptional response to RAI was defined as reduction of tumor volume based on RECIST v1.1. We performed a retrospective case–control study of genomic and transcriptomic characteristics of exceptional responders (ER; n = 8) versus nonresponders (NR; n = 16) matched by histologic type and stage at presentation on a 1:2 ratio.

Results:

ER are enriched for mutations that activate MAPK through RAF dimerization (RAS, class 2 BRAF, RTK fusions), whereas NR are associated with BRAFV600E, which signals as a monomer and is unresponsive to negative feedback. ER have a lower MAPK transcriptional output and a higher thyroid differentiation score (TDS) than NR (P < 0.05). NR are enriched for 1q-gain (P < 0.05) and mutations of genes regulating mRNA splicing and the PI3K pathway. BRAFV600E tumors with 1q-gain have a lower TDS than BRAFV600E/1q-quiet tumors and transcriptomic signatures associated with metastatic propensity.

Conclusions:

ER tumors have a lower MAPK output and higher TDS than NR, whereas NR have a high frequency of BRAFV600E and 1q-gain. Molecular profiling of thyroid cancers and further functional validation of the key findings discriminating ER from NR may help predict response to RAI therapy.

Translational Relevance

Radioactive iodine (RAI) has been traditionally used in the post-operative management of patients with thyroid cancer. In this adjuvant setting, it is not possible to determine whether responses are due to the surgical resection itself or to the effects of radioiodine. To ascertain responsiveness to RAI more rigorously, we studied patients with metastatic thyroid carcinoma who had an objective structural response to this therapy. We evaluated genomic and transcriptomic characteristics of exceptional responders (ER) and nonresponders (NR) to RAI with a case–control design. ER were enriched for mutations that activate MAPK through RAF dimerization (RAS, class 2 BRAF, RTK fusions), whereas NR were associated with BRAFV600E, which signals as a monomer, 1q-gain, and mutations of genes regulating mRNA splicing and the PI3K pathway. Molecular profiling of ER and NR, coupled with functional validation of key findings discriminating ER from NR, may help predict response to RAI therapy.

Radioactive iodine (RAI) therapy has been used to treat patients with distant metastases of differentiated thyroid cancer since the 1940s (1). Early on, it was apparent that patients showing clinical responses were quite infrequent (2). Despite this, the indications for RAI were broadened to encompass adjuvant treatment of patients with differentiated thyroid cancer (DTC) following total thyroidectomy, which was widely adopted worldwide based primarily on large retrospective cohort studies (3). A recent phase III trial of patients randomized to RAI versus no RAI following total thyroidectomy for low-risk DTC (T1, N0, or Nx) showed no benefit of RAI therapy (4), consistent with prior retrospective studies (5, 6). Current consensus guidelines recommend consideration of adjuvant RAI therapy for patients following surgery for intermediate and high-risk DTCs based on prospective observational or population studies (5, 7). However, its use is undergoing a well-justified reappraisal because determining responses to RAI in the adjuvant setting is challenging, and because of the absence of high-quality evidence for benefit (6). There is also a greater appreciation that the toxicities of RAI, although relatively modest, cannot be disregarded (8, 9). Hence, there is a need to identify patients who are most likely to benefit from these treatments, which often require high administered activities of 131I.

Papillary thyroid cancers (PTC) harbor mutually exclusive activating mutations in BRAF (60%), RAS genes (12%), and gene fusions of receptor tyrosine kinases (RET, NTRK3, and others; 12%; refs. 10, 11). They all constitutively activate the MAP kinase signaling pathway, which is known to silence the expression of genes that govern iodide uptake and organification in thyrocytes (12, 13) in a manner that is inversely proportional to the flux through the pathway (14). The Cancer Genome Atlas study of PTC revealed that expression of genes required for iodine transport and metabolism is profoundly reduced in most, but not all, BRAFV600E-mutant cancers, whereas they are comparatively preserved in RAS-mutant PTCs. BRAFp.V600E signals as a monomer and is unresponsive to negative feedback by ERK, resulting in a higher MAPK signaling output (15). In contrast, mutant RAS proteins constitutively signal through RAF dimers, which are subject to negative feedback by ERK, resulting in an attenuation of the flux through the MAPK pathway. Nonetheless, there is considerable heterogeneity in these genotype-based predictions, particularly in the setting of advanced disease in older patients where the mutation burden is higher.

Molecular analyses of baseline and on-treatment lesional biopsies in patients on a redifferentiation trial with vemurafenib showed promising evidence that key transcriptomic signatures can segregate responders from nonresponders to RAI after MAPK inhibition (16). In addition, certain genomic lesions (e.g., mutations of genes encoding subunits of the SWI/SNF chromatin remodeling complex) lock thyroid cells into a dedifferentiated state that is unresponsive to MAPK pathway blockade (17). Integrating this information with that of patients who responded in an exceptional way to RAI should provide biomarkers that may be useful in clinical practice.

Studies of exceptional responders to targeted therapies have provided key insights into specific vulnerabilities of cancer cells. Salient examples include the discovery that TSC1 and TSC2 mutations confer sensitivity to rapalogs in bladder cancer (18), or the synthetic lethality to combined CHK1 inhibition and DNA-damaging chemotherapy in a patient with small cell carcinoma with ATM deficiency and a RAD50 hemizygous mutation (19). These and other similar reports are appealing because they are supported by a strong mechanistic rationale. There is already a detailed understanding of the transcriptional program required for iodide incorporation and organification, how it is affected by the oncogenic drivers of the disease and the signaling pathways they activate, and its potential reversibility, and yet this information is not used in clinical practice. To explore these questions, we have performed a case–control study of the genomic and transcriptomic characteristics of patients with exceptional responses to conventional RAI as determined by RECIST v1.1 and compared them with a cohort of RAI nonresponders with the goal of identifying biomarkers or mechanisms that predict for response to radioactive iodine.

Study population

We defined exceptional responders to radioactive iodine (ER) as patients whose index lesional volumes decreased after administration of conventional 131I therapy as determined by RECIST v1.1 criteria. We identified these patients by a data-line search of the Memorial Sloan Kettering Cancer Center (MSKCC) electronic medical record (EMR) as well as by individual referrals from colleagues. Among 20,367 unique patients with thyroid cancer treated at MSKCC between 1985 and 2018, 5,863 received RAI therapy within the institution. Of these, 3,363 had undetectable thyroglobulin antibodies (TgAb), rendering serum thyroglobulin (Tg) levels interpretable as a biomarker. Thyroidectomized patients with baseline Tg levels >200 ng/mL usually have distant metastatic disease (20). Altogether, 221 of 3,363 patients had Tg levels >200 ng/mL. Of these, 34 ended up with a Tg <20 ng/mL after RAI, but only 11 of 221 (5%) had a structural response by RECIST v1.1. We consented and collected surgical specimens from 8 ER patients (5 of them identified by the data-line search, and 3 by physician referrals) as well as surgical or biopsy samples from 16 patients with metastatic disease unresponsive to RAI, which were used as comparators in a 1:2 ratio. A nonresponder to RAI (NR) was defined as a patient not showing structural reduction of their tumor or showing structural progression after RAI therapy. ER and NR patients were matched by histopathologic tumor type and by disease stage at presentation blinded to genotype. Supplementary Table S1 details the clinical characteristics of ER to radioactive iodine. No patients with Hurthle cell carcinoma were included in this study.

Samples for molecular profiling

To comprehensively characterize the somatic single-nucleotide variants (SNV) and copy-number alterations in the patient cohort, we performed whole-exome sequencing (WES) of tumor/normal pairs (6/8 ER and 13/16 NR) as well as targeted cancer-exome sequencing using the MSK-IMPACT platform to identify SNV or fusions on a subset of the samples (2 ER, 15 NR; Fig. 2A). One ER sample (ERAI #3, Struma ovarii) failed QC for WES and MSK-IMPACT sequencing, and one sample each from ER and NR group were genotyped previously by Foundation One and Sequenom, respectively. We also performed RNA-seq on five samples each from ER and NR tumors. Altogether 4 to 10 slides per frozen tissue sample with pathologically confirmed tumor cells were scraped and collected in an Eppendorf or AutoLys M tube for nucleic acids extraction for WES and RNA-seq.

WES

FFPE tissue was deparaffinized, digested with a protease solution, and the DNA extracted using DNeasy Blood & Tissue Kit (Qiagen, Catalog No. 69504) or MagMAX FFPE DNA/RNA Ultra Kit (Thermo Fisher Scientific, Catalog No. A31881) on the KingFisher Flex Purification System (Thermo Fisher Scientific). DNA from whole blood or normal thyroid tissue were used as paired normal to filter out germline variants. After PicoGreen quantification and quality control by Agilent BioAnalyzer, 100 to 500 ng of DNA were used to prepare libraries using the KAPA Hyper Prep Kit (Kapa Biosystems KK8504) with eight cycles of PCR. After sample barcoding, 100 to 500 ng of library were captured by hybridization using the xGen Exome Research Panel v1.0 or v2.0 (IDT) according to the manufacturer's protocol. PCR amplification of the post-capture libraries was carried out for 12 cycles. Samples were run on a HiSeq 4000 or NovaSeq 6000 in a PE100 run, using the HiSeq 3000/4000 SBS Kit or NovaSeq 6000 S4 Reagent Kit (200 cycles; Illumina), respectively. Normal and tumor samples were covered to an average of 97X and 183X, respectively. The resulting FASTQ files were aligned and processed using TEMPO (code for processing WES data can be found at https://github.com/mskcc/tempo). Briefly, reads were aligned using Burroughs-Wheeler Aligner (BWA)-MEM (21) to the GRCh37 reference genome. Genome Analysis Toolkit (GATK, RRID:SCR_001876; ref. 22) best practices were used for base recalibration.

Somatic genome variants were called using the union of Mutect2 from GATK version 4.1.0.0 (23) and Strelka2 version 2.9.10 (24). Variants were then filtered out as follows: variant allele frequency <0.05, tumor read depth <20, tumor alternate read count <3, and normal read depth of <10, along with filtering out repeated regions from RepeatMasker (RRID:SCR_012954; ref. 25) and variants that appear at allele frequencies >0.01 in GNOMAD (26). TMB was calculated by TEMPO as the number of nonsynonymous mutations in coding regions covered by the baits. Variants were then annotated using MutSigCV v1.3.5 (https://www.genepattern.org/) and mutations identified as nonsilent or null were used for further analysis.

MSK Integrated Mutation Profiling of Actionable Cancer Target (IMPACT) sequencing

IMPACT sequencing is designed to capture all protein-coding exons and select introns of commonly implicated oncogenes, tumor suppressor genes, and members of pathways deemed actionable by targeted therapies (27). All samples received for MSK-IMPACT testing were accessioned by the CLIA-compliant Molecular Diagnostics Service laboratory. Tumor and patient-matched normal blood samples were matched for simultaneous sequencing using a patient-unique medical record number (MRN). A hematoxylin and eosin (H&E) stained slide was reviewed by a molecular pathologist and annotated for relevant specimen information including tumor type, tumor purity, and whether microdissection of the indicated tumor region was necessary prior to nucleic acid extraction. DNA samples were normalized to yield 50 to 250 ng input and diluted in 55 μL on the Biomek FXP Laboratory Automation Workstation (Beckman Coulter), prior to shearing on the Covaris instrument. Sequence libraries were prepared on the Biomek FXP through a series of enzymatic steps including shearing, end-repair, A-base addition, ligation of barcoded sequence adaptors, and low-cycle PCR amplification (Kapa Biosystems). Tumors and matched normal were combined in pools of 24 to 36 libraries for multiplexed captures using custom-designed biotinylated probes (Nimblegen). Captured DNA fragments were sequenced on an Illumina HiSeq2500 as paired-end 100-base pair reads. Tumor samples were sequenced to a mean unique depth of coverage of 718X.

The SNV found by WES were filtered against the list of 505 genes present in the MSK IMPACT platform. Of these, 58 were found to be mutated in the WES analysis. These, along with the ALK or RET fusions and TERT promoter mutations detected by MSK-IMPACT, were represented in the Oncoprint. The mutational allele frequencies were corrected using the BRAF or RAS driver mutation allele frequency, which were considered to be heterozygous. Genes with corrected allele frequency >10% were included in the Oncoprint. The somatic SNV identified by WES, the MSK-IMPACT list of 505 genes and the genes used in the Oncoprint, are provided in Supplementary Table S2. The chromosomal and cytogenetic locations of genes present in the MSK-IMPACT 505 gene panel are detailed in Supplementary Table S3.

Copy-number analysis

We used the FACETS algorithm to segment the chromosomes into copy-number states from the WES dataset and assign the total and minor copy numbers as well as cancer cell fraction for each segment (28). Each segment was classified as gained or lost according to whether the assigned total copy number exceeded 2 or was lower than 2. The gain or loss was called clonal if the segment cancer cell fraction was at least 75% as large as the sample purity and was called subclonal if it was lower. The resulting gain-loss assignments were shown as a heatmap. For each sample we also computed the percent genome altered: that is, megabases of gain or loss divided by genome length, and compared this between the two cohorts.

RNA sequencing (RNA-seq)

RNA-seq was performed on the subset of samples (five ER and five NR) that had frozen tissue in our biobank. RNA was extracted using the MagMAX FFPE DNA/RNA Ultra Kit (Thermo Fisher Scientific, Catalog No. A31881) on the KingFisher Flex Magnetic Particle Processor (Thermo Fisher Scientific, Catalog No. 5400630) or using RNeasy FFPE Kit (Qiagen, Catalog No. 73504) according to the manufacturer's instructions. After RiboGreen quantification and quality control by Agilent BioAnalyzer, 0.5 to 1 μg of total RNA with DV200 percentages varying from 25 to 60 underwent ribosomal depletion and library preparation using the TruSeq Stranded Total RNA LT Kit (Illumina, Catalog No. RS-122–1202) according to instructions provided by the manufacturer with eight cycles of PCR. Samples were barcoded and run on a NovaSeq 6000 in a PE50 run, using the NovaSeq 6000 S2 Reagent Kit (100 cycles).

For the analysis of gene expression from the TCGA study of PTC, the raw sequencing data (FASTQ files) were downloaded from GDC (https://docs.gdc.cancer.gov/). RNA-seq sequencing reads were aligned against human genome assembly hg19 by STAR 2-pass alignment (29). QC metrics, such as general sequencing statistics, gene feature, and body coverage were then calculated based on the alignment using RSeQC (30). RNA-seq gene-level counts were computed using the R package GenomicAlignments (31) over aligned reads with UCSC KnownGene (32) in hg19 as the base gene model. Union counting mode was used and only mapped paired reads after alignment quality filtering were considered.

Gene-level raw read count values, fragments per kilobase million (FPKM) counts and log 2 expression values were processed with the R package DESeq2 (RRID:SCR_000154; ref. 33). DESeq2 normalized the raw count data by sample-specific size factor and took sex as a covariate while testing for significant differences in gene expression between sample groups with multiple test correction through FDR. The normalized gene counts were used to compute the enhanced thyroid differentiation score (eTDS), which consists of the 16 TDS genes and the top 29 genes positively and 19 genes negatively associated with the TDS in the TCGA cohort, as well as the BRAF-RAS score as described (10, 16, 34). Log2 fold change gene expression values from this study and from the PTC TCGA (10) were used as input for GSEA analysis (35) against Hallmark, GO:BP (C5) and oncogenic signature (C6) gene sets from MSigDB (36) using the R package clusterProfiler (37).

Statistical analysis

Clinical variables were compared with Student t test and Mann–Whitney test for continuous variables and with Chi-squared test or Fisher exact test for categorical variables. Kaplan–Meier curves were generated for survival analyses and the log-rank test was used to compare the curves. A two-tailed P-value less than 0.05 was considered statistically significant for all calculations. STATA 17 or GraphPad-Prism (version 9.0; GraphPad Software, Inc. RRID:SCR_002798) were used for these calculations. The PCA analysis was performed through the R package DESeq2 (33). To assess the variance between the transcriptional landscapes, we summarized the data points into the first two principal components and evaluated the association between them and the covariates of interest. The association between the covariate and the principal components was evaluated using the proportion of variation explained (which is the ratio of the “between sum of squares” to the “total sum of squares” from the two components). We calculated the explained variation for all possible permutations of the covariate through enumeration and ranked the explained variation for the observed covariate pattern among them.

Study approval

Written informed consent was obtained from all patients and the study was conducted in accordance with the Declaration of Helsinki ethical guidelines. The study was reviewed and approved by the Institutional Review Board at MSKCC (IRB #16–016).

Data availability

Data were generated by the authors and available on request.

Cohort selection, clinical data, and outcomes with RAI therapy in patients with metastatic thyroid cancer

The clinical characteristics of ER and the corresponding NR are outlined in Table 1. The two groups were well matched in terms of tumor histopathology and M1 disease at presentation. They did not differ significantly in their median age and sex distribution, although women were overrepresented in the ER cohort. Most patients (87.5%) presented with distant metastases, but neither the peak suppressed thyroglobulin, nor the cumulative dose of RAI administered differed between the groups (Table 1). The nadir thyroglobulin levels, representing the lowest registered value documented after RAI therapy, were lower in ER versus NR (0.8 ng/mL vs. 187 ng/mL, P = 0.005), and the time to disease-progression was significantly longer in ER than that of NR (Table 1; Fig. 1D). The median (interquartile range) follow-up time for ER was 6.3 years (3.6–8.4) and for NR it was 4.6 years (2.8–8.6). Representative CT scans of ER patients before and after RAI therapy documenting the structural responses are shown in Fig. 1 and Supplementary Table S1).

Table 1.

Demographic and clinical characteristics of exceptional responders and nonresponders to RAI.

Exceptional responders (n = 8)Nonresponders (n = 16)P value
Median age (IQR) (years) 58.5 (33.5–65) 60.5 (48–65) 0.82 
Sex (F:M) 6:2 6:10 0.193 
Histopathology 
 PTC, n (%) 4 (50) 8 (50) 
 FTC, n (%) 1 (12.5) 2 (12.5)  
 Struma ovarii, n (%) 2 (25) 4 (25)  
 PDTC, n (%) 1 (12.5) 2 (12.5)  
Distant metastases at presentation   
  Without DM at presentation  
  With DM at presentation 14  
Median peak suppressed TG (IQR) (ng/mL) 845 (46.5–1888) 74 (23.3–3078) 0.92 
Median nadir TG after RAI (IQR) (ng/mL) 0.8 (0.15–44.85) 187 (6–928) 0.005 
Median cumulative dose of RAI (IQR) (mCi) 214 (125–336) 174 (147–258) 0.83 
Median time to disease progression (months) 76 (39–100) 11 (6–36) 0.0004 
Died, n (%) 1 (12.5) 5 (31) 0.621 
Exceptional responders (n = 8)Nonresponders (n = 16)P value
Median age (IQR) (years) 58.5 (33.5–65) 60.5 (48–65) 0.82 
Sex (F:M) 6:2 6:10 0.193 
Histopathology 
 PTC, n (%) 4 (50) 8 (50) 
 FTC, n (%) 1 (12.5) 2 (12.5)  
 Struma ovarii, n (%) 2 (25) 4 (25)  
 PDTC, n (%) 1 (12.5) 2 (12.5)  
Distant metastases at presentation   
  Without DM at presentation  
  With DM at presentation 14  
Median peak suppressed TG (IQR) (ng/mL) 845 (46.5–1888) 74 (23.3–3078) 0.92 
Median nadir TG after RAI (IQR) (ng/mL) 0.8 (0.15–44.85) 187 (6–928) 0.005 
Median cumulative dose of RAI (IQR) (mCi) 214 (125–336) 174 (147–258) 0.83 
Median time to disease progression (months) 76 (39–100) 11 (6–36) 0.0004 
Died, n (%) 1 (12.5) 5 (31) 0.621 
Figure 1.

Radiologic findings and treatment outcomes in exceptional responders and nonresponders to RAI therapy. A, Representative dorsal and ventral 131I whole body scan images of exceptional responder ER#1 showing RAI uptake in the disseminated metastatic lesions. B, Left, chest CT scans of ER#1 at the onset of RAI therapy; red circles circumscribe index metastatic lesions. Right, posttreatment CT scans show structural tumor regression by RECIST v1.1 criteria. C, Left, representative 131I SPECT/CT scans documenting RAI uptake in the spleen (top) and liver (bottom) in exceptional responder ER#3. Right, abdominal MRI scans showing resolution of liver and spleen metastatic lesions post-RAI therapy that were noted at the onset of therapy (red circles). D and E, Kaplan–Meier survival estimates of progression-free survival (D) and disease-specific survival (E) in ER compared with NR.

Figure 1.

Radiologic findings and treatment outcomes in exceptional responders and nonresponders to RAI therapy. A, Representative dorsal and ventral 131I whole body scan images of exceptional responder ER#1 showing RAI uptake in the disseminated metastatic lesions. B, Left, chest CT scans of ER#1 at the onset of RAI therapy; red circles circumscribe index metastatic lesions. Right, posttreatment CT scans show structural tumor regression by RECIST v1.1 criteria. C, Left, representative 131I SPECT/CT scans documenting RAI uptake in the spleen (top) and liver (bottom) in exceptional responder ER#3. Right, abdominal MRI scans showing resolution of liver and spleen metastatic lesions post-RAI therapy that were noted at the onset of therapy (red circles). D and E, Kaplan–Meier survival estimates of progression-free survival (D) and disease-specific survival (E) in ER compared with NR.

Close modal

Somatic SNV, structural variants, and copy-number alterations

The patient selection into the NR and ER cohorts was genotype-agnostic. We identified driver mutations of genes encoding effectors that signal through the MAPK pathway in all but one of the samples, an ER tumor whose degraded DNA was not informative (Fig. 2A). None of the tumors from the 8 patients in the ER group harbored the class I BRAFV600E mutation (of note, the BRAFV600E protein in the ERAI #3 sample was ruled out by IHC), as compared with 7 of 16 NR patients (ER 0% and NR 44%, Fisher exact test P value = 0.052). The oncogenic drivers in the ER cohort included class II BRAF and RAS mutations, as well as RET and ALK fusions. Notably, all the ER tumor drivers constitutively activate the MAPK pathway through RAF dimerization, in contrast to the predominance of BRAFV600E in the NR cohort, which signals as a monomer. The NR tumors harboring class II BRAF or RAS mutations were enriched for mutations in genes encoding for receptor tyrosine kinases, effectors in the mTOR pathway, as well as RNA splicing regulators as compared with ER (P = 0.04). Most epigenetic defects in NR were either truncating or nonsense mutations, including those of the SWI/SNF complex, whereas in two ER cases we detected missense mutations of unknown significance (Supplementary Table S2). TERT promoter and EIF1AX mutations were equally distributed between both groups (P = 0.99 and 0.18, respectively). Moreover, the tumor mutation burden (TMB) was not associated with differential responses to RAI therapy. Pathway analysis on WES data using “Hallmarks” as well as “Canonical Pathways” gene sets showed no specific pathway enrichment in either of the cohorts. This is likely reflective of the fact that thyroid cancers are known to harbor a very low TMB and few driver events.

Figure 2.

Somatic SNV and copy-number alterations in ER and NR. A, Oncoprint of somatic SNV in ER (left) and NR (right). Color key for sample characteristics as well as genetic alterations is shown on the right. Middle panel shows frequencies of alterations in individual genes or pathways in the ER or NR cohorts. B, Heatmap of genome-wide copy-number variations in the ER and NR tumors from the whole exome and MSK-IMPACT sequencing data. The genome was divided into Mb bins. A bin was considered as a clonal gain (red) if the total copy number was ≥2 and the tumor purity was at least 75%, or subclonal (orange) for copy-number gain with tumor purity <75%. Clonal (blue) and subclonal (cyan) losses were called if total copy number was <2. The samples were grouped by response to RAI, and the chromosome boundaries and centromeres are drawn. C, Frequencies of chromosome 1q gain and 22q loss in the ER and NR cohort (*, P < 0.05). Abbreviation: DDR, DNA damage response.

Figure 2.

Somatic SNV and copy-number alterations in ER and NR. A, Oncoprint of somatic SNV in ER (left) and NR (right). Color key for sample characteristics as well as genetic alterations is shown on the right. Middle panel shows frequencies of alterations in individual genes or pathways in the ER or NR cohorts. B, Heatmap of genome-wide copy-number variations in the ER and NR tumors from the whole exome and MSK-IMPACT sequencing data. The genome was divided into Mb bins. A bin was considered as a clonal gain (red) if the total copy number was ≥2 and the tumor purity was at least 75%, or subclonal (orange) for copy-number gain with tumor purity <75%. Clonal (blue) and subclonal (cyan) losses were called if total copy number was <2. The samples were grouped by response to RAI, and the chromosome boundaries and centromeres are drawn. C, Frequencies of chromosome 1q gain and 22q loss in the ER and NR cohort (*, P < 0.05). Abbreviation: DDR, DNA damage response.

Close modal

WES identified widespread copy-number alterations in the NR group with frequent gains in chromosome 1q and 5, and losses in chromosomes 11, 16, and 22q (Fig. 2B; Supplementary Fig. S1), whereas the ER tumors were comparatively quiet. Despite this, there was no significant difference in the percent of the genome altered by CNA between ER and NR (P = 0.97, Wilcoxon-rank sum test). This points to the particular significance and likely functional relevance of the recurrent Ch 1q and 22q CNAs, because there is no widespread genomic instability in these differentiated cancers. The chromosome 1q gain was present exclusively in the NR group (0/6 ER and 8/13 NR, P-value 0.05; Fig. 2C). Chromosome 1q gain has been reported to be associated with distant metastases and poor prognosis in PTC (38).

Transcriptomic characterization of exceptional responders and nonresponders to RAI therapy

We next performed bulk transcriptomic profiling of the thyroid tumors to identify gene expression changes that associate with response to RAI therapy. Unsupervised principal component analysis (PCA) of the RNA-seq data showed that response to RAI therapy had the highest explained variation among all possible (126) permutations. Sex had the third largest among 120 possible permutations and histology the 56th largest among 210 permutations (Fig. 3A; Supplementary Fig. S2A and S2C–S2E). The variance between the individual ER samples was lower compared with NR tumors, which may be due to the constraints imposed by the preserved differentiated state of ER, as compared with the loss of lineage identity and the disparate de-differentiated states in the NR. To determine whether the differential response to RAI therapy relates to differences in their differentiation state, we computed the eTDS, which comprises 64 genes associated with thyroid differentiation (Supplementary Table S4), including genes that play a role in iodine uptake, organification, and thyroid hormone biosynthesis (16). As shown in Fig. 3B, the eTDS was higher in ER compared with NR. The canonical 16-gene TDS was also higher in ER compared with NR (P-value = 0.02). In the TCGA PTC study, the transcriptomes of BRAF-mutant PTC were distinct from their RAS-mutant counterparts. A transcriptional signature consisting of 71 genes, the BRAF-RAS score (BRS), quantifies the extent to which a particular tumor transcriptionally resembles either a BRAFV600E or a RAS-mutant thyroid cancer. The TCGA found that BRAFV600E-like tumors had a higher ERK output score and lower TDS values compared with RAS-like tumors (10). We found that ER tumors were RAS-like, whereas NR were BRAFV600E-like (Fig. 3C). Further, there was a strong direct correlation between BRS and eTDS (Fig. 3D).

Figure 3.

Transcriptomic profiling and signaling pathway correlations between ER and NR. A, Variance in the transcriptional landscapes among individual tumor samples resolved into the first two principal components is shown as a PCA plot. B, eTDS in ER vs. NR computed from RNA-seq data using a set of genes that correlate with thyroid differentiation state. C, Classification of ER and NR tumors as BRAFV600E-like (negative y-axis values) or RAS-like (positive values) using the BRAF-RAS score (BRS). D, Correlation between TDS and BRS in ER and NR (*, P < 0.05; **, P < 0.01).

Figure 3.

Transcriptomic profiling and signaling pathway correlations between ER and NR. A, Variance in the transcriptional landscapes among individual tumor samples resolved into the first two principal components is shown as a PCA plot. B, eTDS in ER vs. NR computed from RNA-seq data using a set of genes that correlate with thyroid differentiation state. C, Classification of ER and NR tumors as BRAFV600E-like (negative y-axis values) or RAS-like (positive values) using the BRAF-RAS score (BRS). D, Correlation between TDS and BRS in ER and NR (*, P < 0.05; **, P < 0.01).

Close modal

Although sex distribution was similar between patients in the ER and NR group (Table 1), there was a sex imbalance in the samples available for transcriptional profiling by RNA-seq (ER: five females, NR: two females and three males; Supplementary Fig. S2B). We therefore performed DESEQ analysis using sex as a covariate to obtain sex-corrected log2 fold change expression values. The eTDS score computed from sex-corrected gene expression values was similar to the score from crude FPKM expression, and both confirmed that ER had ∼4.5-fold higher thyroid differentiated gene expression (Supplementary Fig. S3A). Interestingly, in the TCGA cohort sex did not appear to influence the BRS, whereas in RAS mutant tumors females had a higher TDS score and lower ERK score compared with their male counterparts (Supplementary Fig. S3B–S3D).

To investigate the relationship between thyroid differentiation and chromosome 1q gain, we analyzed the transcriptomic data from the TCGA thyroid cancer cohort. BRAFV600E-1q gain PTC had a lower TDS than BRAFV600E-1q-quiet PTC but did not have a higher ERK transcriptional output (Fig. 4AC). None of the RAS-mutant tumors from the TCGA study had chr 1q gain, consistent with the data in our patient tumor series. To identify genes that are overexpressed in chr1q that may play a role in thyroid differentiation, we mined the RNA-seq data from the TCGA study and compared expression of genes in chr1q in BRAFV600E tumors with chr1q gain with copy-number quiet BRAFV600E tumors. The gene expression changes between these two groups is shown in Supplementary Fig. S4 as a heatmap. Among these, only 15 genes had >1.5-fold expression in tumors with 1q gain (Supplementary Table S5), which suggests that transcriptomic changes associated with chr1q gain and its impact on thyroid differentiation is likely associated with effects on broader genome-wide transcription.

Figure 4.

Impact of chromosome 1q gain on thyroid differentiation and gene expression. A–C, TDS (A), ERK output score (B), and BRAF-RAS score (C) in BRAF or RAS mutant PTCs with chromosome 1q gain, 22q loss, or CN quiet in the TCGA PTC cohort (* P < 0.05; **** P < 0.0001). D and E, GSEA comparing: 1. ER vs. NR using crude (brown) and sex-corrected (green) differentially expressed genes from this study; 2. RAS vs. BRAF PTC from TCGA (orange); 3. BRAF PTC with 1q quiet vs. 1q gain from TCGA (blue) using the Gene Ontology-Biological Processes gene sets (D) and the Hallmark and Oncogenic Signature gene sets (E). The common gene set enrichments among four comparators with P < 0.05 are displayed as a heatmap using normalized enrichment scores. Gene set enrichment results for individual comparisons are shown in Supplementary Table S6.

Figure 4.

Impact of chromosome 1q gain on thyroid differentiation and gene expression. A–C, TDS (A), ERK output score (B), and BRAF-RAS score (C) in BRAF or RAS mutant PTCs with chromosome 1q gain, 22q loss, or CN quiet in the TCGA PTC cohort (* P < 0.05; **** P < 0.0001). D and E, GSEA comparing: 1. ER vs. NR using crude (brown) and sex-corrected (green) differentially expressed genes from this study; 2. RAS vs. BRAF PTC from TCGA (orange); 3. BRAF PTC with 1q quiet vs. 1q gain from TCGA (blue) using the Gene Ontology-Biological Processes gene sets (D) and the Hallmark and Oncogenic Signature gene sets (E). The common gene set enrichments among four comparators with P < 0.05 are displayed as a heatmap using normalized enrichment scores. Gene set enrichment results for individual comparisons are shown in Supplementary Table S6.

Close modal

We next performed gene set enrichment analysis (GSEA) to identify biological pathways associated with differential response to RAI therapy. There was an increase in thyroid differentiated expression gene sets in ER compared with NR, associated with a downregulation of developmentally regulated gene signatures, including an endoderm development gene set (Fig. 4D). Thus, although ER patients retained a more differentiated thyroid tumor consistent with their response to RAI therapy, NR tumors had loss of differentiated gene expression and likely recapitulated their early developmental/progenitor stage. There was excellent consistency between crude and sex-corrected ER/NR gene set enrichments. The developmental and differentiation-related gene set enrichments were preserved in RAS versus BRAFV600E tumors in the TCGA cohort, consistent with ER being RAS-like and NR being BRAF-like (Fig. 3C), as well as in the BRAFV600E-1q-quiet versus BRAFV600E-1q-gain tumors in the TCGA cohort. In addition to this, the GO:BP analysis showed that the BRAFV600E-1q-quiet TCGA tumors had higher activity of genetic programs of cell adhesiveness and wound healing (GO:BP_CELL_CELL_ADHESION_VIA_PLASMA_MEMBRANE_MOLECULES; GO:BP_REGULATION_OF_RESPONSE_TO_WOUNDING: and CELL_JUNTION_ASSEMBLY) as compared with BRAFV600E-1q gain tumors, consistent with the latter having a worse progression-free survival and a more advanced disease stage at presentation (10).

GSEA revealed that ER tumors had lower MAPK activation compared with NR (Fig. 4E). In addition, we found gene set enrichments indicative of downregulation of SWI/SNF function (SNF5_DN.V1_UP) and gain of polycomb repressive complex 2 PRC2 activity (PRC2_EZH2_UP.V1_DN and PRC2_SUZ12_UP.V1_DN) in NR. SWI/SNF chromatin remodeling complexes play a central role in maintenance of the differentiated state across lineages and are antagonistic to PRC2 activity (39, 40). SWI/SNF loss in tumors has been shown to upregulate stem-cell associated program via de-repression of PRC2 activity (40). In the context of thyroid, loss of SWI/SNF causes thyroid dedifferentiation and promotes refractoriness to MAPK inhibition- based redifferentiation therapies (17). These gene set enrichments were also present in RAS versus BRAF tumors and BRAF-1q quiet versus BRAF-1q-gain tumors, suggesting that repression of SWI/SNF function along with a concomitant gain in PRC2 activity is associated with loss of thyroid differentiation and refractoriness to RAI therapy. There was a paradoxical increase in PRC1 activity in ER compared to NR (BMI1_DN.V1_UP and MEL18_DN.V1_UP). Although the mechanisms that govern recruitment of PRC1 and PRC2 to their target sites remain unclear, their binding causes repression of common target genes and the binding sites overlap (41). However, contrary to the hierarchic and synergistic model of PRC1/2 action, it has been reported that these complexes may have opposing roles in HSC/progenitor cells (42).

RAI therapy is commonly given in the adjuvant setting following total thyroidectomy and lymph node dissection for thyroid cancers that are at an intermediate or high risk for recurrence based on the American Thyroid Association (ATA) risk stratification system (43). Clinical response is defined by the absence of structural recurrence as assessed primarily by thyroid ultrasonography, and by very low or undetectable serum thyroglobulin levels. In this context it is not possible to determine whether responses are due to the surgical resection itself or to the effects of radioiodine. Clinical recommendations for RAI adjuvant therapy are based on the best available evidence, which is derived from retrospective studies, with their inherent limitations (5, 7). The current guidelines of the ATA do not consider the evidence for tumor genotype as predictor of RAI responsiveness to be sufficiently validated in the adjuvant context (43).

We reasoned that RAI responsiveness could be ascertained more rigorously in patients with metastatic disease who had an objective structural response to this treatment. Although our EMR screen likely underrepresented the total number of patients with distant metastases, we were struck by the rarity of patients showing a structural response to RAI therapy and, as discussed below, the consistency of their genomic and transcriptomic characteristics.

We found that the driver mutations in the ER tumors activate MAPK in a manner that would predict for partial preservation of the pERK-driven negative feedback, because they all signal via RAF dimers, a key node subject to disruption by pERK upon activation of the pathway. This confers them with constitutively active, yet partially dampened, MAPK signaling output. In contrast, the NR were enriched for BRAFp.V600E, an oncoprotein that generates a higher MAPK pathway flux because it signals as a monomer and is therefore relatively unresponsive to negative feedback (15). The transcriptional readouts of the MAPK pathway output supported the distinction between ER, which were more RAS-like, as opposed to the NR, which were more BRAF-like. The TCGA study of PTCs found a very strong association between the BRS and TDS, a relationship confirmed in our series. Hence, it is likely that the lower output of the MAPK pathway was a key determinant of the preserved differentiation state of ER versus NR tumors.

Although attenuated MAPK pathway output may poise tumors to retain some degree of RAI-responsiveness, it is not sufficient, because 9 of 16 NR tumors harbored RAS or class II/III BRAF mutations. However, these NR tumors were associated with coexisting truncation mutations of tumor suppressor genes involved in epigenetic regulation, RNA splicing, and the DNA damage response pathway. Mutations of genes encoding for ARID1A and ARID1B, key components of the canonical BAF SWI/SNF chromatin remodeling complex, were identified in two of the NR tumors. SWI/SNF complex mutations are enriched in advanced metastatic PTC, PDTC, and ATC (44, 45). In the context of BRAFp.V600E, SWI/SNF complex subunit mutations result in loss of RAI uptake, and insensitivity to the redifferentiation effects of MAPK pathway inhibitors (17). In addition, mutations of genes involved in splicing regulation, particularly the cassette exon splicing regulator RBM10, were enriched in the NR cohort. Recent reports implicate mutations of the TERT promoter in conferring unresponsiveness to RAI (46). In our series, TERT promoter mutations were equally represented in the two groups.

Perhaps the most striking observation was the high prevalence of chromosome 1q gain in the NR. It was found in 8 of 13 (69%) of the NR and none of ER. Moreover, 5 of 7 of the RAS or BRAF class II mutant NR had chromosome 1q-gain. Our reanalysis of the TCGA data revealed that BRAF tumors with 1q-gain had a markedly decreased TDS compared with BRAF-1q-quiet tumors. Consistent with this, the GO:BP analysis showed that the program of thyroid hormone generation was suppressed in BRAF-1q-gain tumors. Although the two groups did not differ in their ERK transcriptional output score they did show significant differences in other MAPK output signatures in the GO:BP analysis. Regardless, it is likely that BRAF-1q-gain tumors impact thyroid differentiation through more complex, MAPK-independent mechanisms. The thyroid bud arises from the definitive endoderm and, upon lineage-specific transcription factor expression, undergoes directed differentiation and assembly into polarized thyroid follicular cells that participate in iodine uptake, oxidation, and incorporation into thyroid hormone (47). BRAF-1q-gain tumors showed downregulation of SWI/SNF function and induction of PRC2 activity, which induces a more compact chromatin state, consistent with the activation of early developmental signatures as well as those of the endodermal lineage, and loss of thyroid specification. Finally, BRAF-1q-gain tumors had decreased signatures of cell–cell adhesion and cell junction assembly and an increase in epithelial-to-mesenchymal transition, which may relate to the high frequency of chromosome 1q gain in advanced, metastatic thyroid cancers (38).

This study was based on a small cohort of patients showing partial or complete structural responses to RAI. This reflects the rarity of these cases, which is why we termed them exceptional responders. Despite that, the case–control design allowed us to derive novel insights into the biological determinants of RAI responsiveness, and of factors that confer refractoriness to this therapy. There was a uniform trend for RAI-responsive tumors to be driven by RAS-like MAPK oncoproteins that activate the pathway while preserving negative feedback inputs that moderate pathway output. BRAFp.V600E mutant metastatic tumors tend to be RAI refractory, but there is a small subset of BRAF-PTCs from the TCGA study that retain expression of thyroid differentiation markers (10). These, however, are small volume tumors that tend not to metastasize and that would not usually be considered for RAI therapy (48). The presence of chromosome 1q gain should be considered as an indicator of RAI-refractory disease, regardless of whether the primary driver is BRAFp.V600E or other oncoproteins in the MAPK pathway. Similarly, loss of function mutations of SWI/SNF subunits and of RBM10 appear to mark tumors that are unlikely to respond to RAI. Small molecule inhibitors of the MAPK pathway can restore RAI responsiveness in a subset of patients with RAI-refractory metastatic disease (16, 49–53). Analysis of responders and nonresponders to these therapies, coupled to functional experiments in mouse models, should help further delineate the molecular landscapes of RAI-avid and RAI-refractory metastatic disease.

N.D. Socci reports grants from NIH during the conduct of the study. J.A. Knauf reports other support from Kura Oncology outside the submitted work. In addition, J.A. Knauf has a patent for thyroid differentiation classifier for patient stratification in thyroid cancer issued; a patent for targeting the adaptive responses of BRAF or RAS-mutant thyroid cancers to RAF/MEK inhibitors with bispecific antibodies to TSHR and/or HER pending; and a patent for treatment of H-RAS-driven tumors issued, licensed, and with royalties paid from Kura Oncology. A.L. Ho reports grants from Bayer, Novartis, AstraZeneca, Verastem, Bristol Myers Squibb, Bioatla, Poseida Therapeutics, OncoC4, Genentech, Astellas, and Celldex; grants and personal fees from Kura Oncology, Merck, Elevar Therapeutics, Eisai, and Ayala; and personal fees from Physician Education Resources, Rgenta, Remix Therapeutics, Clinical Endocrinology Update, Prelude Therapeutics, Exelixis, Cellestia, Inxmed, and Affyimmune during the conduct of the study; in addition, A.L. Ho has a patent for lesional dosimetry methods for tailoring targeted radiotherapy issued to MSKCC. R.A. Ghossein reports grants from NIH during the conduct of the study. J.A. Fagin reports grants from NIH during the conduct of the study as well as grants from Eisai outside the submitted work; in addition, J.A. Fagin has a patent for anti-TSHR multi specific antibodies and uses thereof, application number: 6325769U, issued, and a patent for MSK ref. SK 2014-024-03, title: Treatment of H-Ras-Driven Tumors, application number: 15/305,788, issued. No disclosures were reported by the other authors.

L. Boucai: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, methodology, writing–original draft, writing–review and editing. M. Saqcena: Conceptualization, resources, data curation, software, formal analysis, supervision, investigation, methodology, writing–original draft, writing–review and editing. F. Kuo: Formal analysis, investigation, methodology, writing–review and editing. R.K. Grewal: Resources, software, investigation, visualization, writing–review and editing. N. Socci: Resources, data curation, software, formal analysis, validation, investigation, methodology, writing–review and editing. J.A. Knauf: Conceptualization, resources, data curation, investigation, methodology, writing–review and editing. G.P. Krishnamoorthy: Resources, data curation, investigation, methodology. M. Ryder: Resources, investigation, visualization, writing–review and editing. A.L. Ho: Conceptualization, resources, methodology, writing–review and editing. R.A. Ghossein: Resources, data curation, formal analysis, validation, writing–review and editing. L.G.T. Morris: Resources, data curation, formal analysis, methodology, writing–review and editing. V. Seshan: Software, formal analysis, validation, investigation, methodology, writing–review and editing. J.A. Fagin: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

We would like to thank Dr. Martin I. Surks for contributing a case of an exceptional responder to radioactive iodine. This work was supported by R01 CA50706, RO1 CA 249663, the Jayme and Peter Flowers fund, and the Frank D Cohen fund (to J.A. Fagin). The Byrne fund and Bite Me Cancer from the American Thyroid Association (to L. Boucai), and the NIH/NCI Cancer Center Support Grant P30 CA008748 (PI C. Thompson).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).

1.
Seidlin
SM
,
Marinelli
LD
,
Oshry
E
.
Radioactive iodine therapy; effect on functioning metastases of adenocarcinoma of the thyroid
.
J Am Med Assoc
1946
;
132
:
838
47
.
2.
Stanbury
J
,
Brownell
G
,
Barzelatto
J
,
Correa
J
,
Maisterrena
J
,
Cortazar
J
, et al
.
The metabolic fate of I131 after large therapeutic doses in patients with metastatic carcinoma of the thyroid
.
J Clin Endocrinol Metab
1952
;
12
:
1480
94
.
3.
Mazzaferri
EL
,
Young
RL
.
Papillary thyroid carcinoma: a 10 year follow-up report of the impact of therapy in 576 patients
.
Am J Med
1981
;
70
:
511
8
.
4.
Leboulleux
S
,
Bournaud
C
,
Chougnet
CN
,
Zerdoud
S
,
Al Ghuzlan
A
,
Catargi
B
, et al
.
Thyroidectomy without radioiodine in patients with low-risk thyroid cancer
.
N Engl J Med
2022
;
386
:
923
32
.
5.
Jonklaas
J
,
Sarlis
NJ
,
Litofsky
D
,
Ain
KB
,
Bigos
ST
,
Brierley
JD
, et al
.
Outcomes of patients with differentiated thyroid carcinoma following initial therapy
.
Thyroid
2006
;
16
:
1229
42
.
6.
Schvartz
C
,
Bonnetain
F
,
Dabakuyo
S
,
Gauthier
M
,
Cueff
A
,
Fieffé
S
, et al
.
Impact on overall survival of radioactive iodine in low-risk differentiated thyroid cancer patients
.
J Clin Endocrinol Metab
2012
;
97
:
1526
35
.
7.
Ruel
E
,
Thomas
S
,
Dinan
M
,
Perkins
JM
,
Roman
SA
,
Sosa
JA
.
Adjuvant radioactive iodine therapy is associated with improved survival for patients with intermediate-risk papillary thyroid cancer
.
J Clin Endocrinol Metab
2015
;
100
:
1529
36
.
8.
Hyer
S
,
Kong
A
,
Pratt
B
,
Harmer
C
.
Salivary gland toxicity after radioiodine therapy for thyroid cancer
.
Clin Oncol (R Coll Radiol)
2007
;
19
:
83
6
.
9.
Molenaar
RJ
,
Sidana
S
,
Radivoyevitch
T
,
Advani
AS
,
Gerds
AT
,
Carraway
HE
, et al
.
Risk of hematologic malignancies after radioiodine treatment of well-differentiated thyroid cancer
.
J Clin Oncol
2017
:
JCO2017750232
.
10.
Network
C
.
Integrated genomic characterization of papillary thyroid carcinoma
.
Cell
2014
;
159
:
676
90
.
11.
Kimura
ET
,
Nikiforova
MN
,
Zhu
Z
,
Knauf
JA
,
Nikiforov
YE
,
Fagin
JA
, et al
.
High prevalence of BRAF mutations in thyroid cancer: genetic evidence for constitutive activation of the RET/PTC-RAS-BRAF signaling pathway in papillary thyroid carcinoma
.
Cancer Res
2003
;
63
:
1454
7
.
12.
Mitsutake
N
,
Knauf
JA
,
Mitsutake
S
,
Mesa
C
,
Zhang
L
,
Fagin
JA
.
Conditional BRAFV600E expression induces DNA synthesis, apoptosis, dedifferentiation, and chromosomal instability in thyroid PCCL3 cells
.
Cancer Res
2005
;
65
:
2465
73
.
13.
Chakravarty
D
,
Santos
E
,
Ryder
M
,
Knauf
JA
,
Liao
X-H
,
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
.
14.
Nagarajah
J
,
Le
M
,
Knauf
JA
,
Ferrandino
G
,
Montero-Conde
C
,
Pillarsetty
N
, et al
.
Sustained ERK inhibition maximizes responses of BrafV600E thyroid cancers to radioiodine 1
.
J Clin Invest
2016
;
126
:
4119
24
.
15.
Lito
P
,
Pratilas
CA
,
Joseph
EW
,
Tadi
M
,
Halilovic
E
,
Zubrowski
M
, et al
.
Relief of profound feedback inhibition of mitogenic signaling by RAF inhibitors attenuates their activity in BRAFV600E melanomas
.
Cancer Cell
2012
;
22
:
668
82
.
16.
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
.
17.
Saqcena
M
,
Leandro-Garcia
LJ
,
Maag
JLV
,
Tchekmedyian
V
,
Krishnamoorthy
GP
,
Tamarapu
PP
, et al
.
SWI/SNF complex mutations promote thyroid tumor progression and insensitivity to redifferentiation therapies
.
Cancer Discov
2021
;
11
:
1158
75
.
18.
Iyer
G
,
Hanrahan
AJ
,
Milowsky
MI
,
Al-Ahmadie
H
,
Scott
SN
,
Janakiraman
M
, et al
.
Genome sequencing identifies a basis for everolimus sensitivity
.
Science
2012
;
338
:
221
.
19.
Al-Ahmadie
H
,
Iyer
G
,
Hohl
M
,
Asthana
S
,
Inagaki
A
,
Schultz
N
, et al
.
Synthetic lethality in ATM-deficient RAD50-mutant tumors underlies outlier response to cancer therapy
.
Cancer Discov
2014
;
4
:
1014
21
.
20.
Robbins
RJ
,
Srivastava
S
,
Shaha
A
,
Ghossein
R
,
Larson
SM
,
Fleisher
M
, et al
.
Factors influencing the basal and recombinant human thyrotropin-stimulated serum thyroglobulin in patients with metastatic thyroid carcinoma
.
J Clin Endocrinol Metab
2004
;
89
:
6010
6
.
21.
Li
H
.
Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM
,
2013
, pp
arXiv:1303.3997
22.
McKenna
A
,
Hanna
M
,
Banks
E
,
Sivachenko
A
,
Cibulskis
K
,
Kernytsky
A
, et al
.
The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data
.
Genome Res
2010
;
20
:
1297
303
.
23.
Benjamin
D
,
Sato
T
,
Cibulskis
K
, et al
.
Calling Somatic SNVs and Indels with Mutect2
.
Biorxiv
2019
:
861054
.
24.
Kim
S
,
Scheffler
K
,
Halpern
AL
,
Bekritsky
MA
,
Noh
E
,
Källberg
M
, et al
.
Strelka2: fast and accurate calling of germline and somatic variants
.
Nat Methods
2018
;
15
:
591
4
.
25.
Smit
AFA
,
Hubley
R
,
Green
P
.:
RepeatMasker Open-4.0, 2013 - 2015
.
26.
Karczewski
KJ
,
Francioli
LC
,
Tiao
G
,
Cummings
BB
,
Alföldi
J
,
Wang
Q
, et al
.
The mutational constraint spectrum quantified from variation in 141,456 humans
.
Nature
2020
;
581
:
434
43
.
27.
Zehir
A
,
Benayed
R
,
Shah
RH
,
Syed
A
,
Middha
S
,
Kim
HR
, et al
.
Mutational landscape of metastatic cancer revealed from prospective clinical sequencing of 10,000 patients
.
Nat Med
2017
;
23
:
703
13
.
28.
Arora
A
,
Shen
R
,
Seshan
VE
.:
FACETS: fraction and allele-specific copy number estimates from tumor sequencing
.
Methods Mol Biol
2022
;
2493
:
89
105
.
29.
Dobin
A
,
Davis
CA
,
Schlesinger
F
,
Drenkow
J
,
Zaleski
C
,
Jha
S
, et al
.
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
2013
;
29
:
15
21
.
30.
Wang
L
,
Wang
S
,
Li
W
.
RSeQC: quality control of RNA-seq experiments
.
Bioinformatics
2012
;
28
:
2184
5
.
31.
Lawrence
M
,
Huber
W
,
Pagès
H
,
Aboyoun
P
,
Carlson
M
,
Gentleman
R
, et al
.
Software for computing and annotating genomic ranges
.
PLoS Comput Biol
2013
;
9
:
e1003118
.
32.
Rosenbloom
KR
,
Armstrong
J
,
Barber
GP
,
Casper
J
,
Clawson
H
,
Diekhans
M
, et al
.
The UCSC Genome Browser database: 2015 update
.
Nucleic Acids Res
2015
;
43
:
D670
81
.
33.
Love
MI
,
Huber
W
,
Anders
S
.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
.
34.
Cancer Genome Atlas Research Network
.
Integrated genomic characterization of papillary thyroid carcinoma
.
Cell
2014
;
159
:
676
90
.
35.
Subramanian
A
,
Tamayo
P
,
Mootha
VK
,
Mukherjee
S
,
Ebert
BL
,
Gillette
MA
, et al
.
Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
.
Proc Natl Acad Sci U S A
2005
;
102
:
15545
50
.
36.
Liberzon
A
,
Birger
C
,
Thorvaldsdóttir
H
,
Ghandi
M
,
Mesirov
JP
,
Tamayo
P
.
The Molecular Signatures Database (MSigDB) hallmark gene set collection
.
Cell Syst
2015
;
1
:
417
25
.
37.
Yu
G
,
Wang
L-G
,
Han
Y
,
He
Q-Y
.
clusterProfiler: an R package for comparing biological themes among gene clusters
.
OMICS
2012
;
16
:
284
7
.
38.
Kjellman
P
,
Lagercrantz
S
,
Höög
A
,
Wallin
G
,
Larsson
C
,
Zedenius
J
.
Gain of 1q and loss of 9q21.3-q32 are associated with a less favorable prognosis in papillary thyroid carcinoma
.
Genes Chromosomes Cancer
2001
;
32
:
43
9
.
39.
Alver
BH
,
Kim
KH
,
Lu
P
,
Wang
X
,
Manchester
HE
,
Wang
W
, et al
.
The SWI/SNF chromatin remodelling complex is required for maintenance of lineage specific enhancers
.
Nat Commun
2017
;
8
:
14648
.
40.
Wilson
BG
,
Wang
X
,
Shen
X
,
McKenna
ES
,
Lemieux
ME
,
Cho
Y-J
, et al
.
Epigenetic antagonism between polycomb and SWI/SNF complexes during oncogenic transformation
.
Cancer Cell
2010
;
18
:
316
28
.
41.
Boyer
LA
,
Plath
K
,
Zeitlinger
J
,
Brambrink
T
,
Medeiros
LA
,
Lee
TI
, et al
.
Polycomb complexes repress developmental regulators in murine embryonic stem cells
.
Nature
2006
;
441
:
349
53
.
42.
Majewski
IJ
,
Ritchie
ME
,
Phipson
B
,
Corbin
J
,
Pakusch
M
,
Ebert
A
, et al
.
Opposing roles of polycomb repressive complexes in hematopoietic stem and progenitor cells
.
Blood
2010
;
116
:
731
9
.
43.
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
.
44.
Landa
I
,
Ibrahimpasic
T
,
Boucai
L
,
Sinha
R
,
Knauf
JA
,
Shah
RH
, et al
.
Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers
.
J Clin Invest
2016
;
126
:
1052
66
.
45.
Jin
M
,
Song
DE
,
Ahn
J
,
Song
E
,
Lee
Y-M
,
Sung
T-Y
, et al
.
Genetic profiles of aggressive variants of papillary thyroid carcinomas
.
Cancers (Basel)
2021
;
13
:
892
.
46.
Soe
MH
,
Chiang
JM
,
Flavell
RR
,
Khanafshar
E
,
Mendoza
L
,
Kang
H
, et al
.
Non-iodine-avid disease is highly prevalent in distant metastatic differentiated thyroid cancer with papillary histology
.
J Clin Endocrinol Metab
2022
;
107
:
e3206
16
.
47.
De Felice
M
,
Di Lauro
R
.
Thyroid development and its disorders: genetics and molecular mechanisms
.
Endocr Rev
2004
;
25
:
722
46
.
48.
Boucai
L
,
Seshan
V
,
Williams
M
,
Knauf
JA
,
Saqcena
M
,
Ghossein
RA
, et al
.
Characterization of subtypes of BRAF-mutant papillary thyroid cancer defined by their thyroid differentiation score
.
J Clin Endocrinol Metab
2022
;
107
:
1030
9
.
49.
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
.
50.
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
.
51.
Jaber
T
,
Waguespack
SG
,
Cabanillas
ME
,
Elbanan
M
,
Vu
T
,
Dadu
R
, et al
.
Targeted therapy in advanced thyroid cancer to resensitize tumors to radioactive iodine
.
J Clin Endocrinol Metab
2018
;
103
:
3698
-
705
.
52.
Iravani
A
,
Solomon
B
,
Pattison
DA
,
Jackson
P
,
Ravi Kumar
A
,
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
.
53.
Weber
M
,
Kersting
D
,
Riemann
B
,
Brandenburg
T
,
Führer-Sakel
D
,
Grünwald
F
, et al
.
Enhancing radioiodine incorporation into radioiodine-refractory thyroid cancer with MAPK inhibition (ERRITI): a single-center prospective two-arm study
.
Clin Cancer Res
2022
;
28
:
4194
202
.