Multi-gene co-mutations of BRAF with TERT, PIK3CA, or TP53 are powerful predictors of central lymph node metastasis in papillary thyroid carcinoma

Background: The accurate preoperative prediction of lymph node metastasis (LNM) in papillary thyroid carcinoma (PTC) poses a significant clinical challenge. Although clinicopathological features are commonly utilized, their predictive accuracy remains limited, and the role of multi-gene co-mutations is not fully understood.

Objective: This study aimed to develop and validate an integrated risk model that combines next-generation sequencing (NGS) data with clinicopathologic features for the preoperative prediction of LNM in PTC.

Methods: We retrospectively analyzed 521 patients with PTC. Gene mutations were analyzed using NGS. Independent risk factors for central (CLNM) and lateral (LLNM) lymph node metastasis were identified through univariate and multivariate logistic regression analyses.

Results: The BRAF V600E mutation was the most prevalent (82.15%). Notably, high-risk multi-gene co-mutations —specifically, BRAF V600E co-occurring with TERT, PIK3CA, and/or TP53)—were identified as the strongest independent risk factor for CLNM (odds ratio [OR] = 6.319, 95% confidence interval [CI]: 1.738–22.976, P = 0.005). Other significant risk factors included male sex, age <45 years, bilateral lesions, tumor size >1 cm, lymphovascular invasion (LVI), and extrathyroidal extension,with gross ETE demonstrating the highest ORs (> 21).

Conclusion: Preoperative NGS profiling, particularly the detection of high-risk multi-gene co-mutations, provides a powerful tool for refined risk assessment. This molecularly guided strategy has the potential to inform personalized surgical planning directly, optimizing the extent of lymph node dissection to improve oncologic outcomes while minimizing unnecessary morbidity.

1 Introduction

Papillary thyroid carcinoma (PTC) is the most prevalent endocrine malignancy, accounting for approximately 80% of all thyroid cancer cases (1). While generally associated with a favorable prognosis, a significant clinical challenge is the high frequency of cervical lymph node metastasis (LNM) at diagnosis, particularly involving the central lymph node metastasis (CLNM) and lateral lymph node metastasis (LLNM) compartments (2). Accurate preoperative assessment of LNM status is crucial, as it directly influences the decision to perform and the extent of prophylactic lymph node dissection, thereby impacting both oncological outcomes and surgical-related morbidity (3, 4).

Currently, preoperative risk stratification primarily relies on ultrasonography and clinicopathological features such as tumor size and location (5, 6). However, the sensitivity and specificity of imaging for detecting CLNM, especially in occult metastases, remain suboptimal (7). Furthermore, while certain clinical factors are associated with higher risk, they lack the precision needed for truly individualized surgical planning.

In the era of precision medicine, molecular markers have emerged as promising tools for refining risk assessment. The BRAF V600E mutation, the most common genetic alteration in PTC, has been extensively studied for its association with more aggressive clinicopathological characteristics, including LNM (8, 9). Nevertheless, the predictive power of BRAF V600E as a standalone marker is often debated and appears insufficiently robust to guide surgical decisions alone (10, 11). This limitation highlights the complexity of tumorigenesis and suggests that the aggressive behavior of PTC may be driven not by a single gene, but by the synergistic effects of multiple genetic alterations. For instance, co-occurring mutations in genes such as TERT (promoting cellular immortality), PIK3CA (enhancing invasion), and TP53 (disabling apoptosis) are hypothesized to collaborate with BRAF, leading to a more metastatic phenotype (12, 13). Previous studies, such as those by Vuong et al. (14), have begun to highlight the prognostic superiority of dual mutations (e.g., BRAF+TERT) over single markers.

However, while recent studies have significantly advanced our understanding of specific dual mutations, such as the prognostic impact of concurrent BRAF V600E and TERT promoter mutations (14) and the role of BRAF/TP53 co-alterations (15), the synergistic effect of broader, predefined multi-gene constellations (specifically BRAF co-mutated with TERT, PIK3CA, and/or TP53), their differential impact on LNM patterns (CLNM vs. LLNM), and crucially, their quantifiable value when integrated into a preoperative clinical prediction model, remain underexplored.

Therefore, to address this gap, we conducted a retrospective study utilizing high-throughput next-generation sequencing (NGS) in a large, well-characterized Chinese cohort of 521 PTC patients. We aimed to systematically profile the mutation landscape, investigate the synergistic effects of these specific multi-gene co-mutations, and ultimately integrate these molecular signatures with conventional clinicopathological parameters to construct and validate a quantitative, preoperative prediction model (nomogram). The goal is to develop a clinically translatable, integrated (molecular-clinicopathological) tool that guides individualized surgical decision-making regarding the necessity and extent of lymph node dissection (10, 11).

2 Materials and methods

2.1 Materials

This was a single-center, retrospective, exploratory study. It aimed to investigate the feasibility of integrating multi-gene co-mutation patterns into preoperative risk assessment. It included all consecutive patients diagnosed with PTC by postoperative pathology at Deyang People’s Hospital between June 2022 and January 2025, who had undergone next-generation sequencing (NGS) testing (n=521).

