Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers

A comprehensive investigation of the genomic landscape of papillary thyroid carcinomas (PTC), the most common thyroid malignancy, was recently reported by The Cancer Genome Atlas Network (TCGA Network) (1). These well-differentiated tumors were found to have a low frequency of somatic alterations (2), with the majority harboring mutually exclusive activating mutations in BRAF (60%) and RAS-family genes (13%), as well as fusion oncoproteins, primarily involving receptor tyrosine kinases (RTKs) such as RET, NTRK1 or -3, and ALK. Distinct signaling and transcriptomic consequences were observed between BRAF^V600E-like tumors, which showed higher MAPK transcriptional output and lower expression of genes involved in iodine metabolism, and RAS-like tumors, which had lower MAPK signaling and comparatively preserved expression of iodine-related genes.

The TCGA study excluded poorly differentiated thryoid cancers (PDTCs) and anaplastic thyroid cancers (ATCs) from their analysis in order to focus on a homogeneous histological cohort that would provide sufficient power to identify low-frequency genomic events. Although PDTCs and ATCs account for approximately 5%–10% of thyroid cancers, they represent a major clinical challenge. Patients with PDTC and ATC have a mean survival after diagnosis of 3.2 and 0.5 years, respectively, and account for approximately a third of deaths caused by this disease (3). Virtually all cases are refractory to radioiodine therapy, and traditional chemotherapy and radiotherapy are of marginal benefit (4, 5). Molecularly targeted approaches are being tested in preclinical studies and in early human clinical trials (6, 7). These efforts are constrained by the paucity of information on the genomics of these cancers, which have been investigated primarily through Sanger sequencing of a limited set of candidate genes (8–20). The exception to this is a recently reported whole exome sequencing study of ATC (21), which, although informative, may have undercalled significant mutations — particularly those that are subclonal — because of low tumor purity. This is because ATCs, and to a lesser extent PDTCs, pose a particular challenge for genomic studies due to their extensive infiltration by macrophages (22, 23). This cannot be overcome by microdissection of tissue samples because the macrophages form an interconnected network that envelops the individual tumor cells throughout the tumor specimen. To overcome this, we adopted an ultradeep sequencing strategy using the MSK-IMPACT cancer exome panel, a massively parallel exon capture approach that targets all exons and selected introns of 341 genes frequently altered in human cancer (24). We performed this extensive cancer gene exome sequencing as well as expression profiling in the largest series of PDTC and ATC ever investigated, and we identified a wide spectrum of somatic mutations, genetic fusions, and copy number alterations (CNAs) that clearly delineate profound genomic differences between the two advanced forms of the disease. Moreover, when analyzed in the context of the PTC TCGA study (1), this study provides insights into tumor microevolution, suggesting that PDTCs and ATCs evolve from their well-differentiated counterparts.

Samples, clinical data, and overall approach

One-hundred and seventeen advanced thyroid tumors, including 84 PDTCs and 33 ATCs, met the sequencing quality standards and are reported in this study. Clinicopathological features are summarized in Table 1 and Supplemental Table 1 (supplemental material available online with this article; doi:10.1172/JCI85271DS1). Median age was 58 and 66 years for PDTC and ATC, respectively. Female/male ratios were 1.5:1 (PDTC) and 1.2:1 (ATC), which are distinct from PTC (2.7:1) (25). Most samples were primary tumors: 92/117 (64/84 PDTC; 28/33 ATC), and the remainder from nodal (6) or distant (19) metastases.

Table 1

Summary of the clinicopathological features of the 117 advanced thyroid tumors included in the study

MSK-IMPACT, a next-generation sequencing (NGS) platform for targeted sequencing of 341 cancer genes (24), was performed in all 117 tumors. A subset of 37 (17 PDTCs and 20 ATCs) was also subjected to an mRNA expression array and a CGH platform to validate IMPACT-derived copy number calls, as detailed later. Median tumor purity was 72% and 42% for PDTCs and ATCs, respectively, consistent with the known heavy macrophage infiltration in ATCs and highlighting the need for deeper sequence coverage in these tumors (Figure 1B). Average depth of coverage was 584× for tumors and 236× for paired normal tissues. Average coverage was 500× and 765×, for FFPE and frozen tumors, respectively. Coverage for ATCs was 739×.

Figure 1

Cancer genome alterations in 84 PDTCs and 33 ATCs. (A) Mutation density across the PDTC and ATC cohorts (n = 117), expressed as number of genetic alterations found in 341 genes present in MSK-IMPACT. (B) Clinicopathological features, including sample type, patients’ age (by decade), sex, metastasis site, survival status, tumor purity, cytological phenotype, growth pattern, and PDTC definition (PDTC-Turin vs. PDTC-MSK). Color keys are shown in the right outermost panel. (C–G) Oncoprints of PDTCs (left) and ATCs (right). Middle panel shows percentage of tumors altered for each event; *P < 0.05 between PDTCs vs. ATCs using Fisher’s exact test (see Supplemental Table 3 for extended information). Color key for genetic alterations is shown in the bottom panel. Mutations in drivers and other relevant genes (C); fusion events (D); TERT promoter mutations (E); mutations in TP53 and other tumor suppressor genes (F); and alterations in key pathways and functional groups: PI3K/AKT pathway (includes PIK3CA, PTEN, PIK3C2G, PIK3CG, PIK3C3, PIK3R1, PIK3R2, AKT3, TSC1, TSC2, and MTOR), SWI/SNF chromatin remodeling complex (ARID1A, ARID1B, ARID2, ARID5B, SMARCB1, PBRM1, and ATRX), HMTs (KMT2A, KMT2C, KMT2D, and SETD2), and MMR (includes MSH2, MSH6, and MLH1 genes) (G). See Figure 4 for detailed mutational information.

Somatic mutations

