Abstract
Background/Aim: Previous studies suggest that thyroid cancers in the isthmus are associated with more aggressive behavior and unfavorable patient outcomes compared to lobar tumors. This study aimed to characterize the molecular features of isthmus thyroid cancer.
Materials and Methods: Transcriptomic and proteomic data were retrieved from the Cancer Genome Atlas. We performed 1:2 propensity score matching based on age, sex, tumor subtype, size, extrathyroidal extension, lymph node metastasis, and stage, matching 15 papillary thyroid cancers with the BRAF V600E mutation located in the isthmus to 30 tumors in the unilateral lobar area.
Results: Despite balanced oncogenic drivers, BRAF-mutant isthmus tumors exhibited significantly lower BRAF-RAS scores and higher ERK output scores. Additionally, these tumors had a slightly elevated DNA methylation-based stemness index. Principal component analysis revealed distinct gene expression patterns between the two groups. Differential expression and gene set enrichment analysis indicated enrichment in extracellular matrix remodeling, cell adhesion pathways, and TGF-β signaling. Notable upregulated genes included oligoadenylate synthases, TNR, HAS3, and LIF. Hierarchical proteomic clustering highlighted increased co-expression of several oncoproteins and DNA damage response markers in isthmus tumors.
Conclusion: BRAF-mutant papillary thyroid cancers in the isthmus display distinct molecular characteristics, including increased ERK signaling activity. These findings suggest that topographical location independently influences tumor biology and may account for the more aggressive clinical behavior of isthmus tumors at the molecular level.
Introduction
Thyroid cancer is the most common endocrine malignancy, with its incidence continuously rising (1). Most thyroid cancers are well-differentiated, including papillary, follicular, and oncocytic types. Papillary thyroid cancer accounts for about 90% of newly diagnosed cases. Recent advances in molecular pathogenesis reveal that these tumors arise from acquired driver mutations that constitutively activate the mitogen-activated protein kinase (MAPK) signaling pathway (2). Specifically, thyroid tumorigenesis primarily results from clonal, mutually exclusive mutations in either the BRAF or RAS oncogenes, while a subset of tumors arises from gene rearrangements involving receptor tyrosine kinases (2). During disease progression, additional genetic alterations, such as disruption of p53 function and mutations in the TERT promoter, lead to more advanced and less differentiated thyroid cancer (3).
The thyroid is a butterfly-shaped gland composed of two lobes connected by a narrow strip of tissue called the isthmus. Due to the small volume of tissue in the isthmus, thyroid nodules in this area are uncommon. In a retrospective review of patients with thyroid nodules from six referral centers, only 6% were found in the isthmus (4). A study of 1,973 patients with papillary thyroid cancer revealed that 181 (9%) had tumors located in the isthmus (5). Despite their rarity, thyroid cancers in the isthmus were associated with a higher risk of persistent disease after treatment and shorter disease-free survival (6). This raises an intriguing hypothesis that tumor biology may be influenced by the anatomical topography of the thyroid gland.
The Cancer Genome Atlas (TCGA) is a landmark project that comprehensively maps the genomic profiles of over 20,000 primary cancer samples across more than 30 different cancer types. Recently, researchers from Washington University analyzed TCGA data and found that patients with isthmus thyroid cancers were generally younger and had a higher prevalence of BRAF mutations, along with higher ERK scores, as expected (7). However, this analysis did not account for potential baseline clinicopathological differences. To explore whether isthmus tumors have distinct molecular characteristics compared to lobar tumors, we conducted an extended analysis following propensity score matching. The aim of this study was to investigate whether unique molecular signatures exist in the isthmus.
Materials and Methods
Data sources. Demographic, clinical, pathological, transcriptomic, and proteomic data from the TCGA papillary thyroid cancer (THCA) databases were downloaded from the Genomic Data Commons data portal of the National Cancer Institute (https://portal.gdc.cancer.gov/) or obtained from the initial global analysis publication (8). Primary tumor location and transcriptome data from Illumina HiSeq 2000 RNA sequencing were available for 501 patients. Of these, 372 had proteomics data obtained from the M.D. Anderson Reverse Phase Protein Array Core. Two patients diagnosed with poorly differentiated oncocytic carcinoma and follicular carcinoma were excluded from further analysis.
Five genetic ancestral groups (Europeans, East Asians, Africans, Native/Latin Americans, and South Asians), derived from single nucleotide polymorphism array genotyping and exome sequencing data, were identified in our recent analysis (9). We categorized the oncogenic drivers of the study cohort into BRAF V600E mutation, RAS mutations, fusions, and other mutations or copy number alterations. The mutation status of the TERT promoter was retrieved from our previous study (10).
The thyroid differentiation score (TDS) from the initial TCGA analysis is based on the expression levels of 16 genes associated with thyroid metabolism and function, where higher values indicate a more differentiated tumor status (8). The BRAF-RAS score (BRS) is a 71-gene expression metric that classifies thyroid tumors as either BRAF-like or RAS-like (8). The BRS ranges from −1 to +1, with negative scores indicating a BRAF V600E-like signature and positive scores indicating a RAS mutation-like expression. The ERK output score is derived from prior research that identified 52 genes with rapidly changing expression following MEK inhibition in tumors with the BRAF V600E mutation, reflecting the transcriptional output of the ERK pathway (11). A higher ERK output score indicates increased ERK signaling activity. Additionally, researchers from TCGA evaluated cancer stemness across 33 tumor types and developed stemness indices based on mRNA expression (mRNAsi) and DNA methylation (mDNAsi) (12). Higher values of these indices are associated with biological processes active in cancer stem cells and tumor dedifferentiation, as indicated by histopathological grade.
According to the initial global analysis publication, tumor purity was determined using the ABSOLUTE algorithm, which estimates the fraction of cancer cells based on somatic DNA alterations (13). Since transcriptional profiles are derived from bulk tumor samples, stromal and immune scores have been developed to quantify the presence of stromal and immune cells within the tumor microenvironment (14). The Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) score is calculated by summing the stromal and immune scores, which reflect the total levels of stromal cells and immune cell infiltration in tumor tissues. A high ESTIMATE score implies lower tumor purity. Additionally, we employed several deconvolution algorithms to estimate the fractions of specific cell types. xCell uses curated gene signatures for 64 immune and stromal cell types and applies a spillover compensation approach to provide relative enrichment scores (15). Estimating the Proportions of Immune and Cancer cells (EPIC) infers the proportions of major immune populations using predefined expression profiles (16). Furthermore, Microenvironment Cell Populations-counter (MCP-counter) utilizes specific marker genes to quantify the relative abundance of key immune and stromal cell populations (17).
Propensity score matching. To minimize bias from confounding variables related to tumor location, we performed propensity score matching using a 1:2 nearest-neighbor matching algorithm without replacement. Propensity scores were estimated using a logistic regression model that included age, sex, tumor subtype, tumor size, extrathyroidal extension, lymph node metastasis, and tumor-node-metastasis (TNM) staging. Matching was carried out using the logit of the propensity score as the distance metric, with a caliper width of 0.2 standard deviations of the logit. Since TERT promoter mutation status was missing for a subset of patients, we created a binary indicator for missing TERT promoter mutations and included it in the propensity score model to ensure balanced covariate availability between matched groups. Covariate balance before and after matching was evaluated using standardized mean differences, with values less than 0.1 considered indicative of adequate balance. The optmatch package (version 0.10.8) in R was utilized for the propensity score matching process.
Principal component analysis. After propensity score matching, dimension reduction was performed using principal component analysis (PCA) to highlight strong patterns of gene expression. The transcriptome data were transformed into an orthogonal basis, with the first few principal components capturing the greatest variance. The prcomp() function in R version 4.5.0 was used with the top 500 most variable genes, and a scatter plot was created for visualization.
Identification of differentially expressed genes. Differentially expressed genes were identified using the Bioconductor package DESeq2 (version 1.50.2) in R, incorporating tumor purity as a covariate. The Benjamini-Hochberg correction was applied to adjust p-Values for multiple testing, setting the false discovery rate (q-Value) at 5%. Gene Ontology categories of the differentially expressed genes were analyzed using Enrichr, a free web tool for functional enrichment analysis (18).
Gene set enrichment analysis. To address the limitations of overrepresentation analysis, such as neglecting network structures and overlooking genes that exhibit trends without significant changes, gene set enrichment analysis (GSEA) was employed to eliminate arbitrary cutoffs and improve signal detection (19). Classic GSEA was conducted using GSEA v4.4.0 with MSigDB v2025.1 databases, including Hallmark, the Pathway Interaction Database (PID), and Reactome collections. Gene sets containing fewer than 15 or more than 500 genes were excluded. When multiple probes mapped to the same gene symbol, expression values were collapsed to a single gene using the maximum probe value. Genes were ranked based on the signal-to-noise ratio derived from phenotype labels. Statistical significance was determined using 1,000 phenotype-based permutations. Enrichment significance was assessed through normalized enrichment scores and false discovery rate q-Values. The top pathways with the highest normalized enrichment scores were identified to provide insights into biological processes.
Hierarchical cluster analysis. For the propensity-matched cases with available proteomics data, unsupervised hierarchical clustering was performed on Z-score transformed data using average linkage and Pearson correlation (20). A clustered heatmap of the most variable proteins was generated for visual interpretation.
Statistical analysis. To compare categorical variables between groups, we used Fisher’s exact test or the chi-square test. For continuous variables, medians with interquartile ranges (IQRs) are reported, and differences were assessed using the non-parametric Kruskal-Wallis test or the Mann-Whitney U test (21). A p-value of less than 0.05 was considered statistically significant.
Results
Study cohort characteristics. Of the 499 patients with data on primary tumor location and transcriptome, 22 had tumors located in the isthmus, 214 in the right thyroid lobe, 177 in the left thyroid lobe, and 86 had tumors in both lobes (Table I). Patients with isthmus tumors were more likely to have the tall-cell subtype (p=0.009), to be younger at diagnosis (p=0.036), and to be male (p=0.031). They also had a higher likelihood of having African genetic ancestry and a lower likelihood of having East Asian genetic ancestry (p=0.045). Furthermore, isthmus tumors were less likely to exhibit extrathyroidal extension (p<0.001) but more likely to present with lymph node metastasis (p=0.014). There was no significant difference in tumor size or TNM stages. The frequency of BRAF V600E and TERT promoter mutations was higher in isthmus tumors, although this difference did not reach statistical significance. For isthmus tumors, the prevalence of the BRAF V600E mutation was 68% (p=0.318), and the prevalence of the TERT promoter mutation was 21% (p=0.442).
Clinicopathological features of the study cohort with papillary thyroid cancer.
To minimize bias from confounding factors at clinical presentation, propensity score matching was applied to isthmus tumors and those in the unilateral lobar group. Since the BRAF V600E mutation is the most prevalent genetic alteration in papillary thyroid cancer and oncogenic drivers predominantly influence gene expression (22), our analysis focused on 15 isthmus tumors with the BRAF V600E mutation, excluding other less common oncogenic drivers. A 1:2 matching identified 30 lobar tumors with comparable clinicopathological features (Table II). We verified that all matched covariates achieved an acceptable balance and conducted sensitivity analyses using alternative caliper widths and reduced covariate sets, which yielded consistent results. These 45 cases were subjected to further analyses.
Clinicopathological features of BRAF-mutant papillary thyroid cancer after propensity score matching.
Differentiation and ERK signaling. As shown in Figure 1, there was no difference in TDS between isthmus and lobar tumors (median: −0.76 versus −0.72, p=0.501). Notably, isthmus tumors had significantly lower BRS (median: −0.90; IQR=−0.76 to −0.96) compared to lobar tumors (median: −0.74; IQR=−0.63 to −0.89; p=0.013). ERK output scores were significantly higher in isthmus tumors (median: 18.0; IQR=15.5 to 28.6) than in lobar tumors (median: 12.3; IQR=1.7 to 18.8; p=0.023). These results indicate that BRAF-mutant isthmus tumors exhibit a more pronounced BRAF-like signature and higher ERK signaling activity.
Comparisons of the thyroid differentiation score (TDS), BRAF-RAS score (BRS), ERK output score, and DNA methylation-based stemness index (mDNAsi) between BRAF-mutant papillary thyroid cancers located in the isthmus and lobar regions. Horizontal bars indicate median values. *p<0.05; ns, Not significant (Mann-Whitney U-test).
Isthmus and lobar tumors had similar tumor purity (p=0.522), stromal (p=0.163) and immune (p=0.427) scores, as well as ESTIMATE scores (p=0.324). We examined cell type distributions using three deconvolution algorithms and found no significant differences in the proportions of specific immune and stromal cells between tumors located in the isthmus and those in the lobar regions (Figure 2). These findings suggest that tumor location is not significantly associated with stromal cell composition or immune cell infiltration. The mRNAsi values were similar between the two groups (p=0.665); however, isthmus tumors had slightly higher mDNAsi (median: 0.12; IQR=0.11 to 0.13) compared to lobar tumors (median: 0.11; IQR=0.10 to 0.12; p=0.036).
Estimated proportions of immune and stromal cells in BRAF-mutant papillary thyroid cancers in the isthmus compared to those in the lobar regions. Representative cell types are derived from deconvolution algorithms of xCell (A), EPIC (B), and MCP-counter (C). Horizontal bars indicate median values. ns, Not significant (Mann-Whitney U-test).
Overrepresentation analysis. Principal component analysis revealed that the two groups exhibited distinct yet partially overlapping gene expression patterns (Figure 3A). DESeq2 identified 124 differentially expressed transcripts with a q-value below 0.05 (Figure 3B). Gene Ontology analysis of the differentially expressed genes using Enrichr indicated that the most enriched cellular component was the collagen-containing extracellular matrix (Figure 3C). The top molecular function identified was adenylyltransferase activity, which included OAS1, OAS3, and OASL, all of which were upregulated in isthmus tumors. The most enriched biological process was interleukin-27-mediated signaling, which also included the oligoadenylate synthases. Representative protein-coding genes with the highest fold changes are listed in Table III. Among the upregulated genes, several are recognized for their roles in extracellular matrix remodeling, such as tenascin R (TNR), ovochymase 1 (OVCH1), and hyaluronan synthase 3 (HAS3). Additionally, leukemia inhibitory factor (LIF) is involved in strong pro-tumor cytokine signaling, which helps maintain cancer stemness in solid tumors (23). The upregulation of LIF was accompanied by the downregulation of several differentiation or organization markers, including PRLHR, AZGP1, HOXC8, and MEGF11. This is consistent with the increased mDNAsi, suggesting a further loss of differentiation in isthmus tumors.
Principal component analysis and differential gene expression analysis for BRAF-mutant papillary thyroid cancers in the isthmus versus lobar regions. (A) Principal component analysis of gene expression for the two groups. Each dot or triangle represents a sample, while the circles indicate 95% confidence ellipses. (B) Volcano plot illustrating the relationship between gene expression fold change and significance. Red dots indicate differential expression with an absolute fold change greater than 2 and a q-Value of less than 0.05. (C) Top Gene Ontology categories for differentially expressed genes.
Top ten differentially expressed genes ranked by fold change in BRAF-mutant papillary thyroid cancer located in the isthmus.
Gene set enrichment analysis. In Figure 4, we present the most positively enriched pathways identified through gene set enrichment analysis. Interestingly, the Hallmark, PID, and Reactome databases all highlight dysregulated processes related to cell-cell adhesion and polarity in the isthmus. TGF-β signaling downregulates E-cadherin, and the reduction of E-cadherin, along with tight junction components, results in a loss of epithelial integrity, facilitating cancer cell migration and invasion. These findings are consistent with the differentially expressed genes, which involve the extracellular matrix and cell junction assembly. Additionally, a significant positive enrichment of the p53 pathway was observed in the Hallmark collection. This correlates with the upregulation of APOBEC3A, a potent cytidine deaminase associated with kataegis (localized clusters of hypermutation) and break-associated “mutation showers,” which may contribute to genomic instability. The upregulation of certain stress-responsive genes, such as GPR37 and HSPA6, further reflects this trend.
Gene set enrichment analysis of BRAF-mutant papillary thyroid cancers in the isthmus compared to those in the lobar regions. The top enriched pathways from the Hallmark, Pathway Interaction Database (PID), and Reactome databases, along with the corresponding enrichment plots, are presented.
Proteomics clustering. Among the propensity-matched cases, data from 6 isthmus tumors and 22 lobar tumors were available for proteomics analysis. As illustrated in Figure 5, the protein expression profiles of isthmus tumors were closely clustered in the unsupervised hierarchical analysis, corroborating the results of the principal component analysis of mRNA profiles. Of note, several oncoproteins exhibited increased co-expression in the isthmus tumors, including ABL1 (CABL), MAPK1 (ERK2), ETS1, and LCK. Additionally, increased protein abundance was observed in molecules associated with DNA damage response and apoptosis, including BID, BAX, and phosphorylated CHK1 (Ser345). Downregulation of NOTCH1 and upregulation of SERPINE1 (PAI1) were also prominent. These alterations may be linked to enhanced ERK signaling activity in isthmus tumors.
Clustered heatmap of proteomics data for BRAF-mutant papillary thyroid cancers in the isthmus compared to those in the lobar regions. Unsupervised hierarchical clustering was performed using average linkage and Pearson correlation.
Discussion
A few discrepancies were noted between our findings and the previous report (7). Before propensity score matching, we observed significant differences in sex, age, histological subtype, and lymph node metastasis, while the previous report found no statistical significance in these variables. The primary reason for the discrepancies may be that the previous report grouped unilateral lobar tumors and bilateral tumors together as non-isthmus tumors. This could introduce biases, as basic demographics such as age and sex clearly differ between unilobar and bilateral tumors. Although both the Washington University group and our study did not detect a statistical difference in the frequency of BRAF V600E and TERT promoter mutations, the nearly double prevalence of TERT promoter mutations in isthmus tumors necessitates matched-pair analysis to avoid biased comparisons.
After propensity score matching, both groups exhibited identical distributions of oncogenic drivers and TERT promoter mutations. An interesting finding in this study is that lower BRS and higher ERK output scores were still evident in isthmus tumors. This is supported by the distinct grouping of gene expression observed in the principal component analysis and hierarchical clustering of proteomics data. Although we did not observe differences in TDS after propensity matching, it is noteworthy that isthmus tumors had slightly higher mDNAsi. ERK signaling is well known to promote cancer stemness across various cancers, including thyroid cancer (24). Recent evidence suggests that STAT3 overexpression increases stemness-related gene expression in thyroid cancer cells (25). In this context, it is noteworthy that PTPRT, which dephosphorylates STAT3, was downregulated in isthmus tumors (Table III). The reduced expression of PTPRT may lead to unchecked STAT3 signaling, coupled with the upregulation of LIF, which activates the JAK/STAT3 pathway for stem cell maintenance (23). These aberrant signaling pathways may exacerbate malignant tumor phenotypes in isthmus tumors and provide plausible mechanisms that explain the molecular differences in tumors located in different regions, despite balanced oncogenic drivers.
Our analysis revealed that many differentially expressed genes are associated with the extracellular matrix, while the top enriched pathways involve cell adhesion and extracellular matrix interactions. This is in agreement with clinical observations that loco-regional and distant metastases are more common in isthmus tumors, even among very low-risk or low-risk patients (26). A meta-analysis of 11 retrospective studies indicates that isthmus tumors are linked to a higher risk of central lymph node metastasis (27). However, our deconvolution analysis did not detect differences in the fractions of endothelial cells, specifically microvascular and lymphatic endothelial cells (p=0.109 and 0.166, respectively; data not shown). Overall lymphatic vessel density has been shown to predict nodal metastasis (28). Consistently, a recent study of the Afirma database indicated that isthmus tumors are more likely to harbor BRAF V600E mutations and ALK/NTRK/RET fusions, along with higher ERK activity and follicular mesenchymal transition scores (29). However, the downregulation of immune response regulation observed in the Afirma study was not replicated in the present study. For instance, while IGHG4 was downregulated in the isthmus tumors (Table III), the distribution of B cells was similar between the two groups. A plausible explanation is that the relative immune coldness may result from the higher frequency of BRAF V600E mutations rather than the isthmus location itself. Another possibility is that heightened cancer stemness promotes immune suppression.
Differences in the embryological origins of the isthmus and thyroid lobes have been proposed to explain the aggressive behavior of isthmus tumors (7). Alternatively, we suggest a biomechanical explanation. The thyroid gland is attached to the trachea by the suspensory ligament, also known as Berry’s ligament. During swallowing, the up-and-down movement of the trachea pulls on Berry’s ligament and the nearby isthmus, placing the isthmus under greater stretch tension than the lobar tissue located farther from the trachea. A previous study found that simulated microgravity-related shear stress increased the number of vimentin-positive thyroid cancer cells and elevated levels of extracellular matrix proteins and TGF-β1 (30). TGF-β1 can upregulate the expression of hyaluronan synthases and promote hyaluronan synthesis (31), which alters the stiffness of the extracellular matrix and may enhance stemness and metastasis. Further research is needed to determine whether the molecular characteristics of isthmus tumors are influenced by biomechanical factors such as mechanosensation and mechanotransduction.
A limitation of this study is the absence of an external validation cohort and the relatively small sample size after matching, necessitating cautious interpretation of these findings. Additionally, some information, such as the secretome, is not included in this multifaceted analysis. Current large-scale workflows may not effectively capture extracellular components. For instance, the secretion of sialyl-fibronectin by thyroid tumors correlates with reduced epithelial-mesenchymal transition and improved prognosis (32). Despite these limitations, our convergence of transcriptome and proteomic analyses may provide deeper insights into the molecular characteristics of isthmus tumors, rather than simply relying on clinicopathological comparisons. In thyroid cancer, disturbances in antioxidative stress defense, mismatch repair, non-homologous end-joining, and DNA damage response pathways contribute to both the initiation and progression of the malignant process (33). We recently proposed that the upregulation of TIMELESS may serve as a protective mechanism against accumulating DNA repair deficiencies during dedifferentiation in thyroid cancer (34). Notably, this study found that oligoadenylate synthases and stress-responsive genes such as GPR37 and HSPA6 were upregulated in isthmus tumors, consistent with the increased DNA damage response and apoptosis observed in the proteomic analysis. Following the DNA damage response, ERK can be activated to mediate anti-apoptotic effects and promote cell survival (35). Accordingly, the increased ERK activity in isthmus tumors may reflect a complex adaptive response to genomic instability and cellular stress, potentially offering novel therapeutic options for radioiodine-refractory tumors.
Conclusion
In this study, we employed propensity score matching to control for baseline clinicopathological differences and utilized a multi-omics approach to characterize BRAF-mutant papillary thyroid cancer located in the isthmus. Our findings indicate that isthmus tumors represent a molecularly distinct subgroup, characterized by heightened ERK pathway signaling and more pronounced BRAF-like molecular characteristics. While we propose a biomechanical hypothesis involving chronic stretch tension from tracheal movement as a potential explanation for these molecular differences, the precise mechanisms underlying the topographical influence on tumor biology warrant further investigation.
Acknowledgements
This work was supported by research grants from the National Science and Technology Council of Taiwan (NSTC-112-2314-B-195-009-MY3) and MacKay Memorial Hospital (MMH-11507 and MMH-E-115-08).
Footnotes
Conflicts of Interest
The Authors declare no conflicts of interest.
Authors’ Contributions
Chi-Yu Kuo: Conceptualization, Methodology, Formal analysis, Investigation, Data Curation, Writing - Original Draft. Sheena Yi-Hsin Cheng: Methodology, Data Curation, Writing - Review & Editing. Yi-Chiung Hsu: Methodology, Formal analysis, Writing – Review & Editing. Shih-Ping Cheng: Conceptualization, Investigation, Writing – Original Draft.
Artificial Intelligence (AI) Disclosure
No artificial intelligence (AI) tools, including large language models or machine learning software, were used in the preparation, analysis, or presentation of this manuscript.
- Received December 20, 2025.
- Revision received February 8, 2026.
- Accepted February 9, 2026.
- Copyright © 2026 The Author(s). Published by the International Institute of Anticancer Research.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) 4.0 international license (https://creativecommons.org/licenses/by-nc-nd/4.0).
References
In this issue
Jump to section
Related Articles
Cited By...
- No citing articles found.