Lymph node metastasis‐dependent molecular classification in papillary thyroid carcinoma defines aggressive metastatic outgrowth

Lymph node metastasis (LNM) is the most important prognostic factor and a crucial indicator in the development of treatment strategies for patients with papillary thyroid carcinoma (PTC).1, 2 A vast majority (25%–60%) of patients with PTC undergo thyroidectomy with neck node dissection.3 However, researchers continue to debate whether LNM is leading to over-treatment of patients with PTC.1, 4 Interestingly, we have noted that upon stratifying LNM risk in PTC and clinical outcomes according to distinct molecular characteristics, LNM showing enrichment of genes related to inflammation and epithelial–mesenchymal transition (EMT) were more aggressive than metabolically adapted LNM, and this was applicable in other tumours. A flowchart of this study is displayed in Figure S1. For discovery, we ran differential gene expression analysis based on the negative binomial distribution (DESEQ2) for 97 patients without LNM [LNM (−)] and 195 patients with LNM [LNM (+)] from our PTC patient cohort (Table S1). One hundred four differentially expressed genes (DEGs) were upregulated in LNM (+), whereas the other 140 DEGs were downregulated (Figure 1A). In LNM (+) specimens, DAVID software revealed significant clustering of up-regulated genes related with EMT pathways (Figure S2A–C). Subsequent analysis of DEGs using the K-means unsupervised clustering algorithm was performed to define marked molecular subtypes, wherein K = 2 showed the greatest classification outcome with a cophenetic coefficient value of.736 (Figure 1B, Figures S3 and S4). To develop a predictive scoring model, we employed LASSO regression and 18 gene signatures that effectively deciphered molecular subtypes with an area under curve value of .97 (Figure 1C–E, Figures S5A–C and S6A,B). We then characterised the molecular characteristics of differently defined clusters by evaluating hallmark gene sets of a molecular signature database using single sample gene set enrichment analysis (ssGSEA).5 PTC patients from Cluster 1 in both our cohort and TCGA shared up-regulated patterns of metabolism-related hallmark pathways. PTC patients from Cluster 2 in both our cohort and TCGA exhibited common enrichment of EMT and immune response (Figure 2A, Figure S7). These differences in enrichment were not noted when patient samples were compared according to LNM status (Figure S8A,B). We noticed a dominant proportion of LNM (+) patients (133 [80.6%]) in Cluster 2 in our cohort. Remarkably, 80 N1b (cancer with lateral neck node metastasis) patient samples (48.5%) were allocated in Cluster 2, while 32 (25.2%) N1a (cancer with central neck node metastasis) and 30 (23.6%) N1b patient samples were in Cluster 1. The odds ratio to observe LNM in Cluster 2 was 4.098 (Figure 2B). Validation with an odds ratio of 4.691 was observed for LNM in Cluster 2, which indicates a strong clinical association between Cluster 2 and LNM (Figure S9). We then set out to determine the prognostic significance of our signature genes related with disease progression or recurrence by integrating clinical survival data from the TCGA database. In doing so, we noted that patients with LNM and high LASSO scores in Cluster 2 (upper two thirds of the designated cluster) displayed significant differences in their prognosis, compared to those with low LASSO scores (lower two-thirds of the designated cluster) among LNM patients in Cluster 1 (Figure 2C). Interestingly, aggressive outcomes were also noticed in N1a patients from Cluster 2 (Figure 2D). This is a crucial issue, as it is not clear whether pathological N1a metastasis can serve as a risk factor in patients with PTC.6 Furthermore, we evaluated thyroid differentiation scores and other important PTC-related oncogenic pathways and found corroborating evidence of biological aggressiveness in Cluster 2 (Figure 2E–L, Table S2). Meanwhile, progression-free survival according to LNM status or the presence of BRAF mutation, a unique genomic marker of PTC, did not show prognostic significance (Figure S10A,B). Additionally, we compared demographic and genomic information from PTC patients between different molecular clusters, as well as multivariate linear regression with LASSO score, to seek for any potential clinical relevance (Tables 1 and S3). Notably, BRAFV600E mutation, T stage, degree of LNM, and extrathyroidal extension (ETE) differed significantly between two clusters, and depth of ETE exhibited a significant proportionate relationship with LASSO score, suggesting a greater chance of severe metastasis in Cluster 2 PTC patients (Tables S4–S6). Our gene expression program was extended to actual LNM samples of PTC by applying it to public databases, GSE60542 and GSE151179 (Figure 3A, S11A–F).7, 8 Also, we challenged the applicability of our gene expression program in LNM of other tumour types and incorporated datasets of breast cancer and melanoma from TCGA, GSE56493, and GSE65904.9, 10 Likewise, molecular subtypes were deciphered using our proposed LASSO gene expression program, and identified molecular characteristics were validated in both primary tumour and LNM (Figure 3A). Moreover, LASSO scores of paired primary tumour and LNM samples showed significant relevance, suggesting preservation of molecular profiles after disease progression (Figure S12A–E). Analogous tendencies were observed in poorly differentiated thyroid carcinoma and anaplastic thyroid cancer patients in GSE76039 (Figure S13A,B). To expand the program's clinical significance, we analysed treatment responses to radioiodine therapy among PTC patients in GSE151179. Notably, significant radioactive iodine (RAI) uptake was found in Cluster 1 with markedly lower LASSO scores than patients without radioiodine uptake (Figure 3B,C). In addition, LASSO scores were higher in patients who underwent RAI treatment (RAIT) (Figure 3D) and in an intermediate- or high-dose group, compared to no treatment or low-dose group (Figure S14A, S14B). LASSO scores also showed a clear tendency to increase as suppressed Tg increased among all patients who underwent RAIT (Figure 3E) and patients who underwent high-dose RAIT (Figure S14C). Finally, we observed concordance between the aggressive outcomes noticed in Cluster 2 primary tumours with LNM and LNM samples from other tumour types (Figure 3F–K). Herein, we report two potential distinct LNM molecular subtypes of PTC with different clinical outcomes in terms of LNM and biological features. Proposed signatures were conserved throughout disease progression and validated in other thyroid cancer data and different tumour types. An EMT/inflammation axis, which we suggest as an aggressive LNM mechanism, requires support from additional studies for potential application in clinical settings. We thank Ji Young Kim, Hee Chang Yu, and Hoyoung Kim for their excellent technical support. The authors declare that they have no competing interests. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.