ATCs harbored a higher number of mutations than PDTCs (median ± interquartile range [IQR]: 6 ± 5 and 2 ± 3, respectively, P < 1 × 10^–4) (Figure 1A). The mutation burden in PDTCs was increased compared with the PTCs from the TCGA study (considering mutations in the 341 gene set only): 2 ± 3 and 1 ± 1, respectively (P < 1 × 10^–4). These differences remained highly significant after removing tumors with defects in mismatch excision repair (MMR) genes, which showed a disproportionately higher number of mutations. Mutation burden in PDTCs (expressed as number of mutations below, equal, or above the median) was greater in older patients (47 vs. 58 vs. 64 years, P < 1 × 10^–3) and associated with tumor size (36% vs. 43% vs. 71% > 4 cm, P = 0.04), presence of distant metastasis (8% vs. 29% vs. 57%, P = 2 × 10^–3), and overall survival (19% vs. 25% vs. 46%, logrank P = 0.01) (Supplemental Table 2).

Drivers and frequently altered genes. BRAF^V600E mutations were present in 33% of PDTCs and 45% of ATCs, whereas mutations in NRAS, HRAS, or KRAS occurred in 28% and 24% of PDTCs and ATCs, respectively, and were mutually exclusive with BRAF and gene fusions (Figure 1, C and D, and Supplemental Tables 3 and 4). There was a dichotomy in the distribution of RAS and BRAF mutations of PDTCs according to their histological features. Ninety-two percent of RAS mutations were found in PDTCs fulfilling the Turin definition of the disease (PDTC-Turin; see Methods, ref. 26). By contrast, 81% of BRAF mutations were found in PDTCs defined based only on MSKCC criteria (PDTC-MSK): high mitotic rate and necrosis irrespective of growth pattern (ref. 27 and Figure 1, B and C). BRAF-mutant PDTCs were smaller and had higher frequency of nodal metastases, whereas their RAS-mutant counterparts were larger and had a higher rate of distant metastasis (Supplemental Table 5). BRAF-mutant PDTCs were significantly overrepresented among female patients (P = 0.005).

Neurofibromin 1 gene (NF1) mutations were identified in 3 BRAF/RAS WT ATCs. There were also low-frequency mutations in TSHR and STK11 in both PDTCs and ATCs (Figure 1C). Mutations in PIK3CA and PTEN, which encode key effectors of the PI3K/AKT pathway, were particularly prevalent in ATCs (18% and 15%, respectively, Supplemental Table 3) and were overrepresented with respect to PDTCs (P = 4 × 10^–2 and 6 × 10^–3, respectively). PIK3CA and PTEN showed distinct patterns of co-occurrence in ATCs. All 3 ATCs harboring NF1 mutations also had truncating alterations in PTEN (P = 2 × 10^–3), whereas PIK3CA and BRAF mutations tended to co-occur. All 5 PIK3CA helical domain mutations (E542K or E545K) occurred in ATCs, whereas the single kinase domain mutation (H1047R) was found in a PDTC (Figure 1C).

The EIF1AX-RAS association. Mutations in the eukaryotic translation initiation factor EIF1AX were initially discovered in uveal melanomas (28) and were also reported in 1% of PTCs (6/402), largely occurring in a mutually exclusive manner with BRAF and RAS (1). By contrast, 11% of PDTCs and 9% of ATCs harbored EIF1AX mutations (Figure 1C and Figure 2A), which were strongly associated with RAS (14/15, P < 1 × 10^–4, Figure 2B and Table 2). EIF1AX mutations clustered in two regions: the N-terminal domain, as also observed in uveal melanomas (28) and other tumors (Figure 2A), or at a unique splice acceptor site between exons 5 and 6 (p.A113splice), which was the most prevalent abnormality and which has not been reported in other tumor types. The C-terminal p.A113 splice mutation predicts for alternative usage of a cryptic splice acceptor within exon 6, resulting in a 12–amino acid in-frame deletion. Our analysis of RNASeq data from two cases with this mutation in the PTC-TCGA confirms this prediction (not shown). EIF1AX mutations were associated with larger tumors and predicted for shorter survival in PDTCs (Figure 2C and Supplemental Table 6).

Figure 2

EIF1AX mutations and EIF1AX-RAS co-occurrence in thyroid cancers. (A) Distribution of EIF1AX mutations in thyroid cancers and other tumors (modified from ref. 1). (B) Oncoprints showing the co-occurrence of EIF1AX with RAS mutations in PTC from TCGA (top, n = 401), PDTCs from our study (middle, n = 84), and ATCs from our series and from (21) (bottom, n = 55). (C) Kaplan-Meier graph showing significantly shorter survival in EIF1AX-mutated PDTCs (log-rank P = 0.048). See Supplemental Table 6 for detailed clinical correlations.

Table 2

Contingency analysis of EIF1AX-RAS mutations in advanced thyroid cancers (PDTCs and ATCs)

TERT in advanced thyroid cancer. There was a high prevalence of telomerase reverse transcriptase (TERT) promoter mutations in these advanced thyroid tumors, and their presence was associated with aggressive and metastatic phenotypes (Figure 1E). Together, 40% of PDTCs and 73% of ATCs harbored TERT promoter mutations (49/117 C228T [c.-124G>A]; 8/117 C250T [c.-146G>A]) as compared with 9% of PTCs from TCGA (Figure 3, A and B, and Table 3). Whereas TERT promoter mutations were subclonal in the small subset of PTC that harbored them, they were clonal in PDTC and ATCs (Figure 3C). TERT mutations co-occurred with BRAF/RAS mutations in PDTCs and ATCs combined (P = 4 × 10^–3, Figure 3B and Table 3), consistent with the proposed mechanism whereby the TERT mutations generate de novo binding elements for ETS-family transcription factors activated by MAPK signaling, such as GABPA (29). Survival of ATC patients harboring TERT promoter mutations was markedly diminished (732 vs. 147 days, P = 0.03, Supplemental Table 7), particularly in cancers with coexisting mutations of BRAF or RAS (Figure 3D). TERT-mutated PDTCs developed more distant metastases (56% vs. 20%, P = 0.01) and had a trend toward greater mortality (Supplemental Table 7). RAS-mutant ATCs with TERT (5/8) and EIF1AX mutations (3/8) did not overlap (odds ratio = 0.01; P = 0.02), consistent with alternate pathways toward progression to ATC.