The inclusion and exclusion criteria were as follows:

2.1.1 Inclusion criteria

(1) Postoperative histopathological confirmation of PTC; (2) Availability of formalin-fixed paraffin-embedded (FFPE) tumor tissue from the surgical specimen for molecular testing; (3) Completion of a standardized NGS assay.

2.1.2 Exclusion criteria

(1) History of any other malignant tumor; (2) Prior treatment for thyroid cancer (e.g., radiotherapy, chemotherapy, or targeted therapy); (3) Incomplete clinicopathological or molecular data.

NGS testing was performed postoperatively on FFPE tissue sections from the primary tumor as part of the standard diagnostic workflow. Clinicopathological features and NGS results were retrospectively analyzed.

The pathological diagnosis of PTC and evaluation of all histopathological features were independently reviewed and confirmed by two board-certified pathologists according to the World Health Organization (WHO) classification criteria. Surgical procedures included thyroid lobectomy + isthmusectomy/total thyroidectomy + central/lateral cervical lymph node dissection. The extent of cervical lymph node dissection was performed in accordance with the 2012 American Thyroid Association management guidelines (16).All patients were treatment-naïve at the time of surgery, with no prior history of radiotherapy, chemotherapy, or other targeted therapies before surgery. This study was approved by the Ethics Committee of Deyang People’s Hospital.

Among the 521 enrolled patients, the majority were female (72.4%) and under 45 years of age (65.5%). Most tumors were ≤1 cm in diameter (55.7%), unifocal (77.2%), and of the classical variant (95.6%). Lymphovascular invasion (LVI) was present in 8.6% of cases. The BRAF V600E mutation was prevalent (82.1%), while high-risk multi-gene co-mutations, defined as BRAF V600E co-mutated with TERT, PIK3CA, and/or TP53, were relatively rare (3.8%). We acknowledge that the limited sample size of this key subgroup constrains the statistical power for detailed interaction analyses.

2.2 Targeted sequencing and bioinformatic analysis

Targeted sequencing was performed using an in-house, laboratory-developed test (LDT) for thyroid cancer, which is based on a standardized commercial workflow. Genomic DNA was extracted from formalin-fixed, paraffin-embedded (FFPE) tumor samples using the QIAamp DNA FFPE Tissue Kit (Qiagen, Germany). Sequencing libraries were constructed using the Thyroid Cancer 88-Gene Targeted Sequencing Kit (RuiJing Biosciences, China). This panel targets the full exonic regions of 86 genes and specific intronic regions for the detection of fusions in 6 genes (88 genes in total; for the complete gene list see Supplementary Table 1), covering key driver genes such as BRAF, TERT, RAS, TP53, and PIK3CA. High-depth sequencing was conducted on an Illumina NextSeq 550 platform (Illumina, Canada) to achieve an average sequencing depth of >500×.

Bioinformatic analysis followed the vendor’s (RuiJing Biosciences) standard pipeline. Raw sequencing data underwent adapter trimming and quality filtering. Clean reads were aligned to the human reference genome (GRCh37/hg19), followed by duplicate read marking. Variant calling for single nucleotide variants and small insertions/deletions was performed following the Genome Analysis Toolkit (GATK) best practices. Identified variants were annotated against public databases (dbSNP, ClinVar, and COSMIC) and classified according to established clinical guidelines (e.g., AMP/ASCO/CAP).

Rigorous quality control was applied throughout the process. Key metrics assessed included tumor cell content (minimum 10%), DNA input quantity, library complexity, sequencing quality (e.g., Q30 score >80%), and on-target rate. Only samples passing all predefined QC thresholds were included in the downstream analysis.

2.3 Statistical methods

Statistical analysis was performed using SPSS 25.0 software. Categorical variables were expressed as frequencies (%). Differences between groups were assessed using univariate analysis (Chi-square or Fisher’s exact test, as appropriate). To ensure comprehensive and transparent reporting, all pre-specified clinicopathological and molecular variables were included in this univariate analysis, and their complete results are presented (in Tables 1, 2, 3), regardless of statistical significance.

Table 1

Table 1. Baseline clinicopathological characteristics of the study population (n=521).

Table 2

Table 2. Univariate analysis of risk factors for lymph node metastasis [n (%)].

Table 3

Table 3. Univariate analysis of central and lateral cervical LNM [n (%)].

Variables with a univariate P value < 0.05 were considered candidate predictors for subsequent multivariate binary logistic regression, which was performed using the enter method to identify independent risk factors. Before model fitting, multicollinearity among all candidate variables was assessed using the variance inflation factor (VIF); a VIF value ≥ 5 was considered indicative of substantial collinearity, prompting the retention of only the clinically more relevant variable. Results of the logistic regression are presented as odds ratios (ORs) with 95% confidence intervals (CIs).

The discriminative ability of the final multivariate model was evaluated using the concordance index (C-index), equivalent to the area under the receiver operating characteristic(ROC)curve. The 95% confidence interval for the C-index was calculated using the DeLong method. Model calibration was assessed by comparing predicted probabilities with observed outcomes. Internal validation and overfitting correction were performed using bootstrap resampling with 1000 repetitions, as detailed in the Results section. A two-sided P value < 0.05 was considered statistically significant.

