This study was approved by the Institutional Research and Ethics Committee at Chungnam National University Hospital (CNUH-2020-11-004-001), and Institutional Research and Ethics Committee at Seoul National University Hospital (H-1508-147-700). Forty PTC patients were enrolled from Seoul National University Hospital, and the other patients were enrolled at Chungnam National University Hospital. Informed consent for the study was obtained from the patients.

Fresh frozen tissues from 282 PTC patients who underwent thyroidectomy were analyzed using a massively parallel sequencing method. Normal thyroid tissues (n=155) were obtained from patients who underwent thyroidectomy, with confirmed cases of differentiated thyroid cancer. The clinicopathological characteristics of patients according to histology are shown in Supplementary Table 1.

In addition, a 2,100 Agilent Bioanalyzer (Agilent Technologies, Waldbronn, Germany) was used to estimate the RNA integrity number score. After stranded total RNA-seq library construction, we used the TruSeq Stranded Total RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) based on the manufacturer’s instructions. We sequenced the samples via the HiSeq 4000 sequencing system, yielding 2×100-bp sequencing reads. The raw data were deposited in the Korean Nucleotide Archive (KoNA, https://kobic.re.kr/kona) with the accession ID PRJKA210106.

Sequencing reads were trimmed by trim galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) with a threshold of average sequence quality over 30 and mapped to the human reference genome GRCh38 by STAR aligner (STAR-2.7.9a) [18]. Ensembl gencode (ver. 27) gene annotation was used to profile gene expression from raw read counts, normalization, and quantification by cufflinks at the gene level [19].

We performed single-sample gene set enrichment analysis (ssG-SEA) for functional characterization by samples. For the reference gene matrix transposed (GMT) file, we downloaded GMT files in MSigDB [20]. Then, the gene set variation analysis (GSVA) function in the GSVA package was used with the default option [21]. We used Enrichr for functional classification called GSEA by differentially expressed genes (DEGs) with a threshold P-value below 0.01 from Fisher’s exact test [22].

Sequenced reads were mapped to the human reference genome GRCh38 for fusion detection using STAR fusion with default options [5]. The result of fusion analysis was filtered out through the following three processes: (1) genes were eliminated, such as immunoglobulin genes, uncharacterized genes, and RNA genes; (2) fusions including the same gene or paralog gene were removed; and (3) fusions detected in normal samples were omitted. We performed outlier analysis for fusion candidates. The outlier means that the sample’s level at a given gene was greater than the 75th percentile +1.5 interquartile range. This trend should not be discovered in the normal group.

We used the thyroid differentiation score (TDS) score from the Thyroid Carcinoma cohort from the Cancer Genome Atlas (TCGA-THCA) study [9]. TDS scores is from ssGSEA performed with 16 genes (DIO1, DIO2, DUOX1, DUOX2, FOXE1, GLIS3, NKX2-1, PAX8, SLC26A4, SLA5A5, SLC5A8, TG, THRA, THRB, TPO, TSHR) associated with thyroid function and metabolism from the TCGA study. As a result of calculating the correlation of the TDS scores, correlation values were highly significant (R=0.859, P<2.2e-115).