Figure 3

TERT promoter mutations in thyroid cancers. (A) Location and overall frequency of TERT promoter mutations in PDTCs and ATCs. (B) Oncoprints of TERT promoter mutations vs. BRAF and RAS in (top) PTCs from TCGA (n = 381); (middle) PDTCs (n = 84); and (bottom) ATCs (n = 33). (C) Allelic frequency of TERT promoter mutations in thyroid cancers. Graph shows TERT mutant allelic frequency (MAF) corrected for tumor purity, determined based on allelic fraction of driver mutations (BRAF or RAS) for all three tumor types. (D) Kaplan-Meier survival in ATCs with log-rank P values. Top: TERT-mutant vs WT. Bottom: WT, TERT-mutant with or without BRAF/RAS mutations. See also Supplemental Table 7.

Table 3

Contingency analysis of TERT-BRAF/RAS mutations in thyroid cancers

TP53 and other tumor-suppressor genes. Inactivation of p53 has been considered as a hallmark of advanced thyroid tumors (14, 30). We found that TP53 mutations, although highly prevalent in ATCs, were relatively rare in PDTCs (73% vs. 8%, P < 1 × 10^–4, Figure 1F), which contrasts with previous reports (refs. 11, 31, and Supplemental Table 8) and constitutes a key distinguishing event in the biology of these tumors. ATM, a cell-cycle checkpoint and DNA damage response gene, was mutated in 7% of PDTCs and 9% of ATCs. ATM-mutated tumors had a higher mutation burden: median was 5 and 2 in ATM-mutant versus WT PDTCs (P = 0.04), and 19 and 5.5 in ATC, respectively (P = 7 × 10^–3). A higher mutation burden is consistent with the loss of the canonical function of this checkpoint kinase, which is activated in response to DNA double strand breaks and is required for appropriate DNA repair (32). Infrequent truncating mutations were also found in RB1, NF2, and MEN1.

Novel genes and pathways altered in advanced thyroid tumors. Mutations of genes encoding members of the PI3K/AKT/mTOR pathway were seen more frequently in ATCs than PDTCs (39% vs. 11%, P = 1 × 10^–3). Besides PIK3CA and PTEN, mutations of PIK3C2G, PIK3CG, PIK3C3, PIK3R1, PIK3R2, AKT3, TSC1, TSC2, and MTOR were also present (Figure 1G, Figure 4, A and E, and Supplemental Figure 1).

Figure 4

Pathways and novel functional groups mutated in advanced thyroid tumors. Expanded oncoprints of genes belonging to the indicated functional categories, as defined in Figure 1G. Samples are divided by tumor type (ATC or PDTC) within each panel. Only altered cases, out of 117 tumors, are shown. Missense, truncating, and in-frame mutations are represented as green, black, and brown squares, respectively. (A) PI3K/AKT/mTOR pathway (includes PIK3CA, PTEN, PIK3C2G, PIK3CG, PIK3C3, PIK3R1, PIK3R2, AKT3, TSC1, TSC2, and MTOR); (B) SWI/SNF chromatin remodeling complex (ARID1A, ARID1B, ARID2, ARID5B, SMARCB1, PBRM1, and ATRX); (C) HMTs (KMT2A, KMT2C, KMT2D, and SETD2); and (D) MMR (MSH2, MSH6, and MLH1). (E) Percentage of tumors altered for each functional category and tumor type.

Genes encoding components of the SWI/SNF chromatin remodeling complex were mutated in 36% of ATCs and 6% of PDTCs (P = 1 × 10^–4). This is the first report of mutations in ARID1A, ARID1B, ARID2, ARID5B, SMARCB1, PBRM1, and ATRX genes in advanced thyroid tumors. Consistent with evidence that disruption of one protein in this complex is typically sufficient to impair function (33), we observed a pattern of mutual exclusivity of mutations in this category of genes (Figure 4B). The single exception was for an ATC with concurrent missense mutations in ARID1A (minor allele frequency [MAF] = 5%), ARID1B (MAF = 6%), and a frameshift change in SMARCB1 (MAF = 22%), although their frequencies suggest that they are likely subclonal events within a heterogeneous tumor.

Mutations of the histone methyltransferases (HMTs) KMT2A, KMT2C, KMT2D, and SETD2 were found in 24% of ATCs and 7% of PDTCs (P = 0.02) (Figure 4C). Additional mutations in chromatin remodeling and epigenetic regulators other than SWI/SNF and HMTs were also seen, including frequent alterations affecting the histone acetyltransferase CREBBP and sporadic inactivating mutations in other epigenetic players such as EP300, BCOR, and BCL6 (Supplemental Figure 1). Mutations in a few epigenetic regulators (ARID1B, KMT2A, and KMT2C) were also identified in 1%–2% of PTCs from the TCGA, but these events are clearly enriched in advanced thyroid tumors.

Alterations in members of the DNA MMR pathway, including MSH2, MSH6, and MLH1, were found in 12% of ATCs and 2% of PDTCs (Figure 4D). MMR mutant tumors showed a “hypermutator phenotype” (as described in ATCs; ref. 21): median mutation number in MMR-mutant vs. WT ATCs was 16.5 and 5 (P = 1 × 10^–3) and, in PDTCs, 7.5 and 2 (P = 9 × 10^–3), respectively.

Other genes and functional categories were mutated in a small proportion of PDTCs and ATCs (Supplemental Figure 1 and Supplemental Table 4), i.e., RTKs such as EPHA3 (3 mutations, exclusively in ATCs), EGFR, FLT1 (VEGFR1), FLT4 (VEGFR3), and KDR (VEGFR2), as well as in all four members of the NOTCH family (NOTCH1–4). Finally, there were infrequent mutations (mutated in at least 2 ATCs or 3 PDTCs) in DIS3, FAT1, POLE, RBM10, RAD54L, RECQL4, and SF3B1.