3 Results

3.1 Baseline characteristics of the study population

This study included 521 patients with papillary thyroid carcinoma. Their baseline clinicopathological characteristics and the associations with central lymph node metastasis (CLNM) are summarized in Table 1. In univariable analysis, younger age (<45 years), male sex, larger tumor diameter (>1 cm), bilateral lesions, specific pathological subtypes, extrathyroidal extension (ETE), lymphovascular invasion (LVI), and high-risk co-mutations were all significantly associated with a higher prevalence of CLNM (all P < 0.05). In contrast, multifocality and the BRAF V600E mutation alone showed no significant association with CLNM status.

3.2 Distribution of gene mutations

Mutations were detected in 94.43% (492/521) of the specimens. Single point mutations or gene fusion mutations were detected in 88.48% (461/521), while multi-gene co-mutations (≥2) were detected in 5.95% (31/521). Among these, 3.84% (20/521) were BRAF combined with TERT/PIK3CA/TP53. No mutations were detected in 5.57% (29/521) of the specimens. The detailed spectrum of gene mutations is presented in Table 4.

Table 4

Table 4. Summary of gene mutation types in PTC patients detected by high-throughput sequencing.

3.3 Risk factors for CLNM/LLNM

3.3.1 Univariate analysis of risk factors for CLNM

CLNM was found in 54.70% (285/521) of patients. The metastasis rate was significantly higher in males, patients aged <45 years, tumor size >1 cm, bilateral lesions, high-risk co-mutations, tall cell/sclerosing subtypes, extrathyroidal extension, and Lymphovascular Invasion (LVI) (Table 2).

3.3.2 Univariate analysis of risk factors for LLNM

LLNM was found in 24.38% (127/521) of patients. The metastasis rate was significantly higher in males, patients aged < 45 years, tumor size >1 cm, multi-gene mutations, bilateral lesions, tall cell/sclerosing subtypes, extrathyroidal extension, and Lymphovascular Invasion (LVI) (Table 2).

3.3.3 Multivariate analysis of risk factors for CLNM/LLNM

Variables with P < 0.05 in the univariate analysis were included in the multivariate logistic regression model. Independent risk factors for CLNM (Table 5) included male sex, age <45 years, bilateral lesions, lymphovascular invasion (LVI), and extrathyroidal extension. High-risk multi-gene co-mutations (BRAF with TERT, PIK3CA, and/or TP53) were the strongest independent predictor (OR = 6.319, 95% CI: 1.738–22.976, P = 0.005). Independent risk factors for LLNM were: age <45 years, bilateral lesions, tumor size >1 cm, Lymphovascular Invasion (LVI), and extrathyroidal extension (Table 6). Multicollinearity diagnostics indicated that all variance inflation factors (VIFs) for the independent variables in the final models were close to 1 (range: 1.002–1.008), confirming the absence of substantial multicollinearity.

Table 5

Table 5. Multivariate analysis of CLNM.

Table 6

Table 6. Multivariate analysis of LLNM.

3.3.4 Development of the predictive nomogram

Based on the multivariate logistic regression analysis, we developed an integrated molecular-clinicopathological prediction nomogram to quantify the individual probability of CLNM (Figure 1). This tool, which we hereafter refer to as “the nomogram,” integrates the six identified independent predictors: sex (female as reference), age (≥45 years as reference), bilateral lesions, high-risk multi-gene co-mutations, lymphovascular invasion (LVI), and extrathyroidal extension (ETE).

Figure 1

Figure 1. Nomogram for preoperative prediction of central lymph node metastasis (CLNM) in papillary thyroid carcinoma. To use the nomogram, for each variable, locate the patient’s value on the corresponding axis and draw a line upward to the ‘Points’ axis to determine the score. Sum all points, locate the total on the ‘Total Points’ axis, and draw a line downward to the ‘Risk of CLNM’ axis to obtain the individualized probability.

We performed rigorous internal validation to evaluate the nomogram’s predictive performance. The model demonstrated acceptable discriminative ability, with a concordance index (C-index) of 0.712 (95% CI: 0.668–0.755),which is equivalent to the area under the receiver operating characteristic curve (AUC). The 95% confidence interval was calculated using the DeLong method. Internal validation using bootstrap resampling with 1000 repetitions further confirmed its robustness and calibration. The bootstrap-corrected calibration curve (Figure 2) showed good agreement between the nomogram-predicted probabilities and the actual observed outcomes, indicating satisfactory calibration of the model validated through the bootstrap resampling procedure.

Figure 2

Figure 2. Calibration curve for the CLNM prediction nomogram assessed by bootstrap resampling (n=1000 repetitions). The solid line represents the agreement between the nomogram-predicted probability and the actual observed frequency of CLNM after bootstrap correction. The dashed line represents the ideal reference. Close alignment indicates satisfactory model calibration.

3.4 Stratification power of high-risk co-mutation in microPTCs

