This study included patients with PTMC (14 N0 and 13 N1b) who underwent initial thyroid surgery at the Pusan National University Hospital, Busan, Korea. The biospecimens and data used for this study were provided by the Biobank of the Pusan National University Hospital, a member of the Korea Biobank Network. Tumor and the matched normal fresh frozen tissues were collected through several steps by pathologists. In detail, after gross examination, non-necrotic portions were excised from the resected tumor specimens by pathologists and were immersed in liquid nitrogen immediately. Matched normal thyroid tissue was obtained 2 cm away from the tumor. Sample adequacy for tumor content was assessed at the time of procurement and again on processed blocks by an experienced pathologist (K.U.C.). All tumor samples used in this study passed the purity assessment, and the average tumor cell content was 65%. Written informed consent was obtained from all participants prior to the survey, and the protocol of this study was approved by the Institutional Review Board of the Pusan National University Hospital, Busan, Korea (No. H-1804-021-065).

DNA was extracted from fresh frozen tissue samples using the QuickGene DNA tissue kit S with a QG-Mini80 instrument (Kurabo, Osaka, Japan). SureSelect sequencing libraries were prepared according to the manufacturer’s instructions (Agilent SureSelect^XT Human All Exon V5, Agilent, Santa Clara, CA, USA) by using the Bravo automated liquid handler. Sequencing was performed using an Illumina Novaseq 6000 system following the provided protocols for 2×100 sequencing at DNA Link Inc. (Seoul, Korea).

Quality assessment for the generated WES was performed using FastQC version 0.11.5 [14]. Burrows-Wheeler Aligner version 0.7.12 [15] was used for mapping and alignment based on Genome Reference Consortium Human Build 37 (GRCh37). Next, GATK version 3.5 was used for realignment and recalibration [16]. Variant calling for single nucleotide variants and indels was performed using MuTect version 1.1.7 [17] and Strelka version 1.0.14 [18]. Input sequence data from matched pairs for tumor and normal DNA were analyzed using MuTect/Strelka and, after removal of low-quality reads, the presence of variants in addition to expected random sequencing errors was determined [17,19]. Statistical analysis was performed to identify sites that may have a high probability for the occurrence of somatic mutations by calculating a log-odds score with high confidence. The score was compared to a cutoff determined using the log ratio of prior probabilities of the events under consideration. To reduce the risk of identification of false positives and remove artifacts, six different filters were used to screen candidate variant sites: (1) proximal gap; (2) poor mapping; (3) triallelic sites; (4) strand bias; (5) clustered positions; and (6) presence in the controls [17]. Somatic mutations were identified by comparing the sequencing data for matched tumor and normal samples, and we filtered germline mutations using the single nucleotide polymorphism (SNP) database build 134 and 1,000 Genome data. To reduce the occurrence of false positives, we filtered variants showing <5% variant allele frequency and with total coverage under 15×. Furthermore, missense SNVs and loss of function variants were used for subsequent analyses. The identified variants were annotated by Oncotator v1.8.0 [20] for functional information. Information on variant pathogenicity was annotated using ClinVar and Combined Annotation-Dependent Depletion (CADD) [21]. Variants with CADD scores over 30, which are predicted to be the 0.1% most deleterious substitutions possible in the human genome [21], were defined as pathogenic mutations.

Statistical analyses were performed by using the R program version 3.5.1 (R Foundation for Statistical Computing, Vienna, Austria). Continuous variables are presented as medians with interquartile range (IQR) and were analyzed by Student’s t test. Categorical variables are presented as frequencies and expressed as percentages, and they were analyzed by chi-square test. The MafCompare function in the Maftools package was used to determine associations between differentially mutated genes and lateral neck disease, wherein the mutation load for each gene was compared by Fisher’s exact test [22]. To compare the somatic mutation profiles of papillary thyroid carcinoma (PTC) using The Cancer Genome Atlas (TCGA), we obtained WES data including information on primary tumor sizes and LN metastasis from 291 patients [23]. P<0.05 was considered statistically significant.