Mutations of genes encoding members of the WNT signaling pathway — i.e., CTNNB1 (β catenin), AXIN1, and APC — have previously been reported as genetic hallmarks of ATCs (12, 13), with mutation frequencies >60% for CTNNB1 (reviewed in Supplemental Table 8). The deep-sequencing results in our cohort do not replicate these findings. We found a single tumor carrying a missense mutation (p.L347P) in CTNNB1, which is distinct from the CTNNB1–exon 3 gain-of-function hotspot previously reported in this disease. Variants in AXIN1 were found in 2/117 tumors, both without paired normals available and both of which may have represented low-frequency germline polymorphisms. A single truncating mutation in APC (p.Q1529X) was found in an ATC with an unusually high mutation burden. In addition, our results do not support a relevant role for mutations in genes in the apoptosis, Hedgehog, homologous recombination, immune response, insulin-like, JAK-STAT, tricarboxylic acid, nucleotide excision repair, polycomb, ubiquitination, or TGFβ pathways in PDTC and ATC tumorigenesis. With respect to mutations reported in single cases in a recent exome sequencing study of ATCs (21), we replicated these findings in the cyclin-dependent kinase inhibitors CDKN1B and CDKN2C (and found an additional truncating mutation in CDKN2A), as well as in other genes such as ERBB2, PTCH1, and DAXX (Supplemental Figure 1), but not in TRAF7 and NCOR1.

Gene fusions

Chromosomal rearrangements involving genes known to be translocated in thyroid tumors were frequent events in PDTCs (14%) but were absent in ATCs and did not overlap with BRAF, RAS, TSHR, or STK11 mutations (Figure 1D). RET/PTC rearrangements were detected in 5 PDTCs and involved the most common RET partners CCDC6 and NCOA4. Translocations leading to PAX8-PPARG fusions were observed in 3 PDTCs, whereas fusions involving ALK gene were detected in another 3 tumors. The kinase domain of ALK was recombined with three different upstream partners, including the known STRN and EML4, as well as CCDC149, a novel ALK fusion partner, which is a coiled-coil family gene located on chromosome 4 (fusion included CCDC149 exons 1–10 and ALK exons 20–29). The 11 PDTCs harboring gene fusions occurred in younger patients (49 vs. 58 years, P = 0.04; Supplemental Table 5).

A single ATC without known driver mutations carried a t(15;19)(q13;p13.1) translocation involving the NUT gene (NUTM1, NUT midline carcinoma, family member 1) and BRD4 (bromodomain containing 4), resulting in an in-frame NUT-BRD4 fusion (NUT exons 1–2 and BRD4 exons 14–20). It was detected in an ATC with areas of PDTC in a 34-year-old woman who underwent total thyroidectomy and laryngopharyngectomy, plus radiotherapy, and who is alive 10 years after diagnosis. It clearly represents an outlier from the clinical behavior standpoint, which matches with her unique genetic alteration, involved in large-scale chromatin remodeling (34).

Somatic CNAs

Tumor purity of PDTCs was similar to PTCs (median tumor content 74% and 72%, respectively) whereas it was much lower in ATCs (42%) (Supplemental Figure 2). By correcting for tumor purity, we greatly enhanced our sensitivity for detecting CNAs even in most of the heavily infiltrated ATCs, as well as in some PDTCs. The ability of IMPACT to call arm-level CNAs was explored in 37 tumors that were simultaneously profiled by array-CGH. As seen in Supplemental Figure 3, CNA calls were efficiently replicated in both platforms. In addition, IMPACT identified CNAs that were not detected by array-CGH, particularly in heavily infiltrated tumors.

Whereas the genome of PTCs is largely diploid, CNAs in advanced thyroid tumors were common and widespread (Supplemental Figure 4). Interestingly, CNAs were more frequent in the ATCs and PDTCs that lacked driver mutations (Figure 5A). A similar subset of PTCs from the TCGA analysis also possessed a high prevalence of CNAs in the absence of driver mutations (1). Overall, 8 arm-level recurrent CNAs (losses of 1p, 8p, 13q, 15q, 17p, 22q, and gains of 1q and 20q; Figure 5A) were present at higher prevalence, an effect that persisted at a more conservative copy number threshold (see Methods).

Figure 5

Recurrent MSK-IMPACT–derived CNAs found in 84 PDTC and 33 ATC. Representation of arm-level regions recurrently gained or lost in PDTCs and/or ATCs. CNAs were corrected for tumor purity in each sample with known driver mutations (see Methods and Supplemental Figure 2). (A) IGV representation of the altered chromosomal regions, with approximate locations shown on the top panel (genome build hg19), expressed as red (gain) or blue (loss), with shading intensity proportional to the log-ratio (lr) values. Samples are grouped by tumor type and sorted by genetic driver alteration: BRAF, RAS, fusions (RET/PTC, PAX8-PPARG, and ALK), or none/unknown. Color key and annotations are shown on the left. (B) Frequencies of the indicated CNAs in PDTCs and ATCs. Copy number gains (red) or losses (blue) were defined using two lr thresholds: ±0.1 (lighter shading) and ±0.4 (darker shading). Asterisks denote significant differences expressed as Fisher’s exact test P values for ±0.4 threshold: PDTC, 0.06 for 1p loss; ATC, < 2 × 10^–4 for 8p loss, 17p loss, and 20q gain. (C) Kaplan-Meier survival curves for PDTCs harboring chromosome 1q gain (left, log-rank P values for ±0.1 and ±0.4 thresholds are 0.03 and 0.06, respectively) and for ATCs with 13q loss (middle, P = 0.07 and 0.02) or 20q gain (right, P = 0.01 and 0.06).

Most of the CNAs occurred in a tumor type– and gene context–specific fashion (Figure 5B). Loss of 1p was marginally more frequent in PDTCs (P = 0.06), whereas 8p and 17p losses and 20q gains were far more frequent in ATC genomes (P < 2 × 10^–4 for all three). Chromosome 1p, 13q, and 15q losses were enriched in PDTCs without known driver mutations (χ² P = 0.03, 2 × 10^–3, and 3 × 10^–4, respectively), whereas loss of 22q was strongly associated with RAS-mutated PDTCs as compared with BRAF tumors (P = 1 × 10^–3) (Figure 5A).