Given that 55.7% (290/521) of the cohort consisted of microPTCs (tumor size ≤1 cm), and considering the evolving clinical paradigm towards active surveillance for such tumors, we specifically evaluated the predictive value of the high-risk multi-gene co-mutation within this subgroup.

Among the 290 patients with microPTCs, 11 (3.8%) harbored the high-risk co-mutation (BRAF/TERT/PIK3CA/TP53). The CLNM rate in these high-risk microPTCs was 81.8% (9/11), which was significantly higher than the rate of 46.2% (129/279) observed in microPTCs without the co-mutation. This difference corresponded to a >5-fold increase in the odds of CLNM (Odds Ratio [OR] = 5.23, 95% Confidence Interval [CI]: 1.21 - 22.66; P = 0.027, Fisher’s exact test).

This subgroup analysis confirms that the high-risk co-mutation serves as a powerful discriminator of biological aggressiveness even within microPTCs, highlighting its potential clinical utility for risk stratification in the context of active surveillance decision-making.

3.5 Risk factors for multiple lymph node metastases (>5 nodes) in central and lateral cervical regions

3.5.1 Univariate analysis of risk factors for multiple lymph node metastases (>5 nodes)

LNM occurred in 60.46% (315/521) of patients, of which 71.43% (225/315) had 1–5 nodes involved, and 28.57% (90/315) had >5 nodes involved. The rate of multiple lymph node metastases (>5 nodes) was significantly higher in patients with tumor size >1 cm, bilateral lesions, tall cell/columnar cell/sclerosing subtypes, extrathyroidal extension, and lymphovascular invasion (Table 3).

3.5.2 Multivariate analysis of risk factors for multiple lymph node metastases (>5 nodes)

Variables with P < 0.05 in the univariate analysis were included in the multivariate logistic regression model. Independent risk factors for multiple lymph node metastases (>5 nodes) were: tumor size >1 cm, bilateral lesions, and gross extrathyroidal extension (Table 7), consistent with previous studies (7).

Table 7

Table 7. Multivariate analysis of multiple LNM (>5 nodes) in central and lateral cervical regions.

4 Discussion

4.1 Summary of principal findings

This study establishes a comprehensive risk model for preoperative prediction of lymph node metastasis (LNM) in papillary thyroid carcinoma (PTC) by integrating next-generation sequencing (NGS) data with detailed clinicopathological characteristics. Our analysis unequivocally identifies high-risk multi-gene co-mutations (BRAF with TERT, PIK3CA, or TP53) —as the most powerful independent predictor for central LNM (CLNM) (OR = 6.319, 95% CI: 1.738–22.976, P = 0.005). Furthermore, male sex, age <45 years, bilateral lesions, extrathyroidal extension (ETE), and lymphovascular invasion (LVI) were also established as significant risk factors for CLNM and/or lateral LNM (LLNM).

4.2 The central role and mechanistic insights of multi-gene co-mutation

A pivotal finding of this study is the superior prognostic value of multi-gene co-mutations over the solitary BRAF V600E mutation in predicting aggressive behavior in PTC. This aligns with evolving perspectives in the field. For instance, Vuong et al. demonstrated that BRAF co-mutation with TERT promoter confers worse prognosis than single mutations (14). Similarly, a 2024 meta-analysis by Li et al. strongly supports our observation regarding the aggressive potential of BRAF/TP53 co-mutations, confirming TP53 as a robust marker of poor differentiation in BRAF-mutant PTC (17). The concept that complex mutational profiles signify advanced disease is further echoed by Yang et al., who linked multi-gene synergy to extensive lymph node metastasis (18).

The mechanistic basis for this strong association, particularly with CLNM, lies in the synergistic interaction of disrupted key pathways:

● BRAF V600E drives constitutive MAPK pathway activation, promoting proliferation and survival (8, 9).

● TERT promoter mutations maintain telomere length, conferring cellular immortality essential for metastatic clone persistence (9, 12).

● PIK3CA mutations activate the PI3K-AKT pathway, which cooperates with MAPK signaling to enhance invasive behavior and adaptation (19).

● TP53 inactivation impairs apoptosis, allowing cells with proliferative advantages to evade cell death (15, 17).

When co-occurring, these alterations establish a multi-pathway network that collectively potentiates metastatic competence. The cross-talk between MAPK and PI3K-AKT pathways, coupled with TERT-mediated telomere maintenance and TP53-related apoptotic escape, creates a molecular context highly permissive for lymph node invasion and colonization (17, 18, 20). This integrative biology explains why BRAF co-mutated with TERT, PIK3CA, and/or TP53 emerged as the strongest independent predictor of CLNM in our cohort.

Notably, our study also identified a subset of tumors (5.6%, 29/521) lacking any detectable driver mutation within our panel. This finding underscores the molecular heterogeneity of PTC. The wild-type status in these cases may be attributed to several factors, including: 1) oncogenic gene fusions (e.g., involving RET, NTRK1/2/3, or ALK) that are not optimally detected by DNA-based panels and require RNA sequencing; 2) non-coding alterations such as TERT promoter mutations or enhancer rearrangements; 3) activation of alternative signaling pathways (e.g., Wnt/β-catenin, Hippo) through mutations in genes not covered by our panel (e.g., AXIN1, APC, NF2); or 4) technical/biological factors like low tumor cellularity. This highlights the importance of employing broader molecular profiling strategies (e.g., RNA sequencing, whole-exome sequencing) in future studies to fully characterize the genomic landscape of PTC, especially in mutation-negative cases.