We carried out immune-subtyping from xCell (https://xcell.ucsf.edu/) [6]. The abundances of 64 various cell types for PTC samples were computed using fragment per kilobase of transcript per million (FPKM). We divided the patients into three major clusters using unsupervised K-means clustering using pheatmap package from R 4.1.2 (https://cran.r-project.org/). These groups have a significant difference of ImmuneScore from xCell. Thus, we named them immune-hot, -intermediate, and -cold. ImmuneScore is a sum of 10 immune cell types such as B-cells, CD4+ T-cells, CD8+ T-cells, dendritic cells, eosinophils, macrophages, monocytes, mast cells, neutrophils, and NK cells. It is an experimental score and known to better optimize tumor microenvironment. In three subtypes, we compared various immune abundances and RNA expression (Fisher’s exact test, P<0.01). To discover sets of coexpressed and hub genes in the network, we utilized CEMiTool [23]. This tool integrates the identification and analysis of coexpression gene modules, evaluating which modules contain genes that are overrepresented by specific geneset terms with P-values below 0.01 from Fisher’s exact test.

We collected 282 PTC samples (Supplementary Table 1). The clinicopathological characteristics included capsular invasion, extrathyroidal extension (ETE), central and lateral lymph node metastasis (C_LNM and L_LNM), lymphovascular invasion (LVI), recurrence, age, and sex. The average age of the PTC patients was 49 years. The mean tumor size of the PTC samples was 1.83 cm.

We quantified FPKM values for 15,690 genes through a standard RNA-sequencing analysis pipeline. First, we identified DEGs in tumor samples compared to normal samples, with a fold-change >1.5 and P<0.01, and we selected a total of 2,307 genes as DEGs (Fig. 1A). There were 1,039 downregulated genes and 1,268 upregulated genes. The top five highly upregulated genes in PTC were ZCCHC12, FN1, CHI3L1, SLC34A2, and DCSTAMP.

We performed gene set analyses with DEGs using the Pathways-KEGG 2021 Human and MSigDB Hallmark gene sets in Enrichr (Fig. 1B). The upregulated gene set terms in PTC versus normal samples included cytokine receptor interactions, p53 signaling, and complement and coagulation in KEGG, and the epithelial-mesenchymal transition (EMT), the inflammatory response, and allograft rejection in Hallmark. The down-regulated gene set terms included thyroid hormone synthesis and amino acid metabolism in KEGG and the estrogen response, ultraviolet response, and hypoxia in Hallmark. Interestingly, many immuneassociated gene-sets were significantly enriched, such as the EMT (including ECM1, SPARC, and LOXL1) and the inflammatory response (including ICAM4, MMP14, and TLR2).

We also found that immune-related signaling was upregulated in PTC. When patients with and without LNM were compared, immune-related signaling (e.g., nuclear factor kappa-light-chain-enhancer of activated B cells [NF-κB], tumor necrosis factor-alpha [TNF-α], the epithelial cell score, Wnt signaling, CD4+ T-cells, and memory B-cells) was exclusively higher in patients with LNM. In addition, patients with recurrence had higher scores for CD8+ T-cells and Th1 cells than patients without recurrence (Fig. 1C).

We delineated pathway maps showing the expression of all genes with KEGG terms (Supplementary Fig. 1). In PTC tumors, we significant upregulation of many components of the cytokine receptor interaction pathway, including CCR1-2, CCR4-5, and CCR7 in the CC-subfamily and TGFBR1, BMPR2, and ACVR1 in the transforming growth factor-beta family.

Attempts to apply immunotherapy to thyroid cancer are continuing, and in addition to the discovery of many immune checkpoint inhibitors (ICIs) such as PD-L1, CTLA4, IDO1, and HAVCR2, various immune signaling pathways, such as TNF-α signaling and the IL-6/JAT/STAT3 pathway have been discovered [24,25]. To investigate the role of immune-escape signaling in advanced PTC, we investigated the expression of genes related to immunotherapy and immune subtyping using ImmuneScore. We performed an immune signature analysis using xCell from RNA sequencing data in PTC (Methods). We divided PTC into three groups using K-means clustering with 64 immune-related scores. Based on the immune scores, we classified patients into immune-hot (n=35), immune-intermediate (n=107) and immune-cold (n=140) groups. We identified significantly different immune-related scores and expression of the top 50 genes across the three immune groups (Fig. 2A). In the immune-hot group, B-cells, CD4+ or CD8+ T-cells, dendritic cells, and ImmuneScore values were higher. In the immune-intermediate group, astrocytes, fibroblasts, hepatocytes, and preadipocytes were enriched. In the immune-cold group, the TDS, common lymphoid progenitors, myocytes, osteoblasts were elevated (Fig. 2A, top panel). Analyses of the associations between clinical data and the immune-subtype showed significant differences in ETE status, C_LNM, and L_LNM (Fig. 2A, middle panel, Supplementary Table 2). Interestingly, the immune-hot group revealed poor prognostic indicators (tumor size, ETE, c_LNM, and L_LNM). We selected significantly different genes from each group and performed gene-set enrichment analysis. In the immune-hot group, genes related to interferon-gamma, such as ISG20, CASP8, and PTPN6, showed significant differences (P<0.001). In the immune-intermediate group, genes related to the EMT (e.g., LAMA2, MMP2, VEGFA, and BMP1; P<0.001) and in the immune-cold group, genes related to Myc (e.g., XRCC6, XPOT, NPM1, and CDK4; P<0.001) were significantly enriched (Fig. 2A, bottom panel). An analysis of the distribution of the TDS across the three immune groups showed that the TDS of the immune-cold group was significantly higher and the TDS of the immune-hot group was significantly lower than those of the other groups (Fig. 2B). In addition, we extracted the distribution of the TDS across immune-subtypes and genes that showed high correlations with the TDS (Fig. 2C). As the TDS decreased, the proportion of immune-hot tumors increased and the proportion of immune-cold tumors decreased. Conversely, as the TDS score increased, the proportion of immune-hot tumors decreased and the proportion of immune-cold tumors increased. Moreover, genes such as B3GNT9 (r=0.80), TSHR (r=0.77), MMP15 (r=0.77), LIPG (r=0.76), and CD24 (r=0.76) showed very high correlations with the TDS.

In the immune-hot group with a low TDS score, the expression of genes related to immune-escape signaling (e.g., CTLA4, IDO1, LAG3, and PDCD1) was significantly higher than in the other groups (Fig. 2D). Finally, we performed module and co-expression analyses across the three immune groups (Fig. 2E). In the immune-hot group, genes such as ITGA4, UBC, NFKB1, and PRKCA belonging to the interferon-gamma module acted as hub genes with many interactions, in the immune-intermediate group, MMP2 was a hub gene, and in the immune-cold group, XPOT, TCEA1, and IARS were identified as hub genes. These data suggest that immune subtypes, including immune-hot, -intermediate, and -cold groups, were significantly associated with TDS, and prognostic clinical parameters were associated with the immune-hot signature.

We performed fusion analysis from RNA-sequencing data with STAR fusion (Methods), identified 3,160 fusions, filtered out low-significance fusions, and finally obtained 1,129 fusions. All significant fusions were found in PTC samples. Each fusion count is shown in order of frequency (Fig. 3A). FBXO25 fusion was found in 22 samples, followed by RET, CA13, NTRK1, and MET fusions. We performed fusion outlier analysis with three fusion genes (FBXO25, RET, and CA13) to determine driver fusions (Methods). No significant upregulation was observed for CA13 and FBXO25 fusion samples (Supplementary Fig. 2A). The expression of RET in all fusion-positive tumor samples was higher than that in fusion-negative samples, and CCDC6, NCOA4, and DLG5 were partner genes of RET fusion (Fig. 3B). We investigated how the mitogen-activated protein kinase (MAPK) and phosphatidylinositol 3-kinase (PI3K) pathways, which are known to be downstream signaling pathways impacted by RET fusion, are differentially regulated [26-28] by different fusion partners of RET. DLG5 was excluded from further analysis, as it was found in only one sample. We divided samples into RET-CCDC6, RET-NCOA4, and RET-WT and selected DEGs (Fisher exact test, P<0.05). In the RET-CCDC6 and RET-NCOA4 fusions, 213 and 1,260 genes were significant, respectively, and among them, genes belonging to the MAPK and PI3K signaling pathways were selected and described (Fig. 3C). In the RET-CCDC6 fusion group, CREB1 and PIK3R1 were upregulated, and PPP2R2D and MAP2K7 were downregulated. In the RET-NCOA4 group, more genes belonging to the PI3K and MAPK signaling pathways were differentially expressed than in the RET-CCDC6 fusion group. When the association with clinical factors according to the two partner genes was analyzed, ETE and L_LNM were significantly more frequent in RET-NCOA4 patients than in the RET-CCDC6 group (Supplementary Fig. 2B). We then selected marker genes that can distinguish RET-CCDC6 and RET-NCOA4, and discovered FGF13 in RET-CCDC6 fusion and AKT1 and CCND2 in RET-NCOA4 fusion (Supplementary Fig. 3, Fig. 3D). The three genes were validated in the TCGA-THCA cohort (Fig. 3E). Thus, we found that RET fusion differentially regulated the PI3K and MAPK signaling pathways according to the partner genes of RET (Fig. 3F) [29].

One of the unmet medical needs in treating thyroid cancer is the identification of biomarkers that predict recurrence or distant metastasis [30]. Clinical factors such as age, sex, tumor size, extracapsular invasion, ETE, C_LNM and L_LNM, and laterallymphovascular invasion are currently used to predict recurrence or distant metastasis [31,32]. Using the Fisher exact test for each clinical factor, we found that recurrence occurred more frequently in patients with large tumors, ETE, and L_LNM (Fig. 4A). For these three significant risk factors, we divided samples into recurrence and non-recurrence groups. Most recurrence samples had all three risk factors (n=24, 59%), and most nonrecurrence samples did not have any of the three risk factors (n=88, 36%) (Fig. 4B). However, interestingly, we found three no-risk-factor samples in the recurrence group and 42 all-risk-factor samples in the non-recurrence group. We named these two groups No_RF and All_RF to distinguish them from other patterns of recurrence. We then performed statistical tests to identify genes that could predict recurrence in the No_RF group and identified the homeobox d9 (HOXD9) gene, which was found to have a significant relationship with disease-free survival in the TCGA database using GEPIA (http://gepia.cancer-pku.cn/) (Fig. 4C and D).

To understand the relationship between HOXD9 expression and recurrence, we divided the patients into HOXD9-high and HOXD9-low groups by the median value of HOXD9 expression, compared the ratio of recurrence, and confirmed that the HOXD9-high group had significantly more frequent recurrence than the HOXD9-low group (P<0.001) (Fig. 4D). Furthermore, the HOXD9-high group had larger tumors and a higher frequency of ETE than the HOXD9-low group (Supplementary Fig. 4A and B). We found that HOXD9 was significantly related to recurrence, in addition to other factors known to cause recurrence (Fig. 4E). Additionally, gene-set enrichment analysis revealed that the NF-κB signaling pathway was the most significant gene set among the selected genes at a threshold of P<0.05 from the Fisher exact test and a fold change >1.5 (Fig. 4F and Supplementary Fig. 4C).

We then analyzed transcription factors regulating the expression of HOXD9 using a transcription factor database (Methods) and validated the genes in the TCGA-THCA cohort (P<0.001) (Supplementary Fig. 4D). Transcription factor 3 (TCF3) and enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2) are known to bind to the promoter region of HOXD9, and the correlation with HOXD9 was also significant (Fig. 4G). Thus, we suggest that TCF3 and EZH2 may promote overexpression of HOXD9 associated with recurrence. Therefore, we suggest that HOXD9 may be a novel prognostic marker of recurrent PTC that acts through the NF-κB signaling pathway (Fig. 4H).