Three out of 8 of the recurrent CNAs in the genomes of PDTCs and ATCs were associated with outcome. Patients with PDTCs with gains in chromosome 1q had worse survival (log-rank P = 0.06) (Figure 5C). In ATCs, 13q losses and 20q gains were associated with shorter survival (log-rank P = 0.02 and 0.06, respectively).

Gene expression profiling: signaling and differentiation

We used a subset of fresh-frozen specimens of 17 PDTCs and 20 ATCs to derive insights into the gene expression profiles of advanced thyroid cancers. The 37-tumor dataset was representative of the main driver genetic alterations described and included 13 BRAF, 12 RAS, 5 with alterations in other drivers (NF1, NCOA4-RET, CCDC149-ALK, NUT-BRD4, and STK11), and 7 without known driver mutations (Supplemental Table 9).

A principal component analysis efficiently separated both entities based on their global gene expression (Figure 6A). Only 2 tumors clustered out of their tumor type group: an ATC from a metastatic specimen (the single one in this subset of 20 ATCs) and a BRAF-mutant PDTC that was heavily infiltrated with macrophages (tumor purity = 46%, high M2-macrophage signature score, see also Figure 7A).

Figure 6

Principal component analysis (PCA) and BRS of 17 PDTCs and 20 ATCs. (A) Two-dimensional PCA discriminates PDTCs (squares) from ATCs (circles). Color-coding for driver alterations is shown in B. Asterisks represent ATC and PDTC outliers; see text. (B) Heatmap generated by applying the 67-gene BRS signature to advanced thyroid tumors. Expression values are displayed as Z-scores after scaling the values of each gene across the 37 samples. The 26 most informative genes are shown; the complete 67-gene signature is shown in Supplemental Figure 5. Samples are sorted by ascending BRS score: (BRAF^V600E-like on the left and RAS-like on the right) and annotated for tumor type and driver alteration. (C) Detailed comparison of driver mutation vs. BRS in BRAF- and RAS-mutant PDTCs and ATCs. Paradoxically, RAS-mutant ATCs are primarily BRAF^V600E-like (Mann-Whitney U test, P = 0.003). Box plots were generated using the Tukey method: horizontal lines within each box represent median values, box heights symbolize the IQR (IQR = Q3–Q1); Q3 and Q1 quartiles correspond to the top and bottom boundaries of the box, respectively; whiskers represent values up to 1.5 times IQR greater than Q3 (top: Q3 + 1.5 × IQR) or smaller than Q1 (bottom: Q1 – 1.5 × IQR).

Figure 7

M2 macrophage signature and TDS of 17 PDTCs and 20 ATCs. (A) Unsupervised clustering based on a 68-gene M2-macrophage signature in advanced thyroid tumors (35). Expression values are displayed as Z-scores after scaling the values of each gene across the 37 samples. The 11 most discriminatory genes (variance greater than 2) are shown; the complete 68-gene signature is shown in Supplemental Figure 6. ATCs clearly cluster apart from PDTCs consistent with their extensive macrophage infiltration. (B) Relative expression of the 16 genes of the TDS in 20 ATCs and 17 PDTCs, compared with 9 PTCs from He et al. (39) evaluated with the same mRNA array platform. ATCs have low TDS values for virtually all TDS genes, whereas PDTCs are comparable to PTCs. The 16-gene TDS signature discriminates ATCs and PDTCs (see unsupervised clustering in Supplemental Figure 7). (C) Correlation plots between TDS and BRS in PTCs from TCGA (top) and PDTCs and ATCs (bottom). Trend lines, Pearson’s correlation coefficients (r) and associated P values are shown in the graphs. TDS and BRS are positively correlated in PTCs (r = 0.74, P < 0.0001) and PDTCs (r = 0.72, P < 0.01); i.e., RAS-like tumors tend to be more differentiated than BRAF-like cancers. This relationship is lost in ATCs, which are profoundly undifferentiated (r = –0.43, P = 0.06). (D) Comparison of TDS values in BRAF- and RAS-mutated PDTCs and ATCs. Whereas ATCs are undifferentiated regardless of their driver alteration (Mann-Whitney U test, P = 0.21), BRAF-mutated PDTCs show a decrease in TDS compared with their RAS-mutant counterparts (P = 0.06). Box plots from B and D were generated using the Tukey method: horizontal lines within each box represent median values; box heights symbolize the IQR (IQR = Q3–Q1); Q3 and Q1 quartiles correspond to the top and bottom boundaries of the box, respectively; and whiskers represent values up to 1.5 times IQR greater than Q3 (top: Q3 + 1.5 × IQR) or smaller than Q1 (bottom: Q1 – 1.5 × IQR). Values outside these limits are considered outliers and are represented by dots.

MAPK signaling: the BRAF-RAS score. As a consequence of their genomic simplicity, the key oncogenic drivers of PTCs are associated with distinct biological, signaling, and gene expression properties. Pertaining to this study, the recently published TCGA analysis of 390 PTCs showed clear gene expression differences between BRAF and RAS tumors, which were used to construct a BRAF-RAS score (BRS) (1). We aimed to evaluate whether these driver-dependent gene expression characteristics persist in advanced cancers, which harbor a more complex cancer genome. Sixty-seven out of the 71 genes in the BRS were present in the mRNA array and were assessed in our PDTCs and ATCs with known driver mutation status. We found a high correlation between BRS values and BRAF/RAS mutation status (Figure 6, B and C). All 13 BRAF^V600E-mutated PDTCs and ATCs were BRAF-like. However, although RAS-mutant PDTCs were strongly RAS-like, RAS-mutant ATCs were BRAF-like (P = 3 × 10^–3), suggesting that a high MAPK transcriptional output is a characteristic property of ATCs, regardless of the driver mutation (Figure 6C and Supplemental Figure 5).