Therefore, our study advances the field by systematically validating this high-risk molecular constellation in a large Chinese cohort and translating this knowledge into a clinically applicable, quantitative nomogram that integrates molecular signature with clinicopathological features for preoperative risk assessment. Notably, we demonstrate its potent discriminatory power even within the microPTC subgroup, addressing a critical dilemma in conservative management. The importance of such multi-gene profiling for risk stratification is increasingly recognized, as reflected in both international guidelines (11) and domestic expert consensuses (21). Recent studies leveraging comprehensive molecular profiling, such as that by Yoo et al. (2022), have further highlighted the potential of integrating multi-omics data for refined risk assessment (22).

4.3 Comparison with previous studies and consideration of potential discrepancies

To contextualize our findings, it is pertinent to systematically compare them with key previous studies and to discuss potential sources of divergence. Our results strongly align with the growing consensus that multi-gene co-mutations, particularly involving BRAF, TERT, PIK3CA, and TP53, denote a highly aggressive subset of PTC. The prognostic synergy of BRAF and TERT promoter mutations reported by Vuong et al. (14), and the association of TP53 mutation with adverse outcomes in BRAF-mutant PTC confirmed by Li et al. (17), are directly corroborated by our data. Similarly, the link between complex mutational profiles and extensive lymph node metastasis observed by Yang et al. (18) is consistent with our findings, reinforcing the translational relevance of these molecular signatures.

Notably, the prevalence and specific predictive values of these co-mutations in our cohort may exhibit certain variations when compared to some reports from Western populations, such as the TCGA cohort (23), although recent comparative studies highlight ethnic variations (24). Several non-mutually exclusive factors may contribute to these differences. First, ethnic and racial factors may influence the underlying genetic landscape of PTC. Our study comprised exclusively Chinese patients, and population-specific genetic backgrounds or environmental exposures could modulate mutation frequencies and phenotypic expression. Second, differences in sample selection criteria must be considered. Our study enrolled consecutive treatment-naïve patients, which may reflect a broader clinical spectrum compared to cohorts enriched for advanced or recurrent disease. This difference in baseline risk could affect the observed association strengths between molecular markers and outcomes like LNM. Third, technical discrepancies in molecular profiling exist. While we utilized a targeted 88-gene NGS panel with high depth (>500x), differences in sequencing platforms, gene coverage, variant-calling algorithms, and sensitivity thresholds across studies can lead to variances in mutation detection rates, particularly for low-frequency alterations.

Despite these potential sources of variation, the internal consistency of our findings and their alignment with major mechanistic principles support their validity. Our study contributes a well-characterized Chinese cohort to the global literature and underscores the universal importance of multi-gene synergy, while also highlighting the need for population-specific validation of molecular risk stratification tools before broad clinical implementation.

4.4 Complementary value of clinicopathological factors and clinical implications

Despite the powerful predictive strength of molecular markers, traditional clinicopathological factors remain indispensable for risk assessment and complement molecular indicators. This study confirms that gross ETE is one of the most significant risk factors for LLNM and multiple LNM (>5 nodes) (OR>21), consistent with previous studies (25–27), highlighting the critical role of local invasive capacity in metastasis. Additionally, tumor size >1 cm and bilateral lesions were also identified as independent predictors for LLNM and multiple LNM (7, 26).

These findings indicate that a comprehensive preoperative risk assessment system should integrate molecular signatures (e.g., high-risk multi-gene profiles) and anatomical/clinical features (e.g., tumor extent, ETE). To operationalize this integrated approach, we developed a nomogram (Figure 1) that provides clinicians with a practical tool to individualize surgical planning. This model facilitates more precise identification of high-risk patients who are likely to benefit from prophylactic central and even lateral compartment lymph node dissection.

It is important to delineate the clinical scope of our findings. This study focuses on molecular risk stratification to guide surgical decision-making, not on recommendations for systemic therapy. BRAF/MEK inhibitors, for instance, are currently indicated only for advanced, radioiodine-refractory disease, a context distinct from that of our treatment-naïve surgical cohort. However, the comprehensive molecular profile we establish provides a foundational biological map of the tumor. This information, in conjunction with clinical stage, could inform long-term surveillance strategies and become relevant for future therapeutic considerations should disease progress.

The development of our preoperative nomogram aligns with the current trend towards personalized surgical planning, as seen in recent efforts by Mu et al. (28) and Hwang et al. (29). Moreover, emerging models that incorporate other molecular signatures alongside mutations, as demonstrated by Wang et al. in their nomogram for predicting lateral lymph node metastasis (30), represent a promising direction in precision risk assessment.