The baseline characteristics of the 14 N0 and 13 N1b PTMCs are listed in Table 1. All patients had classical PTC. The median age of the patients at the time of initial thyroidectomy was 50 years (IQR, 46 to 61), and 23 (85%) patients were female. The median primary tumor size was 0.9 cm (IQR, 0.7 to 1.0), with no significant differences between the N0 and the N1b group (P=0.871). Microscopic and gross extra-thyroidal extension (ETE) were identified in 16 (59%) and two (7%) patients, respectively. Six (22%) patients had multifocal PTC, and lymphovascular invasion was present in two (7%) patients.

All patients underwent prophylactic central compartment node dissection, whereas therapeutic central and lateral neck dissection was performed in patients with clinically apparent lateral cervical LN metastasis. The median number of removed LNs was 4 (IQR, 2 to 8) in N0 and 35 (IQR, 32 to 50) in N1b PTMCs. Among the removed LNs, the median number of metastatic LNs was 7 (IQR, 3 to 15), with a maximum size of 12 mm (IQR, 7 to 12). Among patients with N1b, extranodal extension of metastatic foci from the metastatic LNs was detected in five cases.

We sequenced 27 matched DNA sample pairs from primary tumors and normal thyroid tissues. WES was performed to 100× coverage at the raw data level, and a median target coverage depth of 67× (range, 53× to 91×) in primary tumors and 70× (range, 50× to 99×) in matched normal thyroid tissues was followed (Supplemental Table S1). After calling variants with MuTect and Strelka, with highly sensitive algorithms and specific cutting-edge mutation calling methods, we identified 313 and 281 significant somatic variants in N0 and N1b tumors, respectively. The mutation burden was not significantly different between N0 and N1b tumors, with a median of 23 (range, 12 to 46) and 24 mutations (range, 12 to 50) in N0 and N1b PTMCs, respectively (P=0.918) (Supplemental Fig. S1). The top 20 genomic alterations across cases are shown in Fig. 1A, Supplemental Table S2. The most frequently mutated genes were DPP6, and PGS1, each occurring in four cases. Mutations in PGS1 SLC4A8, DAAM2 and HELZ were detected only in N1b disease, and the K158Q mutation in PGS1 was significantly enriched in N1b PTMCs (four cases; Fisher’s exact test P=0.041). BRAF mutations, all represented by the V600E mutation, were detected in three PTMCs (11%), i.e., two N0 (14%) tumors and one N1b (8%) tumor. Somatic alterations in genes included those that are a part of the receptor tyrosine kinase (RTK)-RAS signaling pathway; mutations in BRAF, CBLC, IRS1, NF1, and NTRK2 in N0 and BRAF, MAPK1, and RASGRP2 in N1b were most frequently detected in patients with PTMC (Fig. 2). The RTK-RAS and NOTCH pathways were impacted in three N1b patients. In the TCGA dataset, the BRAF V600E mutation was the most frequent alteration (54%), and was detected more frequently in N1b (66%) than N0 PTCs (49%, P=0.002) (Fig. 1B). The NRAS mutation was identified in 11% of patients with PTC, and the frequency of detection was higher in N0 (13%) than N1b tumors (3%, P=0.013). The BRAF V600E mutation was confirmed via polymerase chain reaction (PCR)-based Sanger sequencing (Supplemental Fig. S2).

To identify mutations implicated in the aggressive PTMC phenotype, we evaluated pathogenic mutations using ClinVar and the CADD scores (>30). We identified a median of 3 (range, 1 to 11) pathogenic mutations per case among N0 PTMCs and 4 (range, 2 to 7) among N1b PTMCs (Supplemental Tables S3, S4). Among them, BRAF, USH2A, CFTR, and PHIP mutations were classified as pathogenic based on both ClinVar and CADD scores (Table 2). In N0 PTMC, the BRAF V600E mutation (two cases) and a nonsense mutation in USH2A (one case) were detected. In N1b PTMC, the V600E mutation in BRAF (one case), a missense mutation in CFTR (one case), and nonsense mutations in USH2A (one case) and PHIP (one case) were detected. The CFTR V561E mutated N1b PTMC showed microscopic ETE with histological features of Hashimoto’s thyroiditis. Another single N1b PTMC harbored a nonsense mutation in PHIP (p.Y235*) and showed microscopic ETE and lymphovascular invasion. Both tumors were multifocal PTMCs. In the TCGA dataset, one N0 PTC harbored a missense mutation in CFTR (p.L702I) and there were no tumors that showed mutations in PHIP.