Macrophage infiltration. ATCs are known to be extensively infiltrated with macrophages (22, 23). Although it is assumed that these are M2 macrophages, which promote tumorigenesis, this has not been proven. We applied a previously characterized signature (35) of 78 genes overexpressed in M2 macrophages (68 of which were represented in our array) to the 37 tumors (Figure 7A and Supplemental Figure 6) and found that it was sufficient to discriminate ATCs from the great majority of PDTCs, which are less prone to macrophage infiltration. The only 3 PDTCs that clustered with ATCs had a lower median estimated tumor purity than the rest of the PDTCs, suggestive of significant stromal contamination (46% vs. 84%, respectively; median purity in ATCs was 36%). Consistent with this, the prevalence of likely clonal mutations in ATCs, such as BRAF and TERT, is higher in this study than that described in prior reports (refs. 10, 17–20, 36–38, and Supplemental Table 8), probably due to IMPACT’s deeper coverage.

Thyroid differentiation score. Loss of expression of thyroid differentiation markers is one of the hallmarks of advanced thyroid cancers and has profound consequences for the clinical management of these patients, who are usually refractory to radioiodine therapy due to loss of expression of NIS (sodium iodide symporter, SLC5A5) and other genes required for iodine incorporation. The TCGA analysis of PTCs used a thyroid differentiation score (TDS) consisting of 16 genes involved in iodine metabolism and thyroid specification to investigate driver-dependent effects on these parameters. We compared the TDS in PDTCs and ATCs with 8 PTCs profiled with the same platform (39). Overall, PDTCs and PTCs did not differ greatly, whereas ATCs had profoundly suppressed mRNA levels for TG, TSHR, TPO, PAX8, SLC26A4, DIO1, and DUOX2 genes (Figure 7B and Supplemental Table 9). None of the tumor types expressed THRB, DUOX1, SLC5A5, or SLC5A8. Unsupervised clustering based on the TDS discriminated ATCs from PDTCs (Supplemental Figure 7), with the exception of 3 PDTCs. These clustered with ATCs and corresponded to patients who died of the disease, two of whom were the only BRAF^V600E-mutated PDTCs in this subset.

Finally, we evaluated the relationship between TDS, BRS, and the driver mutations of these tumors. As shown in Figure 7C, TDS and BRS were correlated in PDTCs, as tumors with low TDS corresponded with BRAF-like specimens, whereas this correlation was completely lost in ATCs. Similarly, BRAF-mutated PDTCs showed a marked decrease in TDS when compared with their RAS-mutated counterparts, whereas ATCs were homogeneously undifferentiated (Figure 7D).

Most PDTCs and ATCs are thought to arise from preexisting PTCs based on their frequent co-occurrence in the same tumor specimen, where they consistently share a driver mutation (36, 40). Our results, analyzed in the context of the PTC TCGA study, provide insights into tumor microevolution and lend support to this model of tumorigenesis. Particularly compelling in this regard is the fact that mutations in the TERT promoter, which are known to activate its transcription, display a stepwise increase in frequency along the spectrum of disease progression (9% in PTCs, 40% in PDTCs, and 73% in ATCs) (1, 17–20). Interestingly, TERT mutations are subclonal in the few PTCs that harbor them, whereas they are clonal in PDTCs and ATCs, pointing to selection during tumor evolution, possibly by inducing cell immortalization. TERT promoter mutations also track with virulence of advanced disease. They are significantly associated with BRAF or RAS mutations, consistent with the proposed mechanism by which mutations in the promoter, by generating novel consensus motifs for the ETS family of transcription factors, promote TERT overexpression in cells with constitutive activation of MAPK signaling.

By using a platform with a depth of sequencing optimized to identify mutations in tumors known to be associated with abundant stromal contamination, primarily by tumor associated macrophages (TAMs) (22, 23), we identified key genetic lesions that distinguish PDTCs from ATCs (i.e., TP53, TERT, and genes encoding effectors in the PI3K pathway). This includes genetic defects that implicate functional programs not previously associated with thyroid cancer, such as the SWI/SNF complex, HMTs, and others.

With respect to the main driver alterations, BRAF mutations were less prevalent in advanced tumors compared with PTCs, whereas RAS mutations were more frequent. Rearrangements commonly seen in radiation-induced and, to a lesser extent, in sporadic PTCs (RET/PTC, PAX8-PPARG and ALK fusions) were present in a subset of PDTCs but absent in the ATCs we sampled (9, 41, 42). NF1 mutations were only found in ATCs in our series. The TCGA study of PTCs showed that BRAF- and RAS-mutant tumors exhibited profound differences in their clinical and histological characteristics and in their gene expression profile. The BRS is a 71-gene panel that distinguishes BRAF^V600E from RAS-mutant PTCs. It was highly correlated to the transcriptional output of the MAPK pathway, which was highest in BRAF-mutant cancers. This is explainable because ERK activation in RAS-mutant cells induces a negative feedback that disrupts RAF dimerization, thus attenuating pathway output. By contrast, BRAF^V600E signals as a monomer and is unresponsive to this constraint, resulting in a greater flux through the pathway (1, 43). We found that these sharp demarcations between BRAF- and RAS-mutant disease persisted in PDTC but were largely lost in ATC. PDTCs that met the standard histological definition of that entity (Turin proposal, ref. 26) were strongly associated with RAS mutations. By contrast, those defined based on the presence of high mitotic rate and necrosis irrespective of other characteristics (27) were markedly enriched for BRAF. They also had distinct clinical behaviors: BRAF-mutant PDTCs primarily developed locoregional nodal metastases, whereas RAS-mutant PDTCs presented with distant metastases. The BRS tracked with the underlying driver mutation in PDTCs but not in ATCs. This was also true for a score derived from a gene set consisting of mRNAs encoding proteins required for the differentiated function of thyrocytes (the TDS). The greater genomic complexity of ATCs may account for blurring the association of gene expression with the nature of the underlying driver mutation, in particular because of the higher frequency of mutations of genes encoding chromatin modifiers or genes that activate parallel pathways. Interestingly, even ATCs with RAS or other mutations tend to be BRAF-like, as defined by the BRS. In addition, ATCs are extensively infiltrated by TAMs. Accordingly, nonhierarchical analysis of a gene set that defines M2 macrophages clearly separated ATCs from PDTCs. It may be that the greater cellular heterogeneity of ATCs may account in part for the attenuation of the oncogenic driver effects on gene expression.

