Predictive analysis of BRAF V600E mutation and central lymph node metastasis in papillary thyroid carcinoma

Introduction

Thyroid carcinoma (THCA) is the fastest-growing malignancy globally, with papillary thyroid carcinoma (PTC) being the most prevalent subtype (1,2). Although PTC generally has a favorable prognosis due to its indolent nature, lymph node metastasis (LNM) in the neck can occur early, most often first affecting the central compartment. Central LNM (CLNM) is frequently identified only postoperatively, termed occult CLNM (3). In patients with PTC, cervical LNMs markedly increase the risk of recurrence and distant metastasis, making accurate preoperative assessment critical for determining the need for prophylactic central neck dissection (pCLND).

International guidelines recommend therapeutic lymph node dissection for patients with clinically evident central compartment metastasis (cN1 stage). However, the role of pCLND in patients with clinical lymph node negative (cN0) remains controversial. Japanese guidelines support routine pCLND, citing benefits in precise postoperative staging, treatment guidance and potential recurrence reduction (4). Conversely, the 2015 American Thyroid Association (5), 2019 European Society for Medical Oncology (6) and 2022 Chinese Anti-Cancer Association (7) guidelines advise pCLND only for cN0 patients with high-risk features, such as T3-T4 tumors, multicentricity, family history, childhood radiation exposure, or lateral cervical LNM. The 2022 National Comprehensive Cancer Network guidelines do not recommend routine pCLND (8). Although pCLND may decrease cervical lymph node recurrence by 34%, postoperative complication rates can reach 17.7% (9), highlighting the need for individualized preoperative assessment and predictive biomarkers for CLNM.

Currently, pCLND decisions rely heavily on surgical experience. Neck ultrasound, the primary imaging modality, has a sensitivity of only 15–40% for detecting CLNM, per the 8th edition of the AJCC manual (10–12). Computed tomography can be used as a supplement, however the indolent progression of PTC often limits detectable imaging changes associated with CLNM (13,14). Fine-needle aspiration cytology (FNAC) remains the most direct, accurate, and cost-effective preoperative diagnostic tool, enabling detection of molecular biomarkers such as BRAF, RAS and TERT mutations. The BRAF V600E mutation, present in 40–80% of cases, is associated with aggressive PTC features. While some studies suggest it may predict CLNM in patients with cN0 (15), current evidence does not support using BRAF V600E alone to guide pCLND. Integrating additional indicators may enhance predictive accuracy and inform surgical planning (16,17).

To elucidate the relationship between BRAF V600E, clinicopathological features and CLNM in PTC, the present study analyzed gene expression profiles and clinical data from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases. It aimed to develop a risk scoring model for CLNM based on differentially expressed genes (DEGs) stratified by BRAF V600E status. Furthermore, patients with cN0 PTC treated at the General Surgery Department of Shunde Hospital of Guangzhou University of Chinese Medicine (GUCM) were included, with preoperative CLNM tissues collected via ultrasound-guided FNAC and postoperative thyroid tissue analyzed for BRAF V600E. The present study sought to support personalized preoperative assessment, optimize surgical decisions, reduce unnecessary lymph node dissection and improve postoperative quality of life.

Materials and methods

Public data collection and processing

mRNA expression profiles and clinical data for THCA patients were obtained from TCGA database. Level 3 HTSeq-FPKM data were normalized to transcripts per million reads. For external validation, GSE60542 and GSE29265 datasets were retrieved from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Tumor mutation analysis

The prevalence of BRAF mutations across cancer types was initially assessed using cBioPortal (http://www.cbioportal.org). Somatic mutation data from TCGA were analyzed with the R package ‘maftools’ to visualize mutation frequencies and types, focusing on THCA. After confirming BRAF V600E as the predominant mutation in THCA, samples harboring this mutation and their associated clinical data were downloaded from cBioPortal for further analysis.

DEGs' analysis

Data of patients with THCA obtained from TCGA were stratified into BRAF V600E mutant and wild-type groups. Differential expression analysis was performed using the R package limma 3.52.2 (18). DEGs were defined by an adjusted P-value <0.05 and |log2-fold-change (FC)|>1.5.

Functional enrichment analysis

To explore the biological significance of DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using the ClusterProfiler R package 4.4.4 (18).

Construction and validation of the diagnostic risk model

A diagnostic risk model for predicting CLNM in THCA was established using DEGs from both BRAF V600E mutant and wild-type cases via Lasso-Cox regression with the glmnet package 3.0 (glmnet.stanford.edu/). The model was formulated as:

where ‘expi’ denotes gene expression and ‘coefi’ represents the corresponding gene risk coefficient. Model performance was evaluated using Receiver Operating Characteristic (ROC) and Precision-Recall (PR) curves. External validation was performed using the GSE60542 and GSE29265 datasets.

It should be noted that genes such as CST6 and LOX were included solely as part of the predictive model, and their downstream mechanistic roles in CLNM were not experimentally investigated in the present study, representing a limitation.

Immune infiltration analysis

The relative enrichment of 24 immune cell types in THCA was assessed using single-sample Gene Set Enrichment Analysis (ssGSEA) via the R package GSVA 3.19 (19). Spearman's correlation analysis was performed to evaluate associations between the risk score and immune cell infiltration. P-values from the 24 parallel tests were adjusted using the Benjamini-Hochberg false discovery rate (FDR), with significance defined as FDR<0.05. Differences in immune infiltration between high- and low-risk groups were analyzed using the Wilcoxon rank-sum test. Additionally, immune, stromal and ESTIMATE scores were calculated using the ESTIMATE algorithm (20).

The correlation between the risk score and expression of interleukins, chemokines and immune checkpoints was further examined. Immunophenoscore (IPS) was used to predict potential responses to immunotherapy (anti-PD-1 and anti-CTLA-4) based on gene expression profiles, leveraging data from The Cancer Immunome Atlas (https://tcia.at/).

Construction and validation of the nomogram

A binary logistic regression model was constructed using the glm function in R to predict overall survival probability. The RMS package was used to develop and visualize the nomogram. Calibration curves and restricted cubic spline plots were employed to evaluate the nomogram's performance. The concordance index (C-index) quantified its discriminative ability and decision curve analysis was performed to assess net clinical benefit.

Retrospective analyses of pathological samples

Clinical data and ultrasound-guided FNAC results were collected from 32 patients with cN0 PTC who underwent preoperative assessment at Shunde Hospital of GUCM between August 2022 and November 2023 (Foshan, China). Inclusion criteria were: i) FNAC-confirmed or highly suspected PTC with TI-RADS 5 nodules; ii) No LNM on preoperative ultrasound/CT, classified as cN0 per Kowalski criteria; iii) Central lymph node diameter >0.5 cm. Exclusion criteria included prior thyroid surgery, history of thyroid disease treatments (including I-131 therapy), and incomplete clinical data. FNAC was performed by an experienced physician, and the specimens were sent to Guangzhou Da'an Clinical Laboratory for BRAF V600E mutation detection using ARMS-PCR (cat. no. 56404 LOT:026ABB01; Wuxi Shenrui Bio-Pharmaceuticals Co. Ltd.).

Postoperative immunohistochemical (IHC) analysis of BRAF V600E mutation

Clinical and pathological data from 222 patients with PTC treated between January 2022 and November 2023 were collected. Inclusion criteria were: i) age ≥18 years; ii) histologically confirmed PTC; and iii) first diagnosis of THCA. Exclusion criteria included incomplete data, non-primary THCA, or presence of other malignancies. Collected variables included age, sex, histopathological subtype, tumor location and size, central lymph node status, IHC results, and surgical approach. CLNM was independently assessed by two blinded pathologists to minimize observer bias. IHC was used to detect BRAF V600E protein expression, with positivity determined by staining intensity and cytoplasmic localization.

Ethics statement

The present study was approved by the Ethics Committee of Shunde Hospital of Guangzhou University of Chinese Medicine (approval nos. KY2022066 and KY2023127; Foshan, China). All procedures were performed in accordance with relevant guidelines and regulations.

Statistical analysis

Group differences were evaluated using the Wilcoxon test. Associations between categorical variables were assessed using the Chi-squared test or Fisher's exact test, as appropriate; Fisher's exact test was applied when more than 20% of the cells in a contingency table had an expected count of less than 5. Spearman correlation analysis was used to examine associations between risk score and immune infiltration. P<0.05 was considered to indicate a statistically significant difference, with multiple testing corrections applied using the Benjamini-Hochberg FDR method.

Results

Characteristics of the TCGA cohort

A total of 508 TCGA patients pathologically diagnosed with PTC with available clinical and RNA-sequencing data were analyzed. Among these, 490 samples had documented BRAF mutation status and were included in further analyses. Clinicopathological characteristics of the cohort are summarized in Table I.

Table I. Clinical pathological characteristics of patients with papillary thyroid carcinoma from TCGA.

Table I.

Clinical pathological characteristics of patients with papillary thyroid carcinoma from TCGA.

		BRAF V600E
Characteristic	Overall, no. (%)	Mut, no. (%)	Non-mut, no. (%)	P-value
n	508	285	205
Sex (%)				0.822
Female	371 (73.2)	208 (73.2)	152 (74.1)
Male	136 (26.8)	76 (26.8)	53 (25.9)
Age, n (%)				0.772
≤45	239 (47.1)	152 (53.5)	107 (52.2)
>45	268 (52.9)	132 (46.5)	98 (47.8)
Pathologic T stage, n (%)				0.001
T1	144 (28.5)	17 (6)	5 (2.5)
T2	167 (33.1)	79 (27.9)	82 (40.2)
T3	171 (33.9)	76 (26.9)	64 (31.4)
T4	23 (4.6)	111 (39.2)	53 (26)
Pathologic N stage, n (%)				<0.001
N0	231 (50.5)	155 (59.2)	63 (35.2)
N1	226 (49.5)	107 (40.8)	116 (64.8)
Pathologic M stage, n (%)				0.735a
M0	283 (96.9)	170 (97.1)	104 (96.3)
M1	9 (3.1%)	5 (2.9)	4 (3.7)
Pathologic stage, n (%)				<0.001
Stage I	285 (56.4)	38 (13.4)	15 (7.4)
Stage II	52 (10.3)	18 (6.4)	32 (15.7)
Stage III	113 (22.4)	152 (53.7)	122 (59.8)
Stage IV	55 (10.9)	75 (26.5)	35 (17.2)
Ethnicity (%)				0.478
Asian and Black or African American	80 (19.3)	45 (18.3)	32 (21.2)
White	334 (80.7)	201 (81.7)	119 (78.8)
Histological type, n (%)				<0.001
Classical	359 (77.9)	232 (93.9)	111 (56.3)
Follicular	102 (22.1)	15 (6.1)	86 (43.7)
Residual tumor (%)				0.113*
R0	389 (87.4)	215 (85)	163 (91.6)
R1	52 (11.7)	35 (13.8)	14 (7.9)
R2	4 (0.9)	3 (1.2)	1 (0.6)
Primary neoplasm focus type, n (%)				0.551
Multifocal	228 (45.9)	149 (53)	111 (55.8)
Unifocal	269 (54.1)	132 (47)	88 (44.2)
Neoplasm location, n (%)				0.158
Bilateral	86 (17.2)	111 (39.5)	97 (48%)
Isthmus	22 (4.4)	102 (36.3)	70 (34.7)
Left lobe	178 (35.5)	52 (18.5)	29 (14.4)
Right lobe	215 (42.9)	16 (5.7)	6 (3)
Thyroid gland disorder history (%)			0.047
Lymphocytic Thyroiditis	72 (16.1)	12 (4.8)	12 (6.6)
Nodular Hyperplasia	69 (15.4)	31 (12.4)	35 (19.2)
Normal	281 (62.7)	36 (14.3)	34 (18.7)
Other, specify	26 (5.8)	172 (68.5)	101 (55.5)

{ label (or @symbol) needed for fn[@id='tfn1-ol-31-1-15397'] } Categorical data are presented as n (%).

a P-values were obtained using Fisher's exact test (Monte Carlo method, based on 10,000 sampled tables) due to the presence of expected cell counts <5.

BRAF genetic variations in the TCGA cohort

The mutational landscape of BRAF was analyzed across 10,967 samples from 32 cancer types using cBioPortal. THCA exhibited the highest BRAF mutation frequency, accounting for 59.6% of cases (Fig. 1A). Examination of THCA mutations, illustrated by a waterfall plot (Fig. 1B), showed that missense mutations were predominant (77.73%). Analysis of 284 mutation sites across amino acids 0–766, including one duplicate, confirmed that BRAF V600E mutations were exclusively missense (Fig. 1C). Chi-squared analysis revealed a significant association between BRAF V600E and CLNM. The mutation group exhibited a CLNM rate of 59.1%, significantly higher than 35.7% in the non-mutation group (Fig. 1D). A Sankey diagram (Fig. 1E) further illustrated the strong link between BRAF V600E mutations and higher CLNM incidence, whereas the non-mutation group displayed minimal CLNM.

Figure 1. Landscape of BRAF mutations and CLNM in PTC. (A) Pan-cancer distribution of BRAF mutations across 32 tumor types in TCGA, with THCA showing the highest frequency. (B) Waterfall plot of BRAF mutation types in THCA, highlighting predominance of missense mutations. (C) Lollipop plot showing amino acid mutation positions in THCA; V600E is the most frequent variant. (D) CLNM incidence comparison between BRAF V600E mutant and wild-type patients (P<0.001). (E) Sankey diagram illustrating the relationship between BRAF V600E mutation status and CLNM occurrence. CLNM, central lymph node metastasis; PTC, papillary thyroid carcinoma; THCA, thyroid carcinoma; TCGA, The Cancer Genome Atlas.

Identification of DEGs in THCA

A total of 229 DEGs were identified between BRAF V600E mutant and wild-type groups, including 57 upregulated (24.9%) and 172 downregulated genes (75.1%), using thresholds of adjusted P<0.05 and |log2-FC|>1.5 (Fig. 2A and B; Table SI).

Figure 2. Differential gene expression and functional enrichment in BRAF V600E PTC. (A) Volcano plot showing significantly up- and downregulated genes between BRAF V600E mutant and wild-type tumors. (B) Heatmap of top DEGs stratified by mutation status. (C and D) KEGG and GO enrichment analyses of upregulated genes highlighting thyroid hormone biosynthesis and stress response pathways. (E and F) KEGG and GO enrichment analyses of downregulated genes implicating immune regulation and infection-related pathways. DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

KEGG pathway analysis indicated that upregulated genes were significantly enriched in ‘thyroid hormone synthesis, thyroid hormone signaling and Rap1 signaling pathways (Fig. 2C). GO analysis showed enrichment in thyroid hormone metabolism, hormone generation, and metal ion stress response processes (Fig. 2D). By contrast, downregulated genes in the non-BRAF V600E group were associated with KEGG pathways including viral myocarditis, type I diabetes mellitus, and toxoplasmosis (Fig. 2E). GO analysis highlighted their involvement in regulating monocyte, lymphocyte and leukocyte proliferation (Fig. 2F).

Development and validation of a diagnostic model using DEGs

Based on the identified association between BRAF V600E mutation and CLNM (Fig. 3A and B), the authors evaluated whether the 229 DEGs could serve as a diagnostic signature. Lasso regression was used to construct a risk score model:

Figure 3. Development and validation of a diagnostic risk model for central lymph node metastasis based on BRAF-associated DEGs. (A and B) LASSO regression and cross-validation plots are used to identify the optimal gene signature. (C and D) ROC and precision-recall curves showing model performance in The Cancer Genome Atlas cohort (AUC=0.710). (E and F) External validation using GSE60542 dataset (AUC=0.666). (G and H) External validation using GSE29265 dataset (AUC=0.905), confirming model generalizability. DEGs, differentially expressed genes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; AUC, area under the curve.

This model stratified THCA samples into high- and low-risk groups and achieved a ROC AUC of 0.710 (95% CI: 0.664–0.755), indicating acceptable accuracy (Fig. 3C and D). External validation in GEO datasets yielded ROC AUCs of 0.666 (95% CI: 0.542–0.791) in GSE60542 and 0.905 (95% CI: 0.683–1.0) in GSE29265, demonstrating moderate to high predictive performance (Fig. 3E-H).

Incorporation of CLNM status showed significantly higher risk scores in the N1 group, reflecting elevated CLNM rates in high-risk samples (Fig. 4A). Analysis of BRAF V600E-related signature genes revealed that CST6, CLDN10, LOX, and SCEL were upregulated in the LNM group, whereas MT1G, SOD3, and KCNAB1 were higher in the N0 group (Fig. 4B). Individual gene ROC AUCs for LNM prediction were: CST6 (0.64), MT1G (0.657), CLDN10 (0.647), SOD3 (0.648), KCNAB1 (0.645), LOX (0.661), SCEL (0.657) and C3 (0.549), indicating moderate discriminative ability (Fig. 4C).

Figure 4. Expression and predictive value of model signature genes. (A) Boxplot comparing risk scores between CLNM-positive (N1) and -negative (N0) patients, showing higher scores in the metastatic group (P<0.001). (B) Heatmap of signature gene expression across CLNM subgroups. (C) ROC curves showing individual gene performance for CLNM prediction. ***P<0.001. CLNM, central lymph node metastasis; ROC, receiver operating characteristic.

Association between risk score and immune landscape

The relationship between risk scores and the immune landscape in THCA were evaluated by comparing immune cell composition between low- and high-risk groups. ssGSEA demonstrated a strong positive correlation (R>0.5) between risk score and multiple immune cell types, including dendritic cells (DCs), macrophages, immature DCs (iDCs), Th1 cells, neutrophils, Treg cells, and Th2 cells, with significance maintained after Benjamini–Hochberg FDR correction (Fig. 5A). High-risk samples exhibited higher proportions of these immune cells (Fig. 5B), as well as elevated ESTIMATE, stromal and immune scores (Fig. 5C-E), indicating increased immune infiltration.

Figure 5. Association between risk score and tumor immune microenvironment in THCA. (A) Correlation heatmap of risk score with immune cell infiltration based on single-sample Gene Set Enrichment Analysis. (B) Comparison of immune cell type abundances between high- and low-risk groups. (C-E) Immune, stromal and ESTIMATE scores, indicating higher immune/stromal content in high-risk patients. (F and G) Differential expression of interleukins, chemokines and receptors, mostly elevated in high-risk patients. *P<0.05, **P<0.01, ***P<0.001. THCA, thyroid carcinoma.

Analysis of the tumor immune microenvironment revealed that most interleukins, except IL12A, IL17D and IL34, were upregulated in the high-risk group (Fig. 5F). Similarly, chemokines and their receptors were generally elevated, except for CCL16, CCL25 and CCR10 (Fig. 5G). Immune checkpoint analysis showed higher expression of inhibitory molecules, including BTLA, CD244, CD274, CSF1R, CTLA4, HAVCR2, IL10, KDR, LGALS9, PDCD1LG2, TGFβ1, TGFβR1, TIGIT and VTCN1, in the high-risk group (Fig. 6A), alongside increased levels of most stimulatory checkpoint molecules (Fig. 6B). These findings suggest that high-risk THCA patients exhibit enhanced immune infiltration with both inhibitory and stimulatory components, highlighting potential responsiveness to immunotherapy.

Figure 6. Immune checkpoint expression and IPS comparison between risk groups. (A and B) Differential expression of immune checkpoint inhibitors and stimulators in high- vs. low-risk patients, with upregulation in the high-risk group. (C-F) IPS analyses under four immunotherapy conditions (PD-1+/CTLA4+, PD-1+/CTLA4−, PD-1−/CTLA4+, PD-1−/CTLA4−), indicating enhanced immunogenicity in high-risk patients. IPS, immunophenoscore; PD-1, programmed cell death protein 1; CTLA-4, cytotoxic T lymphocyte-associated protein 4. *P<0.05, **P<0.01, ***P<0.001.

Finally, IPS were higher in the high-risk group across CTLA4-negative PD-1-positive, CTLA4-positive PD-1-negative, and CTLA4-positive PD-1-positive conditions (Fig. 6C-F), indicating a more active and potentially immunotherapy-responsive tumor microenvironment (TME).

Nomogram development for risk score in independent diagnostic analysis

Logistic regression analysis of TCGA clinical data identified several factors significantly associated with increased CLNM risk: age (HR=2.032, 95% CI: 1.298–3.179), T2 stage (HR=1.902, 95% CI: 1.086–3.331), T3 stage (HR=2.709, 95% CI: 1.555–4.719), T4 stage (HR=12.940, 95% CI: 3.309–50.602), and risk score (HR=3.910, 95% CI: 2.518–6.071; Table II).

Table II. Univariate and multivariate analyses of central lymph node metastasis status in patients with papillary thyroid carcinoma.

Table II.

Univariate and multivariate analyses of central lymph node metastasis status in patients with papillary thyroid carcinoma.

Characteristic	Total (n=424)	OR (95% CI) univariate analysis	P-value	OR (95% CI) multivariate analysis	P-value
Age, years
>45	225	Reference		Reference
≤45	199	1.664 (1.133–2.445)	0.009	2.032 (1.298–3.179)	0.002
Sex
Male	112	Reference		Reference
Female	312	0.612 (0.395–0.947)	0.027	0.716 (0.435–1.178)	0.188
Focus type
Unifocal	227	Reference		Reference
Multifocal	197	1.601 (1.090–2.352)	0.016	1.469 (0.898–2.404)	0.126
Location
Right lobe	176	Reference		Reference
Left lobe	152	1.343 (0.868–2.078)	0.186	1.250 (0.763–2.046)	0.376
Bilateral	77	2.281 (1.317–3.953)	0.003	1.723 (0.868–3.423)	0.120
Isthmus	19	1.895 (0.727–4.943)	0.191	1.504 (0.539–4.197)	0.436
T
T1	124	Reference		Reference
T2	131	1.659 (0.998–2.758)	0.051	1.902 (1.086–3.331)	0.024
T3	149	3.176 (1.930–5.228)	<0.001	2.709 (1.555–4.719)	<0.001
T4	20	11.472 (3.180–41.388)	<0.001	12.940 (3.309–50.602)	<0.001
Group
Low	205	Reference		Reference
High	219	4.613 (3.063–6.948)	<0.001	3.910 (2.518–6.071)	<0.001

Using these variables, a nomogram was constructed to estimate the probability of CLNM in patients with THCA (Fig. 7A). Calibration analysis showed a C-index of 0.769 (95% CI: 0.723–0.814), indicating favorable model fit (Fig. 7B). Further validation using a restricted cubic spline plot demonstrated overall significance and indicated a non-significant predominantly linear relationship between predictors and outcome (Fig. 7C). Decision curve analysis confirmed that the combined model of age, T stage and risk score provided the highest net clinical benefit (Fig. 7D).

Figure 7. Construction and validation of a nomogram for predicting CLNM in papillary thyroid carcinoma. (A) Nomogram integrating age, T stage and risk score for individualized prediction of CLNM probability. (B) Calibration plot demonstrating agreement between predicted and observed probabilities (C-index=0.769). (C) Restricted cubic spline analysis confirming linear relationship between risk score and CLNM risk. (D) Decision curve analysis showing net clinical benefit of combined model compared with individual predictors. CLNM, central lymph node metastasis; CI, confidence interval.

Predictive value of preoperative BRAF V600E mutation for CLNM

To evaluate the predictive value of BRAF V600E for CLNM in PTC, clinical data from 36 patients with cN0 PTC were analyzed who underwent preoperative ultrasound-guided central lymph node FNAC at Shunde Hospital, GUCM, between August 2022 and November 2023 were analyzed. After excluding four patients with insufficient samples, 32 patients were included, all of whom underwent total thyroidectomy (TT) and pCLND (Fig. 8A).

Figure 8. Predictive utility of preoperative BRAF V600E testing via FNAC. (A) Flowchart of patient inclusion and analysis for preoperative FNAC-based BRAF V600E detection. (B) Comparison of CLNM incidence in patients with and without BRAF V600E mutation, showing higher but non-significant risk in mutant group. (C) Sankey diagram illustrating CLNM distribution by BRAF mutation status. FNAC, fine-needle aspiration cytology; CLNM, central lymph node metastasis; PTC, papillary thyroid carcinoma.

BRAF V600E mutation was detected in four patients using ARMS-PCR, three of whom had CLNM (Table III). Patients were categorized into non-mutation (n=28) and mutation (n=4) groups. CLNM rates were 32.1% in the non-mutation group and 75.0% in the mutation group (Fig. 8B). Although the mutation group showed higher CLNM incidence, the difference was not statistically significant. A Sankey diagram (Fig. 8C) illustrated that most BRAF V600E-positive patients experienced CLNM, whereas most non-mutation patients did not.

Table III. Correlation between BRAF V600E mutation and CLNM in post-operative in patients with papillary thyroid carcinoma.

Table III.

Correlation between BRAF V600E mutation and CLNM in post-operative in patients with papillary thyroid carcinoma.

	Lymph node metastasis, n (%)
Group	Yes	No
non-Mut (n=28)	9 (32.1)	19 (67.9)
Mut (n=4)	3 (75.0)	1 (25.0)

Postoperative BRAF V600E mutation and CLNM risk

A retrospective analysis of 222 PTC cases was performed to assess the association between BRAF V600E and CLNM risk. After excluding 42 cases lacking BRAF testing or lymph node dissection, 180 patients were analyzed, all of whom underwent TT with ipsilateral pCLND (Fig. 9A).

Figure 9. Postoperative IHC detection of BRAF V600E and CLNM association in PTC. (A) Flow diagram summarizing postoperative cohort selection. (B) Representative IHC images showed negative/weak/moderate/strong signal intensity of BRAF V600E expression (scale bar, 100 µm). (C) Bar chart comparing CLNM rates between BRAF-mutant and wild-type groups; higher incidence observed in mutant group without statistical significance. (D) Sankey diagram showing CLNM distribution according to postoperative BRAF mutation status. CLNM, central lymph node metastasis; PTC, papillary thyroid carcinoma; IHC, immunohistochemistry.

Patients were grouped by BRAF V600E status into non-mutation (n=100) and mutation (n=80) groups (Table IV). Representative IHC images of BRAF V600E expression in four patients are demonstrated in Fig. 9B, with additional details in Table V. CLNM occurred in 28.0% of non-mutation and 33.8% of mutation patients, with no statistically significant difference (Fig. 9C). A Sankey diagram (Fig. 9D) shows a higher proportion of CLNM in the BRAF V600E mutation group compared with the non-mutation group.

Table IV. Clinical pathological characteristics of patients with papillary thyroid carcinoma (n=180) from Shunde Hospital of GUCM.

Table IV.

Clinical pathological characteristics of patients with papillary thyroid carcinoma (n=180) from Shunde Hospital of GUCM.

		BRAF V600E
Characteristics	Overall no. (%)	Mut, no. (%)	Non-mut, no. (%)	P-value
n		80	100
Sex, n (%)				0.480
Female	149 (82.8)	68 (37.8)	81 (45)
Male	31 (17.2)	12 (6.7)	19 (10.6%)
Age, n (%)				0.070
>45	81 (45)	42 (23.3)	39 (21.7)
≤45	99 (55)	38 (21.1)	61 (33.9)
T stage, n (%)				0.408a
T1	173 (96.1)	76 (42.2)	97 (53.9)
T2	6 (3.3)	4 (2.2)	2 (1.1)
T3	1 (0.6)	0 (0)	1 (0.6)
N stage, n (%)				0.405
N1	55 (30.6)	27 (15)	28 (15.6)
N0	125 (69.4)	53 (29.4)	72 (40)
Primary neoplasm focus type, n (%)				0.506
Unifocal	137 (76.1)	59 (32.8)	78 (43.3)
Multifocal	43 (23.9)	21 (11.7)	22 (12.2)
Neoplasm location, n (%)				0.307a
Right	78 (43.3)	29 (16.1)	49 (27.2)
Left	70 (38.9)	36 (20)	34 (18.9)
Bilateral	27 (15)	12 (6.7)	15 (8.3)
Isthmus	5 (2.8)	3 (1.7)	2 (1.1)

{ label (or @symbol) needed for fn[@id='tfn3-ol-31-1-15397'] } Categorical data are presented as n (%).

a P-values were obtained using Fisher's exact test (Monte Carlo method, based on 10,000 sampled tables) due to the presence of expected cell counts <5.

Table V. Clinical data and description of immunohistochemistry results.

Table V.

Clinical data and description of immunohistochemistry results.

Protein	Tissue	Age, years	Sex	Location	Quantity	Intensity
BRAF V600E	PTC	54	Female	-	-	Negative
BRAF V600E	PTC	35	Male	Cytoplasmic	40%	Weak
BRAF V600E	PTC	27	Male	Cytoplasmic	10%	Moderate
BRAF V600E	PTC	55	Female	Cytoplasmic	95%

[i] Strong PTC, papillary thyroid carcinoma.

Discussion

The rising incidence of THCA is largely attributable to improved detection methods, yet overall mortality remains low, and numerous thyroid nodules exhibit indolent behavior (21). Nonetheless, LNM markedly increases recurrence and mortality risk, making its prediction and management a central clinical concern. In classical PTC, CLNM is particularly relevant to prognosis. Reported prevalence of CLNM ranges from 20–90%, while preoperative ultrasound demonstrates limited sensitivity, occasionally as low as 12.1% (22). In total, ~60% of patients with PTC are confirmed to have CLNM during initial surgery despite negative preoperative imaging (23), highlighting the need for reliable preoperative predictive models to better guide surgical decision-making.

Routine pCLND remains controversial. Critics argue that it offers limited survival or recurrence benefit while increasing postoperative complications. A total of ~14% of patients undergoing TT with pCLND develop temporary hypoparathyroidism, and 4% experience permanent hypoparathyroidism (24). The modest reduction in recurrence (0.66%) is counterbalanced by a 1.83% increase in temporary hypoparathyroidism (25), and risks of parathyroid and recurrent laryngeal nerve injury (26). Meta-analyses report permanent hypoparathyroidism at 1.1%, permanent recurrent laryngeal nerve injury at 0.5%, and recurrence at 2.8% in the pCLND group (27), with vocal cord paralysis ranging 3.28–27.8% (28). Taken together, these data highlight that indiscriminate pCLND may cause more harm than benefit, thereby emphasizing the potential utility of preoperative diagnostic models to refine indications for CLND and minimize morbidity.

The present study identified age, T stage and BRAF V600E mutation status as significant predictors of CLNM, consistent with established risk stratification frameworks. The ATA risk system guides postoperative management, particularly for radioactive iodine therapy and follow-up intensity (5). CLNM is a key determinant in assigning intermediate or high ATA risk categories. Ghaznavi et al (29) demonstrated that integrating AJCC staging, ATA risk, and age refines disease-specific survival estimates, especially in younger patients with high-risk features. Building on this, our nomogram, incorporating preoperative variables (BRAF V600E, age and tumor size), provides an additional tool to stratify risk before surgery. Importantly, the model showed favorable calibration and discrimination in the TCGA cohort, with decision curve analysis confirming superior clinical benefit compared with individual predictors. Thus, preoperative BRAF-based models and postoperative ATA risk stratification may be complementary, together enabling more individualized treatment strategies. Importantly, our prognostic model did not include classical genes such as TERT or KRAS, consistent with recent literature (30), emphasizing the unique predictive value of BRAF V600E. Moreover, the LASSO-derived gene signature demonstrated strong performance, with external validation in GSE29265 achieving an AUC of 0.905.

The relationship between age, tumor size and CLNM remains debated. Some studies report no significant association (31), whereas others, including Ahn et al (32), identify age ≤45 years as an independent risk factor, consistent with the present findings. Tumor size is also recognized as a risk factor, though thresholds vary: ≥1 cm per Ahn et al (32) versus ≥0.25 cm per Yan et al (33). Our dataset included few nodules ≥4 cm, which may limit statistical accuracy. These discrepancies highlight the need for larger datasets to clarify cutoff values for more robust clinical applications.

Molecular insights have highlighted the significance of BRAF mutations in PTC. BRAF is the most frequently mutated gene in THCA, with V600E promoting cell proliferation and tumor aggressiveness (34,35). Preoperative detection of BRAF V600E informs surgical planning, recurrence risk stratification and potential I-131 therapy (36–38). The results of the present study confirmed that BRAF V600E-positive patients show higher recurrence and metastasis rates, supporting more intensive follow-up and, where appropriate, more aggressive management.

The link between BRAF V600E and CLNM, however, remains controversial. Xing et al (39,40) suggested that this mutation significantly increases the likelihood of CLNM in patients with PTC. In the present study, preoperative FNAC with ARMS-PCR enabled detection of BRAF V600E in patients with cN0, suggesting potential utility in guiding prophylactic CLND. Further analysis of BRAF V600E-associated DEGs identified LOX and CST6 as potential mediators: LOX mediates extracellular matrix remodeling and metastasis (41–43), whereas CST6 may modulate the TME and metastatic potential (44). These genes could act as downstream or parallel mediators of BRAF-driven MAPK signaling, rationalizing their association with elevated CLNM risk.

Immune microenvironment analysis reinforced the clinical implications of our risk model. High-risk patients with THCA displayed a distinctive ‘inflamed’ immune phenotype, with substantial infiltration of DCs, macrophages, neutrophils, Th1/Th2 cells and Tregs, along with higher stromal and immune scores. This pattern, consistent with Xie et al (45), indicates an activated yet functionally constrained immune context. The broad upregulation of cytokines and chemokines, including immunosuppressive mediators such as IL-10 and TGF-β (46), together with increased expression of multiple immune checkpoint molecules (for example, PD-L1, CTLA-4, TIGIT and HAVCR2) (47), reflects T-cell exhaustion. Notably, IPS analysis suggested that high-risk patients may respond more favorably to immune checkpoint inhibitors, particularly combined PD-1/CTLA-4 blockade, suggesting that despite worse prognosis, they could derive substantial benefit from immunotherapy (48). Thus, the current findings not only delineate metastatic risk but also identify an immunological subset of patients potentially suited for targeted immunotherapies.

Several methodological and cohort-related limitations warrant consideration. Because of the limited volume of FNAC specimens, ARMS-PCR was used for BRAF V600E detection, whereas IHC was applied to postoperative samples owing to its practicality in routine pathology. As IHC is generally less sensitive than molecular methods such as PCR or NGS, this methodological difference may introduce false negatives, and the relatively small FNAC cohort, particularly the limited number of BRAF-positive cases, further reduced statistical power. In addition, the detection rates of CLNM varied across cohorts (TCGA: BRAF+ 59.1% vs. BRAF- 35.7%, P<0.001; FNAC: 75% vs. 32.1%, P>0.05; IHC: 33.8% vs. 28.0%, P>0.05), which may reflect differences in detection techniques, population heterogeneity (Western vs. Chinese patients), and sample sizes; although the nomogram demonstrated acceptable discrimination in TCGA (AUC=0.710), its performance in external validation in GSE60542 was modest (AUC=0.666), underscoring the need for larger and more homogeneous cohorts as well as more sensitive molecular approaches (49). Finally, the nomogram of the present study was designed as a preoperative supplementary tool but was not directly compared with the ATA Risk Stratification System, an important limitation that should be addressed in future prospective studies.

Ultimately, these limitations highlight the need for future studies that incorporate larger and more homogeneous cohorts, employ sensitive molecular detection methods, and evaluate head-to-head comparisons with established clinical tools. Such efforts will facilitate the integration of molecularly informed predictive models into routine practice, thereby improving preoperative risk stratification, surgical decision-making, and individualized management of THCA patients.

In summary, the present study identified age, T stage and BRAF V600E-associated high metastatic risk as independent predictors of CLNM in patients with PTC. The nomogram developed herein provides a practical visual tool for preoperative estimation of CLNM risk, supporting more precise surgical planning. For classical patients with PTC with TI-RADS ≥4a nodules and no radiologic evidence of cervical LNM, integrated assessment combining ultrasound-guided fine-needle aspiration with BRAF V600E mutation testing is recommended. This approach enables individualized risk stratification, informing surgical decision-making, reducing unnecessary procedures, and potentially improving patient outcomes and quality of life.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by Self-funded Science and Technology Projects of Foshan (grant no. 2220001005113).

Availability of data and materials

The data generated and/or analyzed in the present study may be found in the TCGA database or at the following URL: (https://portal.gdc.cancer.gov).

Authors' contributions

JP and ZW contributed to conceptualization and study design and were involved in writing the original draft. YH performed the bioinformatic analysis and curated the data. RP was responsible for specimen collection and statistical analysis. ZF conducted lymph node punctures. JW performed the postoperative immunohistochemical analysis of thyroid and lymph nodes. ZW and YH confirm the authenticity of all the raw data. All authors reviewed the data, provided critical revisions, read and approved the final version of the manuscript.

Ethics approval and consent to participate

The study protocol was approved by the Ethics Committee of Shunde Hospital of Guangzhou University of Chinese Medicine (approval nos. KY2022066 and KY2023127; Foshan, China). Written informed consent was obtained from all participants for the use of their clinical and pathological data in the present research.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References