A retrospective caselist of thyroid cancer patients collected at our Institute was selected based on (i) confirmed diagnosis of follicular cell-derived thyroid tumor and (ii) availability of residual archive material from both tumor and matched non-neoplastic thyroid (NT), included as patient-specific control for DNA sequencing. The obtained cohort included both various PTC histological subtypes and PDTCs. PTCs were classified according to the WHO classification of endocrine tumors (38) and PDTCs according to the Turin proposal (39). Formalin-fixed paraffin-embedded (FFPE) tissue blocks were obtained, and hematoxylin and eosin (H&E)-stained slides were reviewed by an experienced pathologist (PB); when necessary, the areas of interest were manually microdissected. Primary tumor and patient-matched NT were obtained from 57 patients and subjected to nucleic acid extraction and quality control, for a total of 126 processed samples.

The study was approved by the Independent Ethics Committee of Fondazione IRCCS Istituto Nazionale dei Tumori (protocol INT DI-20/12/13-0006020), and informed consent was obtained from patients.

Nucleic acids were extracted from unstained FFPE tissue serial sections consecutive to the pathologically revised H&E. Genomic DNA and total RNA were extracted by GeneRead DNA FFPE kit and by miRNeasy FFPE kit (Qiagen, Hilden, Germany), respectively, using the QIAcube-automated purification system. Extracted nucleic acids were quantified by Qubit 4 Fluorometer using Qubit Assay Kits (Thermo Fisher Scientific, Waltham, MA, USA), and quality was assessed by TapeStation 4200 (Agilent Technologies, Santa Clara, CA, USA) using Agilent ScreenTape Assays. Only patients with adequate DNA quantity (total extracted DNA >200 ng) from both tumor and matched NT were processed for DNA sequencing.

Tumor/NT pairs derived from 32 patients underwent library construction. DNA was fragmented by a Covaris M220 sonicator, and libraries were prepared using Illumina TruSeq Exome Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. A total of 60 samples, corresponding to 30 patient-matched tumor/NT pairs, passed library quality control on TapeStation 4200 (Agilent Technologies) and were submitted to library pooling and sequencing on Illumina NextSeq500 System according to the manufacturer’s standard protocol.