EIF1AX, which encodes for a key component of the translation preinitiation complex (PIC), is mutated in only 1% of PTCs but in approximately 10% of PDTCs and ATCs. The concordance of EIF1AX with RAS mutations is extremely strong, which is distinct from PTC, where these are largely mutually exclusive. The biological consequences of this association are currently unknown. The mutations of EIF1AX cluster in specific N- and C-terminal residues. The C-terminal p.A113 splice mutation is specific to thyroid cancer and predicts for alternative usage of a cryptic splice acceptor within exon 6, resulting in a 12–amino acid in-frame deletion. Our analysis of RNASeq data from two cases with this mutation in the PTC TCGA confirms this prediction (not shown). EIF1A plays a key role in regulating the conformation of the PIC and in scanning for the AUG initiation codon, which is disrupted in distinct ways by N-terminal and C-terminal mutations in yeast (44). Interestingly, EIF1AX mutations are mutually exclusive with alterations in any of the PI3K/AKT/mTOR pathway members, suggesting that they may confer overlapping functional gains. In addition, EIF1AX mutations are predictive of worse survival in PDTCs, providing a potentially useful marker for risk stratification in a heterogeneous disease in need of better prognostic indicators (45, 46).

Both entities differed in overall mutation burden, which was significantly higher in ATCs. Within PDTCs, a higher number of mutations was associated with larger tumors, presence of metastasis, and shorter survival. TP53 mutations, in particular, distinguished ATCs from PDTCs (73% vs. 8%, respectively). CNAs, some of which had been previously reported at very low frequencies (47–51), proved to be distinctive of each tumor type. PDTCs have a greater frequency of 1p losses, whereas 8p and 17p losses, as well as 20q gains, were more common in ATCs. Interestingly, 22q losses were strongly associated with RAS-mutant PDTCs. Loss of the 22q tumor suppressor gene NF2, which encodes for merlin, has been recently implicated in this association. Consistent with this, the combined activation of oncogenic Ras with Nf2 loss leads to development of PDTC in mice. This is because inactivation of the Hippo pathway through merlin loss leads to a YAP-TEAD–dependent transcriptional activation of oncogenic and WT RAS, thus enhancing MAPK transcriptional output and promoting transformation (1, 52).

We compared the results of our targeted cancer gene NGS (IMPACT) approach with a recently reported whole exome sequencing (WES) study of 22 ATCs (21). The greater depth of sequencing achieved by IMPACT (739× vs. 264× in the WES study) may explain the differences observed between both platforms in the frequency of TP53 (73% by IMPACT vs. 27% by WES), BRAF (45% vs. 27%), PIK3CA (18% vs. 9%), and PTEN (15% vs. 0%) and may improve detection of subclonal events. This is particularly relevant in ATCs because of their low tumor purity, which calls into question the suitability of WES as the platform of choice. Moreover, the WES approach failed to detect mutations in members of the SWI/SNF and HMT functional groups that we report here and that likely play a fundamental role in the biology of these tumors. We acknowledge, however, the drawbacks of sequencing a limited set of cancer genes. For instance, RASAL1, USH2A, HECTD1, MLH3, and MSH5, which were rarely mutated in ATCs by WES, were not included in IMPACT. Two of these, MLH3 and MSH5, are MMR genes, a functional group that we find to be disrupted in at least 12% of ATCs.

The biological consequences of the novel drivers (e.g., EIF1AX and RAS) and functional groups (SWI/SNF) in the context of thyroid tumorigenesis remain to be explored. Others recapitulate phenotypes observed in genetically engineered mouse models of advanced thyroid cancers, such as PTEN and TP53 (53), BRAF and TP53 (54), RAS and TP53 (55), RAS and NF2 (52), and BRAF and PIK3CA (56).

The findings reported here provide tools that can be leveraged to improve the molecular diagnosis of these clinical entities, many of which likely have prognostic implications. Particularly relevant is the strong association of PDTC-Turin tumors with RAS mutations, whereas PDTC-MSK tumors were strongly associated with BRAF. In addition, RAS- and BRAF-mutant PDTCs have distinct tropism for metastases, with the former tending to home at a distance, whereas the latter metastasize to locoregional lymph nodes. Moreover, the discovery in well-differentiated tumors of subclonal mutations of genes that we show to be enriched in advanced disease should raise particular concerns. This also opens a path to explore the biology of novel genetic associations that may point to tumor dependencies that can be exploited therapeutically.

Patient tissue samples. PDTC and ATC samples were randomly selected from the pathology department files of the institution from 1986–2015. ATCs were classified according to the last WHO classification of endocrine tumors, whereas PDTCs were defined as follows: (i) according to the Turin proposal, by architectural and high-grade features (mitosis and necrosis), the presence of a solid/nested/insular growth, the absence of nuclear features of PTC, and either convoluted nuclei, mitotic activity ≥ 3 × 10 high power fields (HPF), or tumor necrosis (26); and (ii) as a carcinoma displaying high mitotic activity (≥ 5 mitosis/10 HPF, ×400) and/or tumor necrosis, and showing follicular cell differentiation at the morphological or immunohistochemical level (27). MSK-IMPACT targeted sequencing was performed in all 117 tumors (84 PDTCs and 33 ATCs): 80 from formalin-fixed paraffin-embedded tissues (FFPE) and 37 from fresh-frozen material. The 37 frozen tumors, 17 PDTCs and 20 ATCs, were also expression profiled with Affymetrix U133 plus 2.0 array and with the Agilent SurePrint G3 CGH 1x1M array-CGH platform to validate copy number calls.

Single nucleotide variant and indel calling and filtering. Single nucleotide variants (SNVs) and short indels (<30 bp in length) were automatically annotated by the MSK-IMPACT pipeline, as previously described (24). Tumor samples without paired normals (11/117; 5 ATCs and 6 PDTCs) were compared against pooled normals. All variants were annotated based on the information available in catalog of somatic mutations in cancer (COSMIC; http://cancer.sanger.ac.uk/cosmic), NCBI-dbSNP (http://www.ncbi.nlm.nih.gov/snp), and the 1,000 Genomes Project (http://www.1000genomes.org/). Variants highlighted in this study were subsequently manually reviewed. For the 106 tumors with paired normals, all variants confirmed as somatic were reported, regardless of location and clonality. For tumors compared against pooled normals, MSK-IMPACT automatically called SNVs with reported frequencies <1%. We manually reviewed the 11 tumors fulfilling these criteria as follows: (i) keeping variants reported in COSMIC; (ii) removing variants reported as polymorphic (with an reference sequence code in dbSNP); and (iii) removing variants with allele frequencies that were >10% of the allelic fraction of the driver mutation in the same tumor.

MSK-IMPACT sequencing data is publicly available at the cBioPortal for Cancer Genomics (http://www.cbioportal.org/). Mutation plots were generated using the OncoPrinter (v1.0.1) and MutationMapper (v1.0.1) tools, which are available at the cBioPortal (57, 58).

Chromosomal rearrangements were called for genes whose introns were covered by MSK-IMPACT, which included most of the previously reported fusions in thyroid tumors (with the notable exceptions of NTRK1 and NTRK3).

CNAs and estimation of tumor purity. DNA CNAs were primarily called from IMPACT by comparing sequence reads of targeted regions in tumors relative to a standard diploid normal sample, as described (24). Although IMPACT targets a discrete number of exons in each chromosome arm, it efficiently identified arm-level chromosomal genetic gains and losses, as confirmed by a genome-wide Agilent SurePrint G3 CGH 1x1M array-CGH platform in a subset of 37 advanced thyroid tumors (Supplemental Figure 3), which showed excellent agreement in the copy number calls between the two methodologies.

As the macrophage infiltration of advanced thyroid tumors (particularly ATCs) can impact the sensitivity of CNA detection, we corrected CNA values for each tumor based on tumor purity. Tumor purity was calculated based on the mutant allele frequencies of clonal heterozygous somatic mutations in regions lacking overt CNAs (see also Supplemental Figure 2). For tumors with a BRAF or RAS mutation (heterozygous mutations considered clonal events) (59), purity was calculated by doubling the frequency of the alternate allele (e.g., in a tumor with a BRAF^V600E mutation frequency of 0.45, tumor purity was 90%). For tumors without BRAF or RAS mutations, we used other gene mutation frequencies as purity estimators. In tumors lacking canonical drivers, we used an average of all likely clonal events. We discarded likely subclonal events (e.g., mutations with disproportionately low allelic frequencies compared with other somatic variants). We discarded tumor suppressors (i.e., TP53, NF1) as purity estimators due to possible coexisting loss-of-heterozygosity events. Chromosomal regions harboring the drivers (typically BRAF or RAS) were confirmed to be diploid before using those genes as purity estimators. Altogether, we derived tumor purity data from 75 BRAF/RAS-mutated tumors and used alternative genes for 19 others. Purity could not be estimated with confidence for 23 tumors (22 PDTCs and 1 ATC).

Segmented genome-wide copy number was corrected for tumor purity in the following manner. For all segments where r is the log₂ ratio of segmented copy number, then R = 2^r and α is the estimated tumor purity, the copy ratio in cancer cells of a cellularly heterogeneous sample is therefore: R′ = (R/a) − ([1 – a]/a) for all r > log₂(1 − a), otherwise R′ = (r/min[r]) min (R′).

The final purity-corrected copy number is log₂(R’) for all segments in the specimen. CNAs were determined using the purity-corrected segmented data utilizing one of two thresholds (±0.1 and the more conservative ±0.4 threshold). All analyses were expressed at ±0.4, with the exception of survival analyses. As focal events were few, the analyses described here focused on larger and overt arm-level CNAs. Given the challenge of definitively establishing thresholds of heterozygous loss and homozygous deletions in impure tumors (signal compression), our analysis excluded homozygous CNAs. CNAs were visualized in the Integrative Genomics Viewer (IGV), version 2.3.57 (https://www.broadinstitute.org/igv/) (60).

Gene expression. We performed mRNA expression on 37 tumors using the Affymetrix U133 plus 2.0 array. Expression normalization was performed using the gcRMA method (61). Statistical analyses were performed using the R statistical programming language (62). The lists of genes for the BRS and TDS scores were obtained from the TCGA study on PTC (1). Of the 71 BRS genes, 67 were present on our array and were used to compute the BRS score as described (1). We performed sample-wise scaling on the gene-expression profiles of the 9 PTC tumors, which were derived from a published study (39). The lists of differentially expressed genes between M1 and M2 macrophages (35) were from the MSigDB (63). The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (GEO) (64) and are accessible through GEO Series accession number GSE76039.

Clinical associations. Statistical analyses were performed using STATA/IC (version 12; StataCorp LP). Pathological characteristics were only available for resected ATCs. PDTCs had complete clinical and pathological data. Survival was recorded as of April 2015. All associations reported in PDTCs remained significant when only considering PDTCs based on the Turin criteria definition.

Statistics. Distribution of mutation frequencies in PDTCs versus ATCs was assessed by Fisher’s exact tests. Copy number distribution between multiple groups was evaluated with χ² tests. Statistical analyses and graphic representations of mutation and CNA distribution were performed on GraphPad Prism 6.02 (GraphPad Software). Demographic and clinico-pathologic characteristics were compared using Pearson χ² test for categorical variables and 2-tailed Student’s t test or Mann-Whitney U test for continuous variables. Survival analyses were conducted with Cox proportional hazards models. Proportionality assumptions were tested with Schoenfeld residuals and log-log plots. Kaplan-Meier survival curves were built, and the log-rank test was used to assess for significance of the surviving function. A P value less than 0.05 was considered significant.

Study approval. The study was approved by the Institutional Review Board of MSKCC.