Diverse Oncogenic Fusions and Distinct Gene Expression Patterns Define the Genomic Landscape of Pediatric Papillary Thyroid Carcinoma

Papillary thyroid carcinoma (PTC) is the most common malignancy in adolescents and young adults (1, 2) and incidence has been increasing across all age groups for the past three decades (3). Pediatric PTC (PPTC), diagnosed among individuals 18 years and under, constitutes less than 2% of all PTC diagnoses; thus, the biology and molecular genetics of pediatric disease has been sparsely characterized relative to adults.

Clinical features of PTC in children differ from those in adults, recognition of which led to the development of pediatric-specific guidelines, first published in 2015 (4). Children with PTC present with more advanced disease, larger primary nodules, and more aggressive features than adults (5). Lymph node involvement and distant metastases are more common in the pediatric population (6, 7). Recurrence rates are also higher, reaching up to 40% in children versus 20% in adults (8, 9). These differences may reflect the host environment, but may also reflect distinct molecular mechanisms underlying oncogenesis and invasiveness in children (10).

In 2014, The Cancer Genome Atlas (TCGA) adopted an integrated multiplatform analysis to define the genomic landscape of 496 PTCs (11). This study demonstrated that PTC is primarily driven by mutually exclusive substitutions in BRAF or RAS (which accounted for 59.7% and 12.9%, of tumors respectively) or by oncogenic fusions, most commonly between RET and either CCDC6 (RET-PTC1) or NCOA4 (RET-PTC3), which comprised an additional 5.3%. Rarer substitutions and fusions were identified in a further 12% of PTCs. These findings affirmed the importance of aberrant signaling through the MAPK and PI3K pathways in thyroid oncogenesis. TCGA further described gene expression patterns according to (i) impact on markers of thyroid differentiation; (ii) activation of the MAPK-ERK signaling pathway; and (iii) a novel BRAF-RAS score (BRS) that designated tumors as BRAF-like or RAS-like. A fundamental shortcoming of TCGA, however, with respect to PPTC, was the inclusion of only 3 tumors from individuals less than 18 years and no individuals less than 12 years at diagnosis.

Molecular characterization of pediatric disease has been slower in arriving. Several studies have adopted PCR-based genotyping or targeted next-generation sequencing (NGS) panels to describe PPTC (12–18), summarized nicely by Bauer (19) and Paulson (20), although such analyses are limited to relatively small cohorts and finite numbers of candidate genetic variants. Nevertheless, they identified that the distribution of oncogenic variants in children differs from those in adults, with substantially lower substitution rates (specifically BRAF^V600E and RAS substitutions). More recently two studies have employed targeted next-generation gene panels, identifying oncogenic drivers 69% to 77% of tumors (18, 21). Oncogenic fusions appear to exist in excess, relative to substitutions (13, 15, 21). Finally, there is a significantly greater proportion of PPTC without an identified genomic driver. We recently summarized these studies and found that of 338 tumors analyzed, only 164 (49%) were associated with a genetic alteration (22) while the remainder constituted “dark matter” tumors. This is a meaningful deviation from TCGA, wherein 89.9% of tumors was associated with a solitary oncogenic variant and only 10.1% were dark matter.

The current study sought to identify the oncogenic drivers of pediatric dark matter PTCs using whole-exome sequencing (WES) and transcriptome-wide NGS [RNA sequencing (RNA-seq)], to characterize gene expression patterns in pediatric PTC and to examine these in contrast to adult PTC.

This study was approved by the Institutional Research Ethics Board and all patients provided informed written consent for participation. Thyroid nodules were snap-frozen from 45 consecutive pediatric patients at our institution. Peripheral blood lymphocytes were collected for germline DNA extraction. An additional 12 archival tumors were obtained from the Cooperative Human Tissue Network, although sample quantity or quality limited analysis to genotyping alone for all but 4 of these [3 PTC and 1 follicular thyroid carcinoma (FTC)]. Five histologically benign lesions (1 follicular adenoma, 4 follicular nodular disease) were included for gene expression analysis. DNA and RNA were extracted using the AllPrep Universal Kit (QIAGEN).

All surgical specimens were reported utilizing a standardized synoptic report. To ensure accurate tumor subclassification using the 2017 WHO classification of thyroid neoplasms (23), with the exception of 4 cases, all specimens were rereviewed by 2 pathologists with endocrine pathology expertise (O. Mete and R. Chami). Four samples lacked diagnostic-quality archival material or images, and thus we relied on the original pathologist's report for diagnosis. Histologic variant could not be determined for 3 PPTCs based on available images and original pathology report.

Allele-specific quantitative PCR or quantitative RT-PCR, using Thyroid Cancer Mutation Detection Kit and the Thyroid Cancer Fusion Detection Kit (Entrogen Inc), was performed to genotype all 52 carcinoma samples for 17 of the most common genetic alterations found in adult thyroid carcinoma (alterations in BRAF, HRAS, KRAS, NRAS, RET, and PAX8 genes). DICER1 analysis was previously reported (22). The human telomerase reverse transcriptase (hTERT) promoter was generated as a single 220 bp amplicon using primers FWD-AGCGCTGCCTGAAACTCG and hTERT-REV-CACAGACGCCCAGGACCG and custom TaqMan probes for the C228T and C250T variants and 7deaza-dGTP modification as described (24). Digital PCR was performed on a QX200 droplet digital PCR System (BioRad Inc) according to manufacturer's directions. Samples were prepared with 2x ddPCR Supermix, forward and reverse primers, and partitioned into a median of approximately 14,000 droplets per well in the QX200 droplet generator. After PCR amplification, samples were read on the QX200 droplet reader using BioRad QuantaSoft software.

Exome and transcriptome sequencing were performed by The Centre for Applied Genomics (TCAG) using established protocols.

Sequencing information, raw FASTQ files and clinical characteristics for the 390 TCGA PTC samples were obtained from the TCGA data portal (https://portal.gdc.cancer.gov/). The files were downloaded and processed using the FusionValidator pipeline (described below).

Exome libraries were prepared using Agilent SureSelect Human Exome Library Preparation V5 kit and the Agilent Bravo Automation System followed by paired-end sequencing on an Illumina HiSeq 2500 platform. Genomic DNA (750 ng) was fragmented to 200-bp on average using a Covaris LE220 instrument. Sheared DNA was end-repaired and the 3′ ends adenylated prior to ligation of adapters with overhang-T. Genomic library was amplified by PCR using 10 cycles and hybridized with biotinylated probes that target exonic regions; the enriched exome libraries were amplified by an additional 8 cycles of PCR. Exomic libraries were validated on a Bioanalyzer 2100 DNA High Sensitivity chip (Agilent Technologies) for size and by qPCR using the Kapa Library Quantification Illumina/ABI Prism Kit protocol (KAPA Biosystems) for quantities. Exome libraries were pooled and sequenced with the TruSeq SBS sequencing chemistry using a V4 high throughput flowcell on a HiSeq 2500 platform following Illumina's recommended protocol. Approximately 6 to 8 gigabases of raw paired end data of 126-bases were generated per exome library.

Quality of total RNA samples was determined using the Agilent Bioanalyzer 2100 RNA Nano chip. Concentration was measured by Qubit RNA HS Assay on a Qubit fluorometer (Thermo Fisher). RNA library preparation was performed following the NEB NEBNext Ultra II Directional RNA Library Preparation protocol. Briefly, 400 ng of total RNA were used as the input material and enriched for poly-A mRNA, fragmented into the 200- to 300-base range for 4 minutes at 94°C and converted to double stranded cDNA, end-repaired, and adenylated at the 3′ to create an overhang A to allow for ligation of Illumina adapters with an overhang T; library fragments were amplified under the following conditions: initial denaturation at 98°C for 10 seconds, followed by 15 cycles of 98°C for 10 seconds, 60°C for 30 seconds, and 72°C for 30 seconds, and finally an extension step for 5 minutes at 72°C; at the amplification step, each sample was amplified with a different barcoded adapters to allow for multiplex sequencing. One microliter of the final RNA library was loaded on a Bioanalyzer 2100 DNA High Sensitivity chip (Agilent Technologies) to check for size; RNA libraries were quantified by qPCR using the Kapa Library Quantification Illumina/ABI Prism Kit protocol (KAPA Biosystems). Libraries were pooled in equimolar quantities and paired-end sequenced on 2 lanes of a High Throughput Run Mode flowcell with the V4 sequencing chemistry on an Illumina HiSeq 2500 platform following Illumina's recommended protocol to generate paired-end reads of 126-bases in length.

WES paired reads were aligned to the human genome (hg19/GRCg37) using BWA-MEM (v.07.8). PCR duplicates were identified with Picard MarkDuplicates (v.1.108). Somatic mutations were detected using MuTect2 as part of GATK v.3.5 and Delly v.0.7.1 and annotated with ANNOVAR. Custom filters were implemented to increase specificity of detection as described (25). The mutational burden was calculated by dividing the total number of exonic variants identified by the SSM pipeline by the total capture region of an exome, which is on average 30Mb.

Gene fusions were detected using an in-house pipeline, FusionValidator (25). In brief, transcripts were identified using an approach that integrates five independent fusion algorithm detection tools: Defuse 0.6.2, ChimeraScan 0.4.5, STAR Fusion 0.7.0, MapSplice 2.1.9, and FusionCatcher 0/99.4d_beta (26–30), and then performs a dynamic realignment against a large data set of approximately 1,300 nonneoplastic samples from 43 different tissues (31) to remove recurrent artifacts and read-throughs. Putative fusions were then validated in silico by transcript assembly. Finally, a machine-learning predictive model was applied to all fusions to generate a final predictive score (the result of a function based on a linear combination of predictor variables). Gene fusions were annotated and Circos plots generated using FusionHub (32). Breakpoints and validation scores are summarized in Supplementary Table S1.

RNA-seq reads were aligned on human genome using STAR 2.4.2 and HT-Seq was used to count reads aligning on every single gene, after removing PCR duplicates and reads mapping on ribosomal RNA, miRNA, and small nucleolar RNA. Raw gene counts were then normalized using the trimmed mean of M-value (TMM) method in the EdgeR package, on genes with at least one read per million bases in at least 3 samples.

K-means clustering was performed on normalized gene expression values and two well established methods were used to identify the number of distinct clusters present across the cohort and determine the optimal number of clusters: the Elbow Method (33) and the Gap Statistic (34). For the elbow method the within-cluster sum of squares for k = 1 to k = 15 was extracted and the wss values per cluster were plotted. For the Gap statistic a max k of 10 with 50 bootstraps was used.

Differentially expressed genes in each cluster were identified using the Limma package in R³⁸. Unsupervised hierarchical clustering with principal component analysis was used to exclude batch effects. P values were adjusted for multiple testing using the Bonferroni correction for FDR. Differentially expressed genes were identified between each cluster and the remainder of the cohort and clusters were considered statistically significant if FDR ⇐ 0.05.

To ascertain whether gene expression differed between pediatric and adult PTC, we first performed unsupervised clustering including all 37 tumors from the present cohort as well as 389 tumors from TCGA. We further performed a generalized linear model likelihood ratio test to compare within-variant expression for the most common oncogenic variants, BRAF^V600E (6 pediatric vs. 233 adult tumors), CCDC6/RET (3 pediatric vs. 13 adult), NCOA4/RET (4 pediatric vs. 4 adult), and ETV6-NTRK3 (3 pediatric vs. 5 adult). The latter two analyses lacked statistical power to discriminate gene expression (with only 0 and 4 differentially expressed genes, respectively). P values for the generalized linear model test were adjusted for multiple testing using the Benjamini–Hochberg method for controlling the FDR.

Raw exome data files are available at the European Genome-Phenome Archive, accession EGAD00001007502 and RNA-seq data EGAD00001007499. Source data for gene expression are accessible at https://doi.org/10.5683/SP2/IANCYA.

Thyroid differentiation score (TDS) and ERK score were generated as previously described using Ht-Seq–derived expression scoring (35). The TDS is defined by 16 thyroid-specific genes, associated with thyroid metabolism and function and was determined as the mean of the median-centred rlog value for the 16 genes.

The ERK score was defined by a 52 gene signature derived to reflect MAPK activation levels. Fifty of the 52 ERK genes were identified in our analysis and used to determine ERK score. KIR3DL2 and GTF2A1 L were not identified and were not used in ERK score derivation in the current analysis. Both are associated with ERK upregulation. ERK score was defined as the mean of the median-centered rlog values for the 50 genes. Independently, we derived a MAPK score using a nonparametric, unsupervised method.

Analyses were performed using R software 3.5.3, RStudio 1.1.463, GraphPad Prism 7.0 or SPSS v.26. Comparisons of categorical variables were performed via χ² or Fisher exact test. Comparisons of continuous variables were performed via independent t test. ANOVA tests were conducted when continuous variables between more than 2 groups was required, with Tukey correction used for multiple comparisons. The limit of significance for all analyses was defined as P-value < 0.05 using Bonferroni correction where applicable. Kruskal–Wallis and Wilcoxon rank–sum tests were used to assess significant difference of TDS and ERK scores where indicated.

Fifty-one patients with PPTC, 1 with FTC, and 5 with benign nodular disease (BND) were enrolled over an 8-year period. Baseline characteristics are summarized in Table 1, Supplementary Table S2 and Fig. S1. The median age at diagnosis was 13 years (range 5–17). The female:male ratio was 1.04:1. One individual had been exposed to radiation for myeloablation prior to stem-cell transplant. The classical variant of PTC was most common, comprising 33 of 51 PTCs (65%).

We segregated the cohort into 3 age categories: predominantly prepubertal (0–10 years), peri-pubertal (10.1–14 years), and predominantly pubertal (14.1–18 years). Younger patients were more likely to have locoregional lymph node metastases (P = 0.012), distally metastatic disease (P = 0.001), and angioinvasion (P = 0.0016). There were no differences in tumor size (P = 0.126) among the three age categories (Fig. 1, Table 1).

Targeted genotyping of the 17 most common thyroid cancer–associated genetic alterations identified variants in 31 of 51 (61%) of malignancies. Substitutions in BRAF, NRAS, and gene fusions (RET and PAX8) were identified in 11 (22%), 1 (2%), 18 (35%), and 1 (2%) tumors respectively. The remaining 20 tumors (39%) comprised the so called “dark matter” samples and were evaluated using WES and RNA-seq.

In addition to the 20 dark-matter samples, an additional 5 BND and 12 tumors with known variants (6 BRAF^V600E, 1 NRAS^Q61R, and 5 RET fusion samples) were included for the purposes of gene expression analysis and to serve as positive controls. In total, a cohort of 37 tumors, enriched for tumors without common genetic variants, underwent exome and transcriptome sequencing.

The mutational landscape of these 37 tumors is summarized in Figs. 2A–E. Substitutions accounted for 3 (15%) driver-unknown tumors (Fig. 2B) and gene fusions were identified in 15 (75%; Fig. 2C), and are further described below. There were no pathogenic variants identified in 5 BND, 1 PTC, and 1 FTC (Figs. 2D and E). Altogether, pathogenic substitutions and gene fusions were identified in 30 of 31 (97%) of these next-generation sequenced PPTCs and 49 of 50 (98%) of the overall cohort of PPTCs (excluding the solitary FTC).

WES identified somatic variants in DICER1 in two samples, as previously described (22). Additionally, 1 patient carried a GNAS (c.636G>T:p.Q212H) substitution, which has been previously reported (36), albeit not in PTC. All six of the control BRAF^V600E and one NRAS^Q61R substitution were confirmed by WES. No other candidate driver substitutions or insertion–deletions (indels) were identified in the remaining samples. Thus, WES identified 3 tumors with substitutions that were not identified by the initial genotyping. The mean mutational burden was estimated at 0.5 mutation/Mb.

With one exception, substitutions were mutually exclusive, both with each other and with oncogenic fusions. The aforementioned GNAS^Q212H substitution was identified in a tumor derived from a child with a pathogenic germline PTEN substitution (c.710delA).

Although canonical substitutions in the hTERT promoter (C228T and C250T) have been associated with poor prognosis in adults, particularly in association with BRAF^V600E (37), these hTERT promoter substitutions appear to be rare in children (38, 39). Consistent with prior studies, none of the tumors in this cohort bore either of the canonical hTERT promoter substitutions.

Oncogenic fusion transcripts were identified in 20 of 31 (65%) PPTCs evaluated by RNA-seq and no benign lesions (Fig. 2C, F, and G and Table 2). In addition to the known RET-CCDC6 and RET-NCOA4 transcripts, RNA-seq identified several pathogenic oncogenic fusion transcripts infrequently or not previously reported in adult PTC (AKAP13-RET, TRIM27-RET, TNIP1-RET, ETV6-NTRK3, SQSTM1-NTRK3, TFG-MET, TRIM24-MET, EML4-ALK, and TG-THADA) and two novel oncogenic fusions (STRN-RET and TG-PBF). There were no BRAF fusions. Two of the dark-matter samples expressed RET-CCDC6 and RET-NCOA4 and reflect ‘false negatives’ from the initial qPCR genotyping.

Primary data from the 2014 TCGA project were reanalyzed using an in-house bioinformatic pipeline (25). This pipeline accurately identified all oncogenic fusions reported in the initial TCGA analysis, as well as two additional presumed oncogenic fusion transcripts (TPM3-NTRK1 and RAF1-AGGF1) among 47 tumors that had previously been categorized by TCGA as dark matter.

We evaluated the relationship between oncogenic driver and patient/tumor characteristics. Clinical features were first analyzed according to variant class (fusion vs. substitution). The solitary PPTC without oncogenic variant and 1 FTC were excluded from this analysis. Data are summarized in Table 2. The mean age was 14.0 ± 2.5, 12.0 ± 2.9, and 14.6 ± 2.0 years for patients with tumors driven by fusions, substitutions and for benign lesions, respectively (P = 0.033).

Fusion-driven tumors were more likely to be associated with locoregional lymph node metastases (P = 0.001), while variant class was not significantly associated with differences in the rates of distant metastasis (P = 0.338) or angioinvasion (P = 0.116).

Unsupervised hierarchical clustering analysis of the 1,000 most variably expressed genes identified 3 distinct gene expression clusters, shown in Fig. 3A and B. The presence of 3 unique clusters was confirmed independently by k-means clustering analysis using the Elbow and Gap Statistic methods (Supplementary Fig. S1A and S1B, respectively). These 3 clusters correspond broadly to (i) non-BRAF^V600E substitution-driven tumors; (ii) BRAF^V600E and RET-NCOA4 driven tumors; and (iii) RET-CCDC6, and other fusion-driven tumors (Table 3 and Supplementary Table S3). In contrast to variant class, gene expression cluster was strongly associated with histologic variant (P < 0.001). There were no between-cluster differences in size (P = 0.973), distal metastases (P = 0.333), or angioinvasion (P = 0.834).

To ascertain whether gene expression in pediatric tumors differed markedly from adult PTC, we applied the same clustering analysis to the entire set of pediatric tumors in addition to 389 PTCs from TCGA (Fig. 4A). Globally, pediatric tumors did not segregate from the adult tumors based on gene expression. When we compared adult and pediatric tumors driven by the same variant, however, CCDC6-RET–driven pediatric tumors showed separable gene expression clustering (Fig. 4B), while BRAF^V600E-driven tumors did not (Fig. 4C).

We then assessed the impact of genomic variants on the expression of genes related to thyroid differentiation and function. Utilizing a 16-gene TDS (11), tumors expressing the BRAF^V600E substitution demonstrated the lowest expression of genes associated with the differentiated state, while BND, FTC, and those expressing TG-PBF and PTEN substitution showed highest expression of these genes (Fig. 5A and B). Tumors in Cluster 2 differed significantly from those in Cluster 3 with respect to expression of DUOX2, PAX8, SLC5A5 (sodium-iodide cotransporter), and thyroid hormone receptor alpha (THRA). Pairwise comparisons of gene expression amongst all genes comprising the TDS is shown in Supplementary Fig. S2.

Signaling through the BRAF–RAS–MAPK–ERK pathway was also assessed, as previously described (11), and demonstrated a reciprocal pattern with respect to the TDS, inasmuch as less differentiated tumors demonstrated higher RAS–MAPK–ERK signaling (Fig. 5C). Overall, Cluster 2 and 3 did not differ significantly here, although both were significantly different when compared with Cluster 1 (Fig. 5D).

As there was only a single RAS-associated tumor in this cohort, we were unable to independently derive a BRS. Nonetheless, gauging expression of 67 of 71 genes from the TCGA BRS, allowed segregation into BRAF-like and RAS-like tumors. Tumors associated with Cluster 1 (RAS, PTEN/GNAS, DICER1) demonstrated RAS-like profiles, while (with one exception), Cluster 2 and 3 tumors are enriched in tumors with a BRAF-like profile (Supplementary Table S4).

Finally, we applied Gene Ontology, Gene Set Enrichment Analysis (GSEA), and Gene Set Variation Analysis (Supplementary Fig. S3) using the Kyoto Encyclopedia of Genes and Genomes pathway database (40). The 20 most significantly enriched pathways in each cluster (FDR q-value < 0.01, P-value < 0.0001) were identified by GSEA and are summarized in Supplementary Table S3. Pathways enriched in Clusters 2 and 3 showed significant overlap, including upregulation of pathways related to cell adhesion, extracellular matrix, cytokine–cytokine receptor interaction, MAPK signaling, and p53 signaling, while these were largely downregulated in Cluster 1.

This is the first study to apply unbiased whole-exome and whole-transcriptome analysis to an exclusively pediatric cohort of thyroid tumors. In so doing, it defines a genomic landscape distinct from that in adults, dominated by diverse oncogenic fusions with a substantially smaller proportion of BRAF^V600E-driven tumors.

Oncogenic drivers were identified in 98% of tumors, with the majority of the previous unknowns driven by rare oncogenic-fusion transcripts. These results complement the landmark TCGA project wherein only 3 tumors were derived from individuals under 18 years of age; and to other pediatric series that adopted targeted genotyping approaches for the most common variants in adult PTC and demonstrated that approximately 50% of pediatric tumors (41, 42) lacked identified oncogenic variants.

We found that 12 of 32 (38%) PPTCs were driven by fusions that are not among the common adult PTC drivers. Nevertheless, 10 of 12 fusions incorporated an RTK fused to the 3′end of the transcript, a finding with meaningful therapeutic implications, as discussed below.

Given the diversity of fusion transcripts in PPTCs, we propose that the unbiased nature of full transcriptome RNA-seq establishes it as the de facto gold standard for genomic analysis of pediatric thyroid tumors. With decreasing costs and widespread availability, we propose this as an appropriate ‘universal’ test. Moreover, RNA-seq analysis affords the opportunity to evaluate variations in gene expression associated with oncogenic variants, lending biological context to genomic alterations, while targeted RNA analysis methodologies (whether qPCR or NGS-based) are primarily geared towards detection of known fusions, rather than novel transcripts and expression analysis.

As with prior targeted analyses, we identified a predominance of oncogenic fusion transcripts underlying PPTC, suggesting that pediatric tumors arise primarily from double-stranded DNA breaks, rather than errors in replication, lending insight into the primary oncogenic events and/or failed repair mechanisms that drive these tumors. In addition to the canonical RET-CCDC6, RET-NCOA4, and ETV6-NTRK3, several rare fusion transcripts, uncommon to PTC (although identified in other malignancies) were identified. These include AKAP13-RET (11), TRIM27-RET (43), TNIP1-RET (44), SQSTM1-NTRK3 (45), TFG-MET (46), TRIM24-MET (47), EML4-ALK (48), and TG-THADA (11). In addition, the current study identified 2 novel oncogenic fusion transcripts, STRN-RET and TG-PBF (PTTG1IP).

STRN is located on 2p22.2 and has been reported in the context of PTC when fused to ALK, although it has not been previously reported in association with RET. STRN-ALK can transform cells in vitro and results in tumor formation in mice (49, 50). The STRN-RET fusion retains both the coiled-coil domain of STRN and the tyrosine-kinase domain of RET. As such, we suggest that this novel fusion likely retains the potential to drive tumorigenesis.

TG-PBF, has not been described in human tumors, although overexpression of wild-type PBF, has been associated with tumor induction in xenograft models (51) and a transgenic mouse expressing PBF under the influence of the thyroglobulin promoter demonstrates a pattern of follicular epithelial proliferation and radio-iodine resistance (52). PBF was previously identified as a cancer driver across 34 tumor types (53) and PBF overexpression in human thyroid carcinoma is associated with impaired DNA-damage response and poorer patient outcomes (54). In the solid-variant PPTC described herein, PBF was among the 0.5% most highly expressed transcripts. Thus, although lacking the classic features of an oncogenic (RTK-associated) fusion transcript, the TG-PBF fusion transcript, driving high expression of PBF, may indeed represent a novel oncogenic variant.

Oncogenic variants in well-differentiated thyroid carcinoma are almost invariably mutually exclusive (11, 55). We identified a solitary patient with a germline PTEN variant and a somatic GNAS variant. To our knowledge, alterations in these 2 genes have not previously been described together in PTC, however, Kirschner and colleagues previously demonstrated FTCs developing in mice doubly deficient in PTEN and PRKAR1a (which lies downstream of GNAS activation; refs. 56, 57), emphasizing the impact of convergent activation of the PTEN-mTOR and GNAS-PKA signaling pathways in follicular cell transformation.

FTC is rare in children and this cohort included only a single case. We did not identify a genomic variant in this tumor, although the BRS estimate defined this as a RAS-like tumor, and other series have described enrichment for pathogenic RAS substitutions among FTC (55).

Gene expression analysis identified 3 distinct clusters. Interestingly, carcinomas driven by DICER1, PTEN/GNAS, and NRAS demonstrated expression profiles more similar to BND than to tumors driven by BRAF^V600E or by oncogenic fusions (Fig. 3) and were enriched in immunoglobulin gene expression, as has been recently described (58). Additionally, we were intrigued to find that RET-CCDC6 and RET-NCOA4 fell into separate gene expression clusters. To our knowledge, this is the first report of major differences in gene expression between these two common RET-fusions in humans, and may explain differences in tumor behavior as has been described previously (59). Classical variant histology constituted 0% of Cluster 1 tumors, 100% of Cluster 2 tumors, and 56% of Cluster 3 tumors. Thus, among molecular tests, gene expression appears a better discriminant of lesional architecture than oncogenic variant class (Tables 2 and 3; Supplementary Fig. S4A and S4B).

Cluster 2 tumors (including those driven by BRAF^V600E) demonstrate the lowest expression of markers of differentiation [including genes necessary for iodine transport (SLC5A5) and organification (TPO, SLC26A4)] and highest levels of ERK signaling. This may have meaningful therapeutic ramifications with respect to ability of tumors (and metastases) to concentrate radioactive iodine (RAI), as well as response to RAS–RAF–MEK–ERK inhibition, although clinical validation of these associations is certainly warranted. Conversely, Cluster 1 tumors are functionally well differentiated and are rarely associated with metastatic disease.

In light of the clinical differences between pediatric and adult PTC, we anticipated that distinct gene-expression patterns would distinguish tumors between these the pediatric and adult populations. Although we did not find this to be the case overall (Fig. 4A), when focusing exclusively on tumors driven by CCDC6-RET (Fig. 4B), gene expression in children was separable from that in adults. This was not true among BRAF^V600E-driven tumors (Fig. 4C). While these studies were underpowered to explore these differences in depth, we note that children with CCDC6-RET-driven tumors were younger than those with BRAF^V600E-driven tumors (mean age 10.1 vs. 14.5 years). Clinical differences between younger children and adolescents have been described, as shown in Fig. 1 and prior studies (60, 61), and this finding would support the hypothesis that tumors arising in young (pre-/peri-pubertal) children constitute a biologically and genomically unique category, although further study is clearly necessary. These hypothesis-generating analyses were limited by small sample numbers and merit confirmation with larger sample sizes to establish generalizability.

The prognosis for children with PTC is highly favorable. Nonetheless, 15% to 25% of children present with distally metastatic disease, and of these, only 0% to 16% achieve structural or biochemical remission with current therapies (7, 62–64). In this series, the mean age of individuals with distally metastatic (pulmonary) disease was 10.4 years with 8 of 11 M₁ tumors driven by oncogenic fusions. Although only a small proportion of children develop progressive, RAI refractory and/or symptomatic disease, therapeutic implications for such patients, based on the current study are significant. In addition, 58% of children in this study (and 73% with M₁ disease) had tumors driven by RET, NTRK, ALK, or MET fusions, all of which are targetable with effective systemic therapies: selpercatinib and pralsetinib (RET), larotrectinib and entrectinib (NTRK) and crizotinib (ALK and MET), while 11 of 52 (21%) carried the BRAF^V600E substitution (dabrafenib and vemurafenib). Overall 79% of tumors (and 91% of distally metastatic tumors) were driven by variants with currently available systemic therapies. Thus, we propose that comprehensive molecular analysis be considered for all children with distally metastatic tumors that are either progressive and/or radio-iodine refractory.

There were several limitations to this study. First, given the rarity of pediatric thyroid carcinoma, sample availability was limiting, particularly among the youngest patients. Nonetheless, by prospectively ascertaining high-quality frozen samples, we were able to generate robust genomic data, whereas retrospective analysis of archival samples often limits the quality of data that can be generated using genome-wide approaches. As has been noted in other series, BRAF^V600E-driven tumors are infrequent in children and RAS-driven tumors even more so (41, 42, 65). This precluded the independent derivation of a BRAF-RAS score to characterize tumors in the current study. Many of the fusions identified herein have been previously described and associated with thyroid or other malignancies. For novel fusions, notwithstanding biologic plausibility, functional validation is necessary prior to formally ascribing oncogenic potential. Finally, the majority of patients included in this cohort underwent surgery within the past 5 years and many transitioned care upon reaching adulthood, thus we are unable to associate intermediate or long-term outcomes, including recurrence, with underlying genotype. This, however, remains an important area for future investigation, as theoretically, tumor genotype could aid in prognostication and potentially in treatment planning.

Further association of the genomic drivers of PPTC with salient clinical outcomes including radio-iodine responsiveness and achievement of a state of undetectable disease, is an important next step in this process. Such determinations, will require broader sequencing of pediatric tumors and longitudinal studies relating genotype and gene expression to these outcomes.

This study provides an exome- and transcriptome-wide analysis of genomic alterations in childhood PTC. It highlights the genomic distinctions from adult PTC and raises the possibility that these distinctions explain clinical differences that have long been reported between the two. These results provide meaningful insights into the underlying biology of disease and into treatment options for children with progressive and distally metastatic disease.

J.D. Wasserman reports personal fees from Bayer Inc. outside the submitted work. No disclosures were reported by the other authors.

A. Stosic: Formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. F. Fuligni: Software, formal analysis, supervision, writing–review and editing. N.D. Anderson: Formal analysis, investigation, methodology. S. Davidson: Formal analysis, supervision, methodology, writing–review and editing. R. de Borja: Software, formal analysis, investigation, methodology. M. Acker: Resources, project administration. V. Forte: Resources, writing–review and editing. P. Campisi: Resources, writing–review and editing. E.J. Propst: Resources, writing–review and editing. N.E. Wolter: Resources, writing–review and editing. R. Chami: Formal analysis, investigation, writing–review and editing. O. Mete: Formal analysis, investigation, writing–review and editing. D. Malkin: Conceptualization, resources, writing–review and editing. A. Shlien: Conceptualization, resources, software, formal analysis, supervision, funding acquisition, visualization, methodology, writing–review and editing. J.D. Wasserman: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–original draft, writing–review and editing.

This work was funded by the C¹⁷ Research Network (J.D. Wasserman) and SickKids Foundation (J.D. Wasserman). A. Stosic was the recipient of a Restracomp Graduate Scholarship. A. Shlien is partially supported by an Early Researcher Award from the Ontario Ministry of Research and Innovation, the Canada Research Chair in Childhood Cancer Genomics, and the Robert J. Arceci Innovation Award from the St. Baldrick's Foundation. We are grateful to Sholeh Ghayoori for assistance with the initial stages of this investigation. Anita Villani provided critical review of this manuscript, for which we are also grateful. Tissue samples were provided by the Cooperative Human Tissue Network, which is funded by the National Cancer Institute as well as the Children's Oncology Group (study ARAR15B1). Other investigators may have received specimens from the same subjects. The hTERT promoter assay was developed and validated by Tara Paton and Nunio de Rocha Oliveira Nunes.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.