DNA sequencing data were processed as previously reported (40). Briefly, raw fastq files were quality-controlled with FastQC (41) and paired-end reads were aligned to the reference human genome (hg19) using a Burrows–Wheeler Aligner (BWAMEM, v0.7.12) (42). Duplicate and unmapped reads were identified and removed with Picard software (http://broadinstitute.github.io/picard/) and SAMtools v1.3.1.31. Reads were then post-processed according to Genome Analysis Toolkit (GATK) Best Practices 3.7 which include left alignment of small insertions and deletions, indel realignment, and base quality score recalibration.

Somatic single-nucleotide variants (SNVs) and small indels were called by two different variant callers: MuTect2 (v3.7) and Strelka (v.29.10). To create a high-confidence variant list, only the variants called by both algorithms were considered. Somatic variants (substitutions and indels) were annotated with Oncotator (v1.9.9.0). To remove false positives and polymorphisms, variants were then excluded based on at least one of the following additional filters: (i) read depth <30 in both tumor and NT; (ii) unidirectional call; (iii) alternative allele present in matched NT if the variant was not listed in the Catalogue Of Somatic Mutations In Cancer (COSMIC) database; (iv) C>A/G>T variants with a frequency <0.1 (oxoG artifacts); and (v) variants annotated in polymorphism databases (ExAC (43), NCBI dbSNP (44), and 1000 Genome Project (45)) without a COSMIC annotation. COSMIC Cancer Gene Census (https://cancer.sanger.ac.uk/census) was interrogated to explore the impact of somatic mutations in selected genes.

Gene fusions were assessed on total RNA by Oncomine™ Comprehensive Assay Plus RNA panel (OCA Plus RNA, Thermo Fisher Scientific) that covers more than 1,300 isoforms across 49 known cancer-related fusion drivers. Libraries were prepared using the Ion AmpliSeq Library Kit plus with OCA RNA plus pools and sequenced on an Ion GeneStudio S5 Prime sequencer using Ion 530 chips, Ion 510 & Ion 520 & Ion 530 Kit-Chef, and Ion Chef System (Thermo Fisher Scientific), according to the manufacturer’s instructions. Data were processed by Torrent Suite™ and analyzed by Ion Reporter™ software (5.18 version) with the “Oncomine Comprehensive Plus - w2.2 - Fusions - Single Sample” workflow.

Gene expression profiles were established by microarrays using Clariom S Pico Assay (Thermo Fisher Scientific). Total RNA was reverse transcribed, amplified, fragmented, biotin-labeled, and hybridized to Affymetrix GeneChip Human Clariom S (CLS) Arrays according to the manufacturer’s standard protocols. Washing and staining procedures were performed using the GeneChip Hybridization, Wash and Stain Kit (Thermo Fisher Scientific) on Affymetrix GeneChip Fluidics Station 450. Microarrays were scanned with the GeneChip Scanner 3000 7G system (Thermo Fisher Scientific), and data were obtained using Affymetrix GeneChip Command Console (AGCC) software.

Data were processed using the robust multi-array average (RMA) algorithm on paraffin samples (46); raw Affymetrix CEL file data were background-noise-adjusted, normalized, and log2-transformed using the oligo package and RMA function. Probes were annotated with Bioconductor annotation package clariomshumantranscriptcluster.db, whereas probes not associated with gene symbols and control probes were filtered out. Multiple probes mapping to the same gene were collapsed using the collapseRows function with the “MaxMean” method of the WGCNA package (47). All analyses were performed using RStudio version 4.0.3. Microarray data were deposited in the ArrayExpress repository with the accession number E-MTAB-13222.

Transcriptional subtypes were defined using three gene signatures derived from TCGA (5); the complete gene lists are published (5, 33), and the corresponding expression scores were calculated as previously reported (37). Briefly, the BRAF-/RAS-like signaling gene set comprises 71 genes; 69 out of 71 genes were available on the used CLS array and assessed for score computation. The BRAF-RAS score (BRS) was calculated as reported (48). Negative BRS values were defined as BRAF-like subtype, whereas positive BRS values were defined as RAS-like subtype as reported (5); close-to-zero BRS values were considered as borderline subtype. Thyroid differentiation (TD) and MAPK output gene sets comprise 16 and 52 genes, respectively. TD and MAPK output scores were calculated as mean of log2-transformed and median-centered expression across samples as previously reported (5, 33). Positive and negative TD score values were defined as high and low expressions of thyroid function genes, respectively. Positive and negative MAPK output score values were defined high and low MAPK pathway transcriptional activation, respectively.

The COSMIC database (https://cancer.sanger.ac.uk/cosmic, accessed on 30 January 2023) was interrogated to explore DICER1 (COSMIC gene COSG526495; transcript ENST00000526495.5) and ROS1 (COSMIC gene COSG418; transcript ENST00000368508.7) gene mutations. Thyroid tissue was specifically investigated, and the linked data for each gene were downloaded.

DICER1 and ROS1 mutations were reported in 24 and 8 different studies, respectively. We focused on the 19 and 6 studies reporting at least one mutated sample, respectively. For DICER1, an additional study derived from literature (11), and not included in COSMIC, was also considered.

Data derived from published studies were checked on the original publication to confirm mutation type, tissue histology, and patient-matched samples. Data derived from CGP (Cancer Genome Project) studies were manually curated: (i) samples from CGP study_542 were assigned to TCGA study (5) based on ID matching; (ii) one sample from CGP study_542 (COSMIC ID 2122053) with ROS1 synonymous mutation was excluded; and (iii) duplicated samples from CGP study_676 were excluded.

We investigated a retrospective series of TC patients collected at Fondazione IRCCS Istituto Nazionale dei Tumori of Milan. Caselist selection and processing are described in Material and Methods, Figure 1A and Supplementary Figure 1 . Starting from 126 FFPE tissues derived from 57 patients, and comprising TCs and patient-matched non-neoplastic thyroids (NTs), a final set of 30 tumor/NT pairs (collectively 60 tissues) passed quality control standards and were profiled first by DNA whole-exome sequencing (WES) and then by targeted RNA sequencing and transcriptomics ( Figure 1A ). Clinicopathological features are reported in Table 1 and Supplementary Table 1 . This series includes PTCs (90%) and PDTCs (10%). PTC histological types comprise classical and follicular subtypes, the less frequent tall cell and solid-trabecular subtypes, and two cases with mixed components ( Figure 1B , Table 1 and Supplementary Table 1 ), thus being representative of the histological heterogeneity observed in this tumor type.

DNA WES was established using patient-matched NT as filtering control for somatic mutation calling. We specifically focused on non-synonymous mutations, causing amino acid (aa) changing, and comprising missense, nonsense, and splice-site mutations with aa changing. The identified mutations for each tumor are in Supplementary Table 2 .

The mutation load of our set ( Figure 2A and Supplementary Figure 2 ) was low in agreement with that of thyroid cancers from TCGA ( Supplementary Figure 2A ) and from other TC series (5, 12, 13). The median number of mutations was 6.5 (range 1–28), with the majority of samples (19/30, 63%) harboring less than 10 mutations ( Supplementary Figure 2B ). The top mutated samples (number of mutations ≥15) included all the three tumors with PDTC component/histology ( Supplementary Figure 2B ), in agreement with the higher mutational burden reported in advanced and less differentiated TCs (12, 13). A significant correlation between patient age and mutation load was also observed ( Supplementary Figure 2C ), confirming a lower mutation load in younger patients ( Supplementary Figure 2D ) as previously reported (12).

As reported in the previous section, 11/30 cases did not display mutations in the abovementioned known/putative drivers. We therefore investigated these samples for the presence of gene fusions, well-established drivers in PTC, by a targeted RNA sequencing panel covering 49 cancer-related chromosomal rearrangements (see Material and Methods). In the analysis, we also included DICER1-mutated samples, as the co-occurrence of a DICER1 mutation with a rare gene fusion has been reported (14);and one BRAFV600E sample, included as negative control.

RET and NTRK3 fusions were found in four cases, all lacking mutations in BRAF, RAS, and DICER1 genes, thus confirming to be mutually exclusive with each other and with other known/putative drivers ( Figure 2B ). The identified fusions were RET/PTC1 (fusion partner CCDC6 gene) in three cases and ETV6/NTRK3 in one case, both representing the oncogenic fusions more frequently identified in PTC. RET/PTC1 fusions were found in the PTC classical subtype, whereas ETV6/NTRK3 was found in the follicular subtype; this agrees with the genotype/histological subtype distribution observed in other TC series (5, 11) and with the described association of NTRK fusions with the follicular growth pattern (54). Negative samples were one DICER1-mutated sample, the negative control (BRAFV600E), and three samples, which thus remained with the unknown gene driver. Five samples did not pass assay quality control due to low-quality RNA, thus displaying not available data for this analysis.

Collectively, considering also gene fusions, we identified gene alterations in 23/30 samples (77%), whereas 7 samples remained with unknown driver ( Figure 2E ).

To identify additional potential drivers in our cohort, we then revised somatic mutation data ( Supplementary Table 2 ).

In one patient, a ROS1 gene mutation was detected ( Figure 2C ). ROS1 codes for a receptor tyrosine kinase and is included in Cosmic Cancer Gene Census as containing mutations causally implicated in cancer. While its rearrangement has been found in various tumor types (55), including recent TC case reports (56, 57), its mutation seems a quite rare event in thyroid tumors. In the COSMIC database, 25 ROS1 mutations (2%) are described across 1,215 samples ( Supplementary Table 3 ). No common mutations are reported, and the ROS1_P2130A mutation detected in our patient (VAF 0.18, Supplementary Table 2 ) represents a new alteration for this gene in thyroid lesions. Other drivers (as BRAFV600E or HRAS, Supplementary Table 3 ) co-occur with ROS1 mutations, thus raising the possibility that this gene may not represent a standalone driver in TC. However, in our patient, mutations in known drivers are absent, thus suggesting that other alterations may be involved; due to unavailable data, we cannot confirm the absence of gene fusion (see previous section). Interestingly, in our patient ROS1 mutation co-occurs with ZFHX3_A472E ( Supplementary Table 2 ), and the same co-occurrence is observed in COSMIC ( Supplementary Table 3 ; COSMIC ID 2121935). The meaning and impact of these data remain to be established. Notably, ZFHX3 mutations co-occur with other gene drivers in 1.7% of PTC from TCGA (5).

In our cohort, the ROS1 mutation was found in a PTC with the rare solid-trabecular histology, whereas the few available COSMIC data ( Supplementary Table 3 ) describe ROS1 mutations in both follicular-pattern tumors (as FTC and FTA), as well as in ATC and BRAFV600E PTCs, thus suggesting that the current data are sparse and insufficient to address a genotype/phenotype association for ROS1 mutations in TC.

Revising the other samples with unknown drivers, especially the three with confirmed absence of gene fusions (i.e. sample ID HB08, HD97, HE27; Figures 2B, C ), we found a variable number of somatic mutations (range 1–25, Figure 2A ); mostly were missense mutations ( Supplementary Table 2 ) and affected genes previously reported in other cancer types, such as UBQLN3, PLCG2, and FANCF genes ( Figure 2C ).

In addition, in the HD97 patient, mutations in two genes encoding mitochondrial proteins were identified: the mitochondrial gene ND5 (NADH dehydrogenase subunit 5, p.F429L mutation) and the nuclear gene MTO1 (mitochondrial tRNA translation optimization 1, p.M386L mutation) ( Figure 2C and Supplementary Table 2 ). Of note, mutations in other mitochondrial genes were identified in six other patients of our caselist, and affecting ND4 and ND1 genes (NADH dehydrogenase subunit 4 and subunit 1), involved in the activity of the mitochondrial membrane NADH dehydrogenase (complex I), which catalyzes the electron transfer from NADH through the respiratory chain, and the COX1 gene (cytochrome c oxidase I) involved in the electron transport in mitochondrial respiratory chain complex III and IV ( Figure 2C ). Mutations in mitochondrial genes have been previously identified in thyroid cancer. Nonsense disruptive mitochondrial DNA mutations in complex I subunits have been reported in the oncocytic variant of thyroid carcinoma (58), and more recently in other TC subtypes (59–61). Interestingly, the ND4 gene mutation here identified (ND4_W24*, Supplementary Table 2 ) is nonsense, possibly affecting complex I formation. Furthermore, four of the patients with mitochondrial gene mutations are BRAFV600E positive, in agreement with a previously suggested correlation between BRAF mutation and mitochondrial alterations in TC (62).

Along with mitochondrial genes, other genes were found commonly mutated in both samples with unknown and known/putative driver ( Figure 2C ). These genes belong to different functional categories and cellular processes, such as cell proliferation (CDC27 gene), ERK signaling (HYDIN gene), transcription regulation by DNA binding (BCL11B and ZNF479 genes), extracellular matrix organization (OLFML2A gene), chromatin organization (MSL2 gene), calcium ion binding (HRNR gene), and apoptosis regulation (MTRNR2L8 and MTRNR2L1), or are less characterized genes ( Figure 2C and Supplementary Table 2 ).

To get more information about the pathways affected by the identified mutated genes, we performed a pathway-level analysis ( Supplementary Figure 3 ). The somatic mutations found in our cohort ( Supplementary Table 2 ) were tested on the Hallmark collection derived from Molecular Signature Database (63). We found that several mutated genes fall into biological processes related to (i) cell proliferation; (ii) p53 and apoptosis; (iii) stress responses as hypoxia and DNA damage; (iv) signaling pathways associated with KRAS, MTOR, TNFA, and estrogen receptor; (v) metabolic functions; and (vi) inflammation, consistently with the biological alterations and features typically observed in cancer (64). The vast majority of the genes identified in these pathways, however, are altered in samples with a known/putative driver, whereas very limited information is obtained for the unknown samples.

In our cohort, we then established transcriptional profiles by RNA microarray; two normal thyroids were included as control. We tested our samples with the three TCGA-derived gene signatures related to BRAF/RAS signaling (BRAF-/RAS-like subtypes), thyroid differentiation (TD score), and activation status of the MAPK pathway (MAPK output) (5). This was aimed not only to validate these transcriptional subtypes in samples with recognized drivers, but also, and more importantly, to obtain signaling information for the samples with putative and unknown drivers.

For NT controls and samples with established drivers, we confirmed literature findings. Normal thyroids displayed, as expected, a high TD score (5), low MAPK output, and RAS-like subtype as previously reported (34, 36) ( Figure 2F ). BRAFV600E samples were BRAF-like, with a low TD score and high MAPK output (4) ( Figures 2D, F ). Samples with RET/PTC1 fusion showed results similar to BRAFV600E samples, except for the higher TD score ( Figure 2F ) indicative of a partial preservation of thyroid function, as already reported (5, 35). On the contrary, RAS-mutated samples were confirmed as RAS-like, with higher TD score and intermediate MAPK output. Also, DICER1-mutated samples were RAS-like, consistently with previous reports (5, 14), and with the histological subtypes (follicular and PDTC solid-insular) in which these alterations were found. Of note, they displayed on average TD scores lower than those of RAS and higher than those of BRAFV600E samples, and reduced MAPK output, consistently with previous data (65). The ETV6/NTRK3 sample displayed a borderline RAS-/BRAF-like signaling subtypes, as already reported (5, 11).

Regarding the samples with unknown driver, the one with ROS1 mutation displayed a borderline RAS-/BRAF-like subtype, whereas the remaining samples were RAS-like, consistently with their histology (mostly follicular and solid/insular/trabecular) and displayed intermediate MAPK output and a heterogeneous degree of thyroid differentiation.