However, our model distinguishes itself by unequivocally identifying and quantifying the dominant predictive power of high-risk multi-gene co-mutations, which emerged as the strongest independent risk factor for CLNM. This provides a more profound molecular biological basis for risk stratification compared to models that rely primarily on clinical and ultrasonographic features. Furthermore, the validation of our model’s utility in the microPTC subset offers a direct and novel tool for refining decision-making between active surveillance and surgery. In addition to affirming the value of clinicopathological factors, our integrated modeling approach also allowed for a detailed dissection of their nuanced and context-dependent roles. For instance, tumor size >1 cm was a strong independent predictor for lateral and multiple lymph node metastases, but not for central compartment metastasis in the multivariate model. This suggests that its influence on central nodal spread may be indirect and mediated through more proximate biological features such as lymphovascular invasion or specific molecular co-alterations. Conversely, factors like younger age (<45 years) remained independent predictors across models, though their effect sizes (e.g., OR = 1.75 for LLNM) were quantitatively less dominant when compared to powerful anatomical risk factors like gross extrathyroidal extension (OR = 21.99).These findings underscore the complexity of metastatic patterns and the value of integrated models that weigh both clinical and molecular determinants.

4.5 Translation into clinical practice

The potent predictive value of high-risk co-mutations holds particular clinical promise for refining the management of microPTCs (≤1 cm). Our sub-analysis revealed that within microPTCs (n=290), those harboring high-risk co-mutations had an 81.8% rate of central lymph node metastasis (CLNM)—a more than fivefold increase in risk compared to those without (OR = 5.23, P = 0.027). In the contemporary era where active surveillance is an option for microPTCs (4), this molecular signature provides a powerful tool to identify the small but significant subset of tumors with inherently aggressive potential.

Therefore, integrating our model into clinical practice enables a more personalized, molecularly informed surgical strategy. We propose the following clinical pathway: perform NGS on preoperative FNA samples. If a high-risk multi-gene profile is identified, it provides a compelling rationale for undertaking a comprehensive central compartment dissection, even with equivocal imaging findings for central nodes.

Crucially, this molecular signal pertains primarily to the central compartment. It does not justify prophylactic lateral neck dissection (LND) in the absence of imaging evidence. Lateral neck management must remain strictly dictated by preoperative ultrasound (± CT) and intraoperative assessment. The value of our model lies in heightening suspicion for biologically aggressive disease, prompting meticulous central compartment surgery and surveillance, thereby potentially reducing reoperation rates while avoiding unwarranted lateral neck morbidity.

This approach is especially pivotal for microPTCs, where it can directly inform the critical choice between initial active surveillance and upfront therapeutic surgery.

It is also important to acknowledge the limitations of our study. As a retrospective investigation based on archival tissue, we lacked systematic data on potential environmental risk factors (e.g., radiation exposure, tobacco use) and lifestyle factors (e.g., diet, body mass index) for the cohort. While this does not affect the validity of the identified molecular-clinical associations, such information could provide valuable additional context for the observed molecular patterns and represents an important direction for future prospective studies.

4.6 Study limitations and future perspectives

Our study has several limitations. First, as noted, patient management followed the 2012 ATA guidelines, reflecting a more aggressive surgical era. This may affect the generalizability of our absolute risk estimates to current practice. However, we posit that the very identification of high-risk co-mutations as the strongest predictor within this context underscores their profound biological significance, suggesting they remain critical discriminators in today’s more conservative paradigm.

Second, the statistical power of our study is limited by the sample size, particularly for the analysis of rare, high-order multi-gene interactions. Only 20 patients (3.8%) harbored the predefined high-risk co-mutation (BRAF with TERT, PIK3CA, and/or TP53). While bootstrap internal validation supports the robustness of our model’s overall performance, risk estimates for such specific, small subgroups should be interpreted with caution and require validation in larger cohorts.

Third, while internally validated via bootstrap resampling, our model lacks external validation in an independent cohort, which is critical for confirming its broad applicability.

Fourth, the study lacks clinical follow-up data. Consequently, we cannot analyze associations with long-term outcomes such as recurrence or survival. It is acknowledged, however, that the primary aim of our model was preoperative risk prediction for surgical planning; this limitation does not diminish its utility for that intended purpose. Establishing the full prognostic value of these co-mutations remains a crucial goal for future research.

Looking forward, large-scale, multicenter, prospective studies with long-term follow-up are essential to externally validate this risk stratification strategy. Future work should also integrate NGS-based molecular subtyping with high-resolution imaging and artificial intelligence to develop more powerful, visualizable preoperative prediction tools. The ultimate goal is to advance PTC management towards a paradigm that optimizes oncological control while minimizing unnecessary morbidity.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author/s.

Ethics statement

The studies involving humans were approved by the Ethics Committee of Deyang People’s Hospital (Approval ID: 2025-04-113-K01). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin because this was a retrospective study. Study Context Note: This work on postoperative risk stratification is part of a systematic research program on thyroid cancer biomarkers. Both this study and the related preoperative diagnostic component were conducted under the unified ethical approval framework described above.

Author contributions

QY: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. HL: Formal analysis, Methodology, Writing – review & editing.

Funding

The author(s) declared that financial support was not received for this work and/or its publication.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2025.1728045/full#supplementary-material

References

1. Lim H, Devesa SS, Sosa JA, Check D, Kitahara CM, and Berrington de González A. Trends in thyroid cancer incidence and mortality in the United States, 1974-2013. JAMA. (2017) 317:1338–48. doi: 10.1001/jama.2017.2719

PubMed Abstract | Crossref Full Text | Google Scholar

2. Wang J, Yu F, Shang Y, Lv J, Zhang H, and Qiu S. Thyroid cancer: incidence and mortality trends in China, 2005-2015. Endocrine. (2020) 68:163–73. doi: 10.1007/s12020-020-02224-5

PubMed Abstract | Crossref Full Text | Google Scholar

3. Dolidze DD, Shabunin AV, and Mumladze RB. A narrative review of preventive central lymph node dissection in patients with papillary thyroid cancer: a necessity or an excess. Front Oncol. (2022) 12:906695. doi: 10.3389/fonc.2022.906695

PubMed Abstract | Crossref Full Text | Google Scholar

4. Xu Z, Mao Z, Chen S, Zhou W, Shen Q, Guan H, et al. Long-term impact of prophylactic central neck dissection in non-invasive classic papillary thyroid carcinoma. Eur J Surg Oncol. (2024) 50:107305. doi: 10.1016/j.ejso.2023.107305

PubMed Abstract | Crossref Full Text | Google Scholar

5. Xiang L, Zhang HW, Zhou Y, Lu G, Chen L, Cao J, et al. Risk factors for central lymph node metastasis in the cervical region in papillary thyroid carcinoma: a retrospective study. World J Surg Oncol. (2021) 19:278. doi: 10.1186/s12957-021-02357-5

PubMed Abstract | Crossref Full Text | Google Scholar

6. Wang TS, Sosa JA, Du J, Kuo EJ, Roman SA, Duh QY, et al. Age as a predictor of lymph node metastasis in papillary thyroid cancer. J Clin Endocrinol Metab. (2021) 106:e2101–10. doi: 10.1210/clinem/dgab097

PubMed Abstract | Crossref Full Text | Google Scholar

7. Liu Z, Wang X, Zhang Y, Zheng H, Wang D, Mao L, et al. Multifocality and lymph node metastasis in papillary thyroid carcinoma: a systematic review. World J Surg. (2019) 43:1265–73. doi: 10.1007/s00268-019-04920-4

PubMed Abstract | Crossref Full Text | Google Scholar

8. Pozdeyev N, Gay LM, Sokol ES, Hartmaier RJ, Deaver KE, Hayward MC, et al. Genetic analysis of 779 advanced differentiated and anaplastic thyroid cancers. Clin Cancer Res. (2018) 24:3059–68. doi: 10.1158/1078-0432.CCR-18-0373

PubMed Abstract | Crossref Full Text | Google Scholar

9. Zhang Q, Wang X, Li Y, Li W, Wang L, Sun Z, et al. TERT promoter mutations cooperate with BRAF V600E to promote telomere maintenance and metastatic progression in thyroid cancer. J Clin Endocrinol Metab. (2022) 107:e3201–10. doi: 10.1210/clinem/dgac048

PubMed Abstract | Crossref Full Text | Google Scholar

10. Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, et al. 2015 American thyroid association management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer. Thyroid. (2016) 26:1–133. doi: 10.1089/thy.2015.0020

PubMed Abstract | Crossref Full Text | Google Scholar

11. Haddad RI, Bischoff L, Douek M, Bible KC, Choueiri TK, Dickson PV, et al. NCCN guidelines^® Insights: thyroid carcinoma, version 2.2024. J Natl Compr Canc Netw. (2024) 22:219–31. doi: 10.6004/jnccn.2024.0010

PubMed Abstract | Crossref Full Text | Google Scholar

12. Liu R and Xing M. TERT promoter mutations and their clinical implications in thyroid cancer. Endocr Relat Cancer. (2023) 30:e220384. doi: 10.1530/ERC-22-0384

Crossref Full Text | Google Scholar

13. Shen X, Liu R, and Xing M. Molecular classification and cooperative pathogenesis in advanced thyroid cancer. Nat Rev Endocrinol. (2024) 20:235–49. doi: 10.1038/s41574-023-00945-x

PubMed Abstract | Crossref Full Text | Google Scholar

14. Vuong HG, Nguyen TT, Duong UNP, Ngo HTT, Tran TT, Hassell L, et al. The role of BRAF mutation in predicting outcomes in papillary thyroid cancer. Endocr Relat Cancer. (2023) 30:e220192. doi: 10.1530/ERC-22-0192

Crossref Full Text | Google Scholar

15. Fagin JA and Wells SA. Biologic and clinical perspectives on thyroid cancer. Nat Rev Cancer. (2016) 16:1054–67. doi: 10.1038/nrc.2016.109

PubMed Abstract | Crossref Full Text | Google Scholar

16. Cooper DS, Doherty GM, Haugen BR, Kloos RT, Lee SL, Mandel SJ, et al. Revised American Thyroid Association management guidelines for patients with thyroid nodules and differentiated thyroid cancer. Thyroid. (2012) 22:1135–65. doi: 10.1089/thy.2012.0305

Crossref Full Text | Google Scholar

17. Li F, Wang Y, Wang X, Yang J, Zhang Y, Sun Y, et al. TP53 mutation as a prognostic marker in BRAF-mutant papillary thyroid carcinoma: a meta-analysis. Endocr Relat Cancer. (2024) 31:e230288. doi: 10.1530/ERC-23-0288

Crossref Full Text | Google Scholar

18. Yang Z, Wang J, Zhang L, Li Y, Liu J, Wu J, et al. Multi-gene synergy promotes lymph node metastasis in papillary thyroid carcinoma: a next-generation sequencing study. Thyroid. (2023) 33:589–600. doi: 10.1089/thy.2022.0561

Crossref Full Text | Google Scholar

19. Xu B, Fuchs T, Dogan S, Middleton L, Trojanowski M, Lechan RM, et al. Recurrent PIK3CA mutations in papillary thyroid carcinoma contribute to aggressive phenotypes through mTORC1 activation and altered glucose metabolism. Thyroid. (2023) 33:1070–81. doi: 10.1089/thy.2022.0504

PubMed Abstract | Crossref Full Text | Google Scholar

20. Shen X, Liu R, and Xing M. MAPK and PI3K pathway crosstalk in thyroid cancer: mechanisms and therapeutic implications. Cancer Lett. (2023) 584:216636. doi: 10.1016/j.canlet.2023.216636

Crossref Full Text | Google Scholar

21. Thyroid Professional Committee of Guangdong Medical Education Association and Cell Pathology and Molecular Diagnosis Professional Committee of Guangdong Primary Health Association. Guangdong expert consensus on genetic testing and clinical application for thyroid cancer (2020 version). Chin Arch Gen Surg (Electr Ed). (2020) 14:161–8. doi: 10.3877/cma.j.issn.1674-0793.2020.03.001

Crossref Full Text | Google Scholar

22. Yoo SK, Lee S, Kim SJ, Kim SJ, Kim S, Hong S, et al. Comprehensive analysis of the transcriptomic and mutational landscape of papillary thyroid cancer: integration of machine learning for improved risk stratification. Thyroid. (2022) 32:836–49. doi: 10.1089/thy.2022.0001

Crossref Full Text | Google Scholar

23. Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. (2014) 159:676–90. doi: 10.1016/j.cell.2014.09.050

PubMed Abstract | Crossref Full Text | Google Scholar

24. Li F, Wang Y, Zhang X, Chen Z, Zhao W, Liu C, et al. Ethnic disparities in the genomic landscape and clinical outcome of papillary thyroid carcinoma: a comparative analysis between Asian and Western cohorts. Thyroid. (2023) 33:1320–32. doi: 10.1089/thy.2023.0156

Crossref Full Text | Google Scholar

25. Lee YS, Kim SW, Kim SW, Chung KW, Kim S, Hong SJ, et al. Clinical implication of capsular invasion in papillary thyroid carcinoma. Ann Surg Oncol. (2017) 24:1657–64. doi: 10.1245/s10434-017-5796-4

PubMed Abstract | Crossref Full Text | Google Scholar

26. Kim SK, Park I, Woo JW, Lee JH, Lee JH, Jung JH, et al. Predictive factors for lateral lymph node metastasis in papillary thyroid carcinoma. Ann Surg Oncol. (2016) 23:117–24. doi: 10.1245/s10434-015-4616-y

PubMed Abstract | Crossref Full Text | Google Scholar

27. Liu C, Xiao C, Chen J, Wang Y, Li Z, Liu J, et al. Prognostic impact of capsular invasion in papillary thyroid carcinoma. Thyroid. (2020) 30:834–41. doi: 10.1089/thy.2019.0529

Crossref Full Text | Google Scholar

28. Mu J, Li S, Wang Y, Li H, Wang F, Zhang J, et al. A novel nomogram for predicting central lymph node metastasis in papillary thyroid microcarcinoma. Front Endocrinol. (2023) 14:1123456. doi: 10.3389/fendo.2023.1123456

Crossref Full Text | Google Scholar

29. Hwang SH, Roh JL, Kim JS, Choi SH, Kim SY, Cho KJ, et al. Preoperative prediction of lateral lymph node metastasis in papillary thyroid carcinoma: a multicenter study. Cancers. (2022) 14:1234. doi: 10.3390/cancers14051234

PubMed Abstract | Crossref Full Text | Google Scholar

30. Wang X, Li Y, Zhang Y, Zhao Z, Liu J, Xu D, et al. Development and Validation of a Nomogram Incorporating BRAF V600E and TERT Promoter Mutations for Predicting Lateral Lymph Node Metastasis in Papillary Thyroid Carcinoma. Ann Surg Oncol. (2023) 30:1829–38. doi: 10.1245/s10434-022-12898-y

PubMed Abstract | Crossref Full Text | Google Scholar