From discovery to implementation

Advances in Clinical and Experimental Medicine

Title abbreviation: Adv Clin Exp Med
Journal Impact Factor (JIF 2025) – 2.5
Journal Citation Indicator (JCI 2025) – 0.40
Scopus CiteScore (2025) – 4.2
Index Copernicus Value (ICV 2024) – 161.00
MNiSW – 70 pts
ISSN 1899–5276 (print), ISSN 2451-2680 (online)
Periodicity – monthly

Download original text (EN)

Advances in Clinical and Experimental Medicine

2026, vol. 35, nr 9, September, p. 1619–1634

doi: 10.17219/acem/214710

Publication type: original article

Thematic category: Immunology; endocrinology and metabolism; molecular biology

Language: English

License: Creative Commons Attribution 3.0 Unported (CC BY 3.0)

Download citation:

  • BIBTEX (JabRef, Mendeley)
  • RIS (Papers, Reference Manager, RefWorks, Zotero)

Cite as:


Chen L, Fu C, Chen Y, Zhao X. Neutrophil-related gene signature predicts prognosis and immune microenvironment patterns in thyroid cancer. Adv Clin Exp Med. 2026;35(9):1619–1634. doi:10.17219/acem/214710

Neutrophil-related gene signature predicts prognosis and immune microenvironment patterns in thyroid cancer

Liang Chen1,A,C,E, Chao Fu1,A,D, Yang Chen1,A,D, Xianbao Zhao1,A,F

1 Tumor Department, Yiwu Central Hospital, China

Graphical abstract


Graphical abstracts

Highlights


• Neutrophil-related genes (NRGs) predict prognosis and overall survival in thyroid carcinoma.
• An NRG-based risk model identifies thyroid cancer patients with significantly poorer overall survival.
• Two molecular subtypes of thyroid carcinoma show distinct immune infiltration and tumor microenvironment profiles.
• NRG biomarkers may support personalized risk stratification and precision treatment strategies in thyroid cancer.

Abstract

Background. Thyroid carcinoma (THCA) represents the most prevalent malignancy of the endocrine system. Although the majority of patients exhibit favorable clinical outcomes, a subset still faces significant clinical challenges, including distant metastasis, recurrence, and poor response to therapeutic interventions, all of which substantially compromise quality of life and overall survival (OS).

Objectives. Neutrophils, as pivotal components of the innate immune system, exhibit dichotomous roles in cancer biology, exerting both tumor-promoting and tumor-suppressive functions. However, the prognostic value of neutrophil-related genes (NRGs) and their contribution to the molecular stratification of THCA remain poorly defined. Our aim was to identify and evaluate NRG biomarkers with the potential to inform personalized therapeutic strategies for patients with THCA.

Materials and methods. We first identified differentially expressed genes (DEGs) between THCA tumor and normal samples. Prognosis-associated DEGs were filtered using univariate Cox regression analysis. Subsequently, least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox analysis were employed to identify a core set of prognostic NRGs. A risk model was constructed based on these core genes and validated in independent testing cohorts. Furthermore, consensus clustering was applied to delineate THCA molecular subtypes based on NRG expression profiles, followed by comprehensive analyses of functional enrichment and immune microenvironment characteristics.

Results. The prognostic signature derived from NRGs demonstrated robust and consistent predictive performance across the training and validation datasets. Patients with high NRG-based risk scores exhibited significantly poorer OS and a more immunologically complex tumor microenvironment (TME) than their low-risk counterparts. Consensus clustering further revealed 2 distinct molecular subtypes of THCA, each characterized by divergent immune infiltration patterns and functional pathway enrichment, suggesting biological heterogeneity with potential therapeutic implications.

Conclusions. Our study highlights NRGs as promising prognostic biomarkers in THCA, uncovering distinct molecular subtypes with potential relevance for personalized therapeutic stratification. These findings underscore the clinical value of NRGs and emphasize the need for further mechanistic investigations to advance our understanding of THCA biology and refine precision treatment strategies.

Key words: thyroid neoplasms, neutrophils, prognostic biomarkers, tumor microenvironment, gene expression profiling

Background

Thyroid carcinoma (THCA) is the most prevalent endocrine tumor.1, 2 It is generally categorized into 3 pathological subtypes: differentiated thyroid cancer (DTC), medullary thyroid cancer (MTC), and anaplastic thyroid cancer (ATC).3 Its incidence is influenced by a variety of factors, most notably geographic region4 and sex, with women exhibiting nearly threefold higher incidence rates compared to men.5 Most cases of THCA exhibit an indolent clinical course and are amenable to conventional treatment strategies, including surgical resection and radioactive iodine therapy.6 While surgery is the standard therapy for MTC and is curative only in intrathyroidal disease, ATC and poorly differentiated thyroid cancer (PDTC) exhibit the most unfavorable prognosis.7 A subset of patients presents with aggressive clinical features, such as extrathyroidal extension,8 vascular invasion, or distant metastases, and consequently faces a markedly poorer prognosis.9 These clinical challenges underscore the pressing need to elucidate the molecular mechanisms driving THCA progression and to identify robust biomarkers that can not only predict therapeutic responses but also stratify patients according to their risk, thereby guiding the development of more precise and effective therapeutic strategies.

Within the tumor microenvironment (TME), a heterogeneous and dynamic network composed of malignant cells, stromal elements, and infiltrating immune populations collectively governs tumor initiation, progression, and therapeutic response.10, 11 Among these, neutrophils represent the most abundant phagocytic leukocyte population, exhibiting both pro-tumorigenic and antitumor functions.12 Tumor-associated neutrophils (TANs) have been shown to possess a spectrum of tumor-suppressive capacities. Neutrophil infiltration has been associated with enhanced CD8+ T-cell recruitment and improved prognosis in colorectal cancer, suggesting their role in promoting antitumor immunity. Conversely, TANs can also contribute to malignant progression by secreting pro-inflammatory chemokines such as interleukin-8 (IL-8) and potent angiogenic factors such as vascular endothelial growth factor (VEGF), thereby promoting neovascularization and tumor expansion.13 Neutrophils also release proteolytic enzymes that degrade components of the extracellular matrix (ECM) and facilitate cancer cell invasion and metastasis.14 This functional duality underscores the complex role of neutrophils in cancer biology; however, their prognostic relevance and immunological context in THCA remain incompletely defined.

Objectives

The primary objective of this study was to identify and evaluate neutrophil-related gene (NRG) biomarkers with the potential to inform personalized therapeutic strategies for patients with THCA. Leveraging transcriptomic data from The Cancer Genome Atlas (TCGA), we sought to elucidate the prognostic value of NRG signatures in THCA. To this end, we constructed a prognostic model centered on a core set of NRGs, capable of stratifying patients by risk and predicting clinical outcomes.

Material and methods

Data collection

The mRNA expression profiles, somatic mutation data, and corresponding clinical information for patients with THCA were retrieved from the TCGA database (https://portal.gdc.cancer.gov). To ensure the reliability of the survival analyses, patients with an overall survival (OS) of less than 30 days were excluded, resulting in a final cohort of 502 cases. The cohort was randomly partitioned into a training and a validation set (7 : 3). Neutrophil-related genes were curated based on previously published literature and further refined using gene sets annotated in the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb).

Identification of differentially expressed NRGs

Differentially expressed genes (DEGs) between tumor and normal samples in the TCGA-THCA cohort were identified using the limma package in R (https://bioconductor.org/packages/release/bioc/html/limma.html) (|log2FC| > 1, adjusted p < 0.05). The intersection of the DEGs and NRGs was then computed to obtain differentially expressed NRGs (DENRGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on the DENRGs using the clusterProfiler R package (https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html).

Selection of prognostic features and construction of the THCA prognostic model

Univariate Cox regression was performed on the TCGA-THCA cohort to identify candidate genes (p < 0.05). Least absolute shrinkage and selection operator (LASSO) regression, implemented via the glmnet package in R (https://cran.r-project.org/web/packages/glmnet/index.html), was used to minimize overfitting and select genes by cross-validation with an optimal penalty parameter (lambda). Multivariate Cox regression was then applied to construct a prognostic model using the survival package (https://cran.r-project.org/web/packages/survival/index.html). Risk scores were calculated based on gene expression levels and the corresponding coefficients. Patients were stratified into 2 groups according to the median risk score. Survival curves were plotted, and the model’s predictive accuracy was evaluated using time-dependent receiver operating characteristic (ROC) curves with the timeROC package (https://cran.r-project.org/web/packages/timeROC/index.html). Risk score distributions, survival status, and heatmaps of gene expression were generated for visual comparison between the groups.

Construction of a nomogram for independent prognostic analysis

To assess the independent prognostic value of the constructed risk model, a nomogram was developed using the rms package in R (https://cran.r-project.org/web/packages/rms/index.html), integrating the prognostic risk score and clinical factors to predict the survival probabilities of patients with THCA. Calibration curves were plotted to evaluate the accuracy of the nomogram’s predictions by comparing the predicted and observed survival rates.

Immune infiltration analysis

Malignant Tumours using Expression data (ESTIMATE) scores for patients with TCGA-THCA were calculated using Single Sample Gene Set Enrichment Analysis (ssGSEA) with the estimate and GSVA R packages (https://estimate.r-forge.r-project.org; https://www.bioconductor.org/packages/release/bioc/html/GSVA.html). Immune infiltration levels were assessed using the CIBERSORT algorithm (https://cibersortx.stanford.edu). Immune checkpoint expression across different groups was analyzed and visualized using box plots. The association between immune subtypes and risk groups (high vs low) was evaluated. Immune subtype scoring data for THCA were downloaded from The Cancer Immunome Atlas (TCIA), and differences in scores between subtypes were analyzed.

Gene set enrichment analysis

The gene set enrichment analysis (GSEA) was performed using the GSEA software (https://www.gsea-msigdb.org/gsea/downloads.jsp) to identify enriched pathways between risk groups. Differential analysis between the high- and low-risk groups was conducted using the limma package (|log2FC| > 0.585, p < 0.05). The GO enrichment analysis was performed on differential genes.

THCA subtype analysis

Based on the NRGs from the constructed model, consensus clustering was performed on THCA samples from the TCGA-THCA dataset using the ConsensusClusterPlus package (https://www.rdocumentation.org/packages/ConsensusClusterPlus/versions/1.36.0/topics/ConsensusClusterPlus). A total of 205 patients were assigned to Cluster 1 and 297 to Cluster 2. Kaplan–Meier survival curves were generated using the survminer package (https://cran.r-project.org/web/packages/survminer/index.html) to compare OS between subtypes. Immune infiltration analysis was subsequently conducted for both subtypes, with differential analysis performed using the limma package and GO enrichment analysis.

Drug sensitivity analysis and regulatory network construction

To identify potential therapeutic targets and effective drugs, we utilized the CellMiner database (https://discover.nci.nih.gov/cellminercdb) to screen for antitumor drugs whose sensitivity is significantly associated with prognostic genes. The pRRophetic R package (https://github.com/paulgeeleher/pRRophetic) was used to predict the half-maximal inhibitory concentration (IC50) of different drugs in risk groups.

Statistical analyses

Statistical significance was determined using R v 4.4.2 (R Foundation for Statistical Computing, Vienna, Austria). Spearman analysis was utilized to compute the correlation coefficient. The difference between groups was assessed for statistical significance using Wilcoxon rank sum test. Two-way analysis of variance (ANOVA) test was utilized for comparing differences among the 3 groups. R v. 4.4.2 software was employed for analyzing all bioinformatics data. All Cox regression analyses, including model fitting and assumption verification, were performed using the survival R package. The ROC curves were used to evaluate the predictive performance of candidate genes used to construct predictive models.

Results

Development of a neutrophil-related prognostic model for THCA

The overall research workflow is illustrated in Figure 1. Differential expression analysis between tumor and normal tissues identified DEGs (Supplementary Fig. 1A). Intersecting the DEGs with the NRGs gen set yielded 91 DENRGs (Supplementary Table 1), which were subjected to GO and KEGG enrichment analyses. These genes were significantly enriched in immune-related processes and pathways, including “leukocyte migration”, “cell chemotaxis”, “cytoplasmic vesicle lumen”, and “cytokine–cytokine receptor interaction” (Supplementary Fig. 1B). Univariate Cox regression identified prognostic DENRGs (Figure 2A), followed by LASSO regression to reduce redundancy and prevent overfitting, ultimately selecting 4 candidate genes (Figure 2B,C). Multivariate Cox regression refined the model to 3 signature genes (Figure 2D): Risk score = 0.583 × CCL17 – 0.467 × NPC2 – 0.564 × DPP4. Schoenfeld and Martingale residuals indicated that the proportional hazards and linearity assumptions of the Cox model were satisfied, and all variance inflation factors were below 3 (Supplementary Tables 2,3 and Supplementary Fig. 2), supporting the validity of the analysis. Spearman’s correlation among the 3 genes was low (all p < 0.05, Supplementary Fig. 3), confirming the independence of the predictors.

Patients were stratified into risk groups based on the median risk score. Time-dependent ROC analysis showed that the area under the curve (AUC) values for 1-, 3-, and 5-year survival all exceeded 0.7 in the training, validation, and full TCGA cohorts, suggesting strong predictive performance (Figure 3A). Compared with other transcriptome-based prognostic models, our model achieved relatively high AUC values (Supplementary Fig. 4). Kaplan–Meier analysis demonstrated significantly worse survival in the high-risk group across all cohorts (Figure 3B). Risk score distributions, survival status plots, and heatmaps of gene expression further illustrated clear distinctions between the 2 risk groups (Figure 3C). Among the 3 signature genes, CCL17 was significantly upregulated in the high-risk group, whereas NPC2 and DPP4 were downregulated (Figure 3D). Expression levels of all 3 genes were significantly elevated in tumor tissues compared with normal controls (Figure 3E). Moreover, survival analysis based on gene expression stratification revealed that high NPC2 expression was associated with a better prognosis (Figure 3F).

Independent prognostic value and nomogram construction

Univariate and multivariate Cox regression analyses, including clinical parameters and the risk score, indicated that both age and the risk score were independently associated with OS (Figure 4A,B). Notably, the risk score remained a significant adverse prognostic factor (p < 0.01, hazard ratio (HR) >1), confirming its independence and robustness. A prognostic nomogram integrating clinical variables and the risk score was developed to predict survival probabilities in patients with THCA (Figure 4C). Calibration plots demonstrated strong concordance between the predicted and observed survival outcomes (Figure 4D,E), indicating good predictive performance of the model. Moreover, patients with advanced-stage disease (stage III–IV) exhibited significantly higher risk scores, aligning with the stratified prognostic potential of the risk model (Supplementary Fig. 5).

Immune infiltration landscape between high- and low-risk groups

Immune phenotyping based on the Immunophenoscore (IPS) revealed significantly higher scores across most immune categories in the low-risk group, except for the IPS_CTLA4_pos_PD1_pos subgroup, which showed no significant difference (Figure 5A, Supplementary Table 4). These results suggest enhanced predicted responsiveness to immune checkpoint blockade in the low-risk cohort. We next assessed the expression of immune checkpoint (ICP)-related genes between the 2 groups (Figure 5B, Supplementary Table 5). Most ICP genes, including CD200, CD11, and TNFSF15, were downregulated in the high-risk group. Notably, IDO1 and TNFSF4 were significantly upregulated in the high-risk group, suggesting potential immune evasion mechanisms. Although the high-risk group exhibited a significantly elevated StromalScore, analysis of the IMvigor210 cohort revealed a reduced immune response compared with the low-risk group (Figure 5C, Supplementary Tables 6,7). ssGSEA revealed distinct immune functional differences: aDCs and type II IFN response were enriched in the low-risk group, whereas B cells, dendritic cells (DCs), and T follicular helper (Tfh) cells were more abundant in the high-risk group (Figure 6A,B, Supplementary Table 8). To further dissect immune infiltration, we applied the CIBERSORT algorithm to estimate the abundance of 22 immune cell types. Significant differences were observed in specific populations (Figure 6C, Supplementary Table 9): monocytes were increased in the high-risk group, whereas M0 macrophages and resting mast cells were more enriched in the low-risk group.

Pathway enrichment analysis with GSEA

Gene set enrichment analysis was performed to identify KEGG pathways differentially enriched between the risk groups. In the high-risk group, significantly enriched pathways included cardiac muscle contraction, oxidative phosphorylation, Parkinson’s disease, and propanoate metabolism (Supplementary Fig. 6). In contrast, the low-risk group showed enrichment in pathways such as adherens junction, chronic myeloid leukemia, DNA replication, and lysosome (Supplementary Fig. 7). Subsequent differential expression analysis between the 2 risk groups (|log2FC| > 0.585, adjusted p < 0.05) identified 67 upregulated and 176 downregulated genes. Gene Ontology enrichment analysis revealed that the upregulated genes were primarily associated with response to copper ion and response to cadmium ion (BP), basement membrane (CC), and ECM structural constituent and heparin binding (MF) (Figure 7A). Conversely, downregulated genes were enriched in epidermis development, wound healing, and axonogenesis (BP); collagen-containing ECM and cytoplasmic vesicle lumen (CC); and enzyme inhibitor activity, glycosaminoglycan binding, and peptidase regulator activity (MF) (Figure 7B).

Molecular subtype analysis of THCA based on prognostic genes

Unsupervised consensus clustering based on the expression of the 3 prognostic genes identified 2 molecular subtypes of THCA: Cluster 1 (n = 205) and Cluster 2 (n = 297) (Figure 8A). Kaplan–Meier survival analysis revealed significantly poorer OS in Cluster 1 (Figure 8B). Expression levels of the 3 genes differed markedly between the 2 subtypes (Figure 8C and Supplementary Table 10), with Cluster 1 exhibiting a significantly higher risk score (Figure 8D and Supplementary Table 11). Immune profiling using the CIBERSORT algorithm indicated distinct immune infiltration patterns between the subtypes. Notably, monocytes were enriched in Cluster 1, whereas M0 macrophages were more abundant in Cluster 2 (Figure 8E and Supplementary Table 12). ssGSEA revealed that Cluster 1 exhibited generally higher immune cell infiltration and immune-related functional scores than Cluster 2 (Figure 8F and Supplementary Table 13). Moreover, tumor purity was higher in Cluster 2, whereas stromal and ESTIMATE scores were significantly higher in Cluster 1 (Figure 9A and Supplementary Table 14). Analysis of the relationship between the subtypes and the 6 immune subtypes (C1–C6) from the TCGA Pan-Cancer immune classification revealed that both THCA clusters were predominantly enriched in the inflammatory subtype C3 (Figure 9B). Immune checkpoint gene expression analysis showed significantly higher levels of CD40, CD44, CD200, and TNFSF15 in Cluster 2, whereas CCL19, GZMB, CTLA4, and IDO1 were more abundant in Cluster 1 (Figure 9C and Supplementary Table 15). Differential expression analysis between the 2 subtypes (|log2FC| > 0.585, adjusted p < 0.05) identified 189 DEGs, including 47 upregulated and 142 downregulated genes in Cluster 1. Gene Ontology enrichment analysis of the upregulated genes highlighted terms such as Wnt-protein binding and CCR chemokine receptor binding, whereas the downregulated genes were enriched in enzyme inhibitor activity (Figure 9D,E). Tumor mutational burden (TMB) analysis showed that the TMB of Cluster 1 was significantly higher than that of Cluster 2 (Supplementary Fig. 8).

Drug sensitivity analysis and network construction

We estimated the IC50 values of commonly used chemotherapeutic agents (docetaxel, cisplatin, fluorouracil, vinorelbine, doxorubicin, and paclitaxel) in the risk groups. IC50 was significantly higher in the high-risk group, suggesting reduced drug sensitivity (Figure 10 and Supplementary Table 16). To further explore potential therapeutic targets, we analyzed the correlation between prognostic gene expression and drug sensitivity using data from the CellMiner database. DPP4 expression was positively correlated with sensitivity to Midostaurin and ENMD-2076 precursor (r = 0.523 and 0.494, respectively). CCL17 showed a significant negative correlation with 7-hydroxystaurosporine (r = –0.450) and a positive correlation with PLX-4720 (r = 0.465). NPC2 expression was negatively correlated with des-fluoro-TAK-960 and docetaxel (r = –0.407 and –0.437, respectively) (Supplementary Fig. 9 and Supplementary Table 17).

Discussion

Through comprehensive bioinformatics analysis of neutrophil-related genes, this study uncovers novel insights of THCA. It highlights the intricate relationship between neutrophils and THCA progression. Neutrophils are critical regulators of tumor progression. Our research confirms existing literature by demonstrating that high NRGs scores are associated with poor prognosis in THCA patients.15, 16 The neutrophil-to-lymphocyte ratio (NLR) has prognostic value, with elevated NLR correlating with significantly worse OS.17 Inflammatory stimuli induce the release of neutrophil extracellular traps (NETs), which have been implicated in promoting tumor growth and progression.18 Additionally, SLCO4A1 may effect progression-free survival (PFS) in THCA patients by modulating neutrophil-mediated immune pathways.19

The 3 core NRGs (CCL17, NPC2, and DPP4) are associated with THCA has been established. CCL17, a CCR4 ligand expressed in M2 macrophages and thymus,20, 21 has been implicated in tumor progression across cancers.22, 23, 24, 25 Consistent with previous findings in THCA,26 our data show that high CCL17 expression is associated with poorer survival. NPC2 which encodes a lysosomal cholesterol transporter,27 has been implicated in several malignancies including colon,28 lung cancer,29, 30 and thyroid cancer.31 In THCA, NPC2 knockdown promotes tumor cell apoptosis in vitro. DPP4 (CD26) is a cell surface glycoprotein.32 DPP4 acts as either an oncogene or tumor suppressor depending on context.33 Its high expression is linked to poor prognosis in clear cell carcinoma,34 pancreatic cancer,35 and colorectal cancer,36 while exhibiting tumor-suppressive roles in ovarian, endometrial,37 and prostate cancer.38 DPP4 overexpression in THCA can be used as a potential prognostic marker.39, 40 Compared to models from other transcriptomic studies, our risk model demonstrated superior predictive performance, with AUCs ranging from 0.797 to 0.933. While direct comparisons with genomic markers such as BRAF mutation status are constrained by differences in data type and clinical context, our integrated nomogram – combining the risk score with key clinical parameters – offers a practical and interpretable tool for personalized patient management.

We analyzed DEGs between risk groups to explore how risk scores influence THCA progression. High-risk tumors showed enrichment of oxidative phosphorylation and myocardial contraction pathways, consistent with cardiovascular complications often observed in aggressive THCA under thyrotoxic conditions.41 OXPHOS also promotes resistance to apoptosis and persistence of TH17 cells in the TME.42 Co-upregulation of copper ion response and collagen remodeling genes indicates a pro-metastatic phenotype, as copper drives tumor growth via RTK/ERK signaling and mediates cuproptosis.43, 44 Increased matrix stiffness may trigger matrix autophagy, further supporting tumor progression.45 These findings reflect the connective tissue proliferative stroma and metabolic plasticity observed in advanced THCA.46 In contrast, low-risk tumors preserve adherens junctions and DNA replication fidelity.47 Suppression of axonogenesis and epidermal development genes may reduce perineural invasion.48, 49 These findings provide mechanistic rationale for exploring targeted or repurposed therapies in aggressive THCA. Notably, bergapten may serve as a natural therapeutic agent for papillary thyroid carcinoma (PTC), as it has been shown to inhibit PTC cell growth through suppression of PI3K/AKT and GSK-3β signaling, highlighting the potential of pathway-specific interventions.50

Previous studies demonstrate that dynamic interactions within the TME drive tumor progression.51, 52 In this study, we systematically profiled the immune landscape of THCA by integrating risk stratification and molecular subtyping. Immunophenoscore analysis revealed significantly higher immune scores in the low-risk group across most categories, except for the CTLA4+PD1+ subgroup. This corresponded with higher expression of most immune checkpoint genes, suggesting greater immunotherapy responsiveness in this group. These results are aligned with prior observations that immunotherapy can improve outcomes in THCA patients.53 Moreover, we uncovered distinct immune infiltration patterns linked to our novel molecular subtypes. Cluster 1 exhibited higher infiltration of M0 macrophages, along with elevated CTLA4 expression and higher TMB. These features suggest an immunologically active phenotype potentially responsive to immune checkpoint inhibitors, particularly CTLA4 blockade.54, 55 In contrast, Cluster 2 displayed lower immune infiltration, reduced TMB, and frequent BRAF mutations, consistent with an immune-cold, genomically stable phenotype that may benefit more from mitogen-activated protein kinase (MAPK)-targeted therapies.56 While previous THCA immune profiling studies have highlighted macrophage and T-cell infiltration,39, 57 our subtype-specific immune landscapes – especially the integration of neutrophil-related risk stratification – have not been previously reported. This framework provides novel insights into THCA immunobiology and may aid in identifying patients more likely to benefit from immunotherapy.

High-risk THCA patients exhibited broadly elevated IC50 values, suggesting low-risk group is more sensitive to chemotherapy drugs than the high-risk group. This pattern may reflect intrinsic resistance mechanisms and underscores the limited benefit of cytotoxic regimens in aggressive THCA. To inform alternative strategies, we integrated gene expression with CellMiner drug response data to identify actionable associations. DPP4 expression correlated with sensitivity to midostaurin and ENMD-2076, suggesting vulnerability to multi-kinase inhibition.58 CCL17 expression was associated with opposing sensitivities to 7-hydroxystaurosporine and PLX-4720, reflecting pathway-specific therapeutic windows.59, 60 NPC2 was linked to resistance against des-fluoro-TAK-960 and docetaxel, implying a role in mitotic drug evasion.61 These associations highlight candidate compounds for drug repurposing in high-risk THCA and warrant further translational exploration.

Limitations of the study

Despite the comprehensive analysis of NRGs in THCA, several limitations remain. Owing to current experimental constraints, functional and drug response validation could not be performed and will be pursued in future studies to support clinical translation. Future application of single-cell RNA sequencing and spatial transcriptomics could uncover spatially confined immune niches or cell–cell interactions shaping the TME. Moreover, functional validation of the prognostic genes and subtype-specific immune features through in vitro perturbation and in vivo tumor models will be essential to confirm their mechanistic relevance and therapeutic potential.

Conclusions

This study elucidates the role of NRGs in shaping the THCA TME and proposes a clinically applicable biomarker system for prognostic stratification and therapeutic guidance. The 3-gene NRG signature predicts patient survival and facilitates risk-adapted treatment strategies involving chemotherapy and immunotherapy. Incorporating the signature into a nomogram alongside clinical variables enhances prognostic accuracy and supports its translational applicability. This model may assist in stratifying patients according to treatment intensity and tailoring follow-up schedules based on individualized risk profiles. Furthermore, subtype-specific immune landscapes underscore actionable therapeutic opportunities and may inform individualized approaches to immunotherapy and drug repurposing.

Supplementary data

The supplementary materials are available at https://doi.org/10.5281/zenodo.17730267. The package contains the following files:

Supplementary Fig. 1. Differential gene expression and enrichment analysis.

Supplementary Fig. 2. Martingale residuals of multivariate Cox model.

Supplementary Fig. 3. Spearman correlation analysis of 3 signature genes.

Supplementary Fig. 4. Comparison of AUC values for multiple models.

Supplementary Fig. 5. Risk scores for different cancer stages.

Supplementary Fig. 6. GSEA results of high-risk groups.

Supplementary Fig. 7. GSEA results of low-risk groups.

Supplementary Fig. 8. TMB of subtype.

Supplementary Fig. 9. Correlation plot of drug sensitivity predictions using the CellMiner database.

Supplementary Table 1. Differentially expressed neutrophil-related genes identified between tumor and normal samples.

Supplementary Table 2. Schoenfeld residual test for proportional hazards assumption.

Supplementary Table 3. Variance inflation factor (VIF) for multicollinearity. assessment in the Cox regression model.

Supplementary Table 4. Comparison of immune phenotype scores between risk groups using Wilcoxon test.

Supplementary Table 5. Comparison of immune checkpoint–related gene expression between high- and low-risk groups using Wilcoxon test.

Supplementary Table 6. Comparison of ESTIMATE scores between high- and low-risk groups using Wilcoxon test.

Supplementary Table 7. χ2 test results for the distribution of no response (NR) and response (R) across high- and low-risk groups.

Supplementary Table 8. Wilcoxon test of immune cell composition and immune function between high- and low-risk groups.

Supplementary Table 9. Comparison of immune cell type proportions between high- and low-risk groups using Wilcoxon test.

Supplementary Table 10. Wilcoxon test of signature gene expression between clusters.

Supplementary Table 11. Wilcoxon test of risk scores between clusters.

Supplementary Table 12. Wilcoxon test of immune cell composition between clusters.

Supplementary Table 13. Wilcoxon test of immune cell composition and immune function between clusters.

Supplementary Table 14. Wilcoxon test of ESTIMATE scores and tumor purity between clusters.

Supplementary Table 15. Wilcoxon test of immune checkpoint gene expression between clusters.

Supplementary Table 16. Wilcoxon test of drug response between risk groups.

Supplementary Table 17. Gene–drug correlation analysis.

Data Availability Statement

The datasets supporting the findings of the current study are openly available in Zenodo at https://doi.org/10.5281/zenodo.17734513.

Consent for publication of personal information

Not applicable.

Use of AI and AI-assisted technologies

Not applicable.

Figures


Fig. 1. Flowchart of the research process
THCA – thyroid carcinoma; TCGA – The Cancer Genome Atlas; KEGG – Kyoto Encyclopedia of Genes and Genomes; GO – Gene Ontology; GSEA – gene set enrichment analysis; DEGs – differentially expressed genes; NRGs – neutrophil-related genes; DENRGs – differentially expressed NRGs; ROC – receiver operating characteristic; LASSO – least absolute shrinkage and selection operator.
Fig. 2. Development and validation of the prognostic model. A. Univariate Cox regression analysis of prognosis-related differentially expressed neutrophil-related genes (DENRGs); B. Feature selection using least absolute shrinkage and selection operator (LASSO) regression analysis of candidate neutrophil-related genes (NRGs); C. Coefficient profiles from LASSO regression of the 4 prognostic NRGs; D. Multivariate Cox regression results for the 3-gene NRG signature
HR – hazard ratio; 95% CI – 95% confidence interval.
Fig. 3. Prognostic performance and expression patterns of the neutrophil-related gene (NRG)-based signature. A. Receiver operating characteristic (ROC) analysis in train, test, and The Cancer Genome Atlas (TCGA)-thyroid carcinoma (THCA) cohort; B. Risk score distribution, and survival status stratified by risk score for the training, test and TCGA-THCA cohort; C. The difference in expression of the 3 NRGs between high- and low-risk in TCGA dataset; D. The difference in expression of the 3 NRGs between normal and tumor tissues in TCGA dataset; E. Signature gene expression across risk groups
*p < 0.050; **p < 0.010; ***p < 0.001; ****p < 0.0001; ns – not significant.
Fig. 4. Comprehensive prognostic evaluation of 3 signature genes and nomogram in thyroid carcinoma (THCA). Forest plot of univariate Cox (A) and (B) multivariate Cox regression analysis; C. Nomogram for predicting 1-, 3-, and 5-year survival; D. Decision curve analysis for 1-, 3-, and 5-year overall survival (OS); E. Calibration plots for 1-year, 3-year, and 5-year survival rates
Fig. 5. Immune landscape and stromal characteristics in the risk subgroup. A. Immunophenoscore (IPS); B. Risk-stratified expression patterns of immune checkpoint molecules; C. ESTIMATE algorithm-derived stromal scores comparison
* p < 0.050; ** p < 0.010; *** p < 0.001; ns – not significant. NR – no response; R – response.
Fig. 6. Immune microenvironment and infiltration analysis in risk subgroups. A. The heatmap showed the immune cell infiltration and immune functions; B. Box plots of Single Sample Gene Set Enrichment Analysis (ssGSEA) immune checkpoint scores; C. CIBERSORT algorithm analysis of immune infiltration
* p < 0.050; ** p < 0.010; *** p < 0.001; ns – not significant.
Fig. 7. Enrichment analysis between risk groups. A. Gene Ontology (GO) enrichment of upregulated genes in the high-risk group; B. GO enrichment of downregulated genes in the low-risk group
Fig. 8. Molecular subtyping and associated immune characteristics. A. Consensus clustering analysis heatmap with 2 subgroups (k = 2); B. Kaplan–Meier survival curve among subtypes; C. Differential expression of signature genes across subtypes; D. Boxplot of risk scores between subtypes; E. CIBERSORT algorithm analysis of immune infiltration; F. Boxplots displaying difference of immune cell infiltration profiles between clusters
*p < 0.050; **p < 0.010; ***p < 0.001; ****p < 0.0001; ns –not significant.
Fig. 9. Immune landscape and functional enrichment across molecular subtypes. A. Boxplots of stromal immune, ESTIMATE scores, and tumor purity among subtypes; B. The interaction between immune subtypes and 2 risk groups was depicted as an alluvial diagram; C. Immune checkpoint gene expression among subtypes; D. Upregulated Gene Ontology (GO) enriched bubble chart; E. Downregulated GO enriched bubble chart
*p < 0.050; **p < 0.010; ***p < 0.001; ****p < 0.0001; ns –not significant. C1 – wound healing, C2 – interferon gamma (IFN-γ)-dominant, C3 – inflammatory, and C4 – lymphocyte-depleted; C6 – tumor growth factor beta (TGF-β)-dominant.
Fig. 10. Drug sensitivity prediction in thyroid carcinoma (THCA). Violin plot of half maximal inhibitory concentration (IC50) values in risk groups
*p < 0.050; **p < 0.010; ***p < 0.001; ****p < 0.0001; ns – not significant.

References (61)

  1. Shin E, Koo JS. Cell component and function of tumor microenvironment in thyroid cancer. Int J Mol Sci. 2022;23(20):12578. doi:10.3390/ijms232012578
  2. Kitahara CM, Sosa JA. The changing incidence of thyroid cancer. Nat Rev Endocrinol. 2016;12(11):646–653. doi:10.1038/nrendo.2016.110
  3. Cherifi F, Awada A. Molecular oncology of iodine refractory thyroid cancer current therapies and perspective. Crit Rev Oncol Hematol. 2025;209:104679. doi:10.1016/j.critrevonc.2025.104679
  4. Kim J, Gosnell JE, Roman SA. Geographic influences in the global rise of thyroid cancer. Nat Rev Endocrinol. 2020;16(1):17–29. doi:10.1038/s41574-019-0263-x
  5. Liao D, Yang G, Yang Y, et al. Identification of pannexin 2 as a novel marker correlating with ferroptosis and malignant phenotypes of prostate cancer cells. Onco Targets Ther. 2020;13:4411–4421. doi:10.2147/OTT.S249752
  6. Schlumberger M, Leboulleux S. Current practice in patients with differentiated thyroid cancer. Nat Rev Endocrinol. 2021;17(3):176–188. doi:10.1038/s41574-020-00448-z
  7. Lorusso L, Cappagli V, Valerio L, et al. Thyroid cancers: From surgery to current and future systemic therapies through their molecular identities. Int J Mol Sci. 2021;22(6):3117. doi:10.3390/ijms22063117
  8. Guan Z, Wang H, Tian M. A cuproptosis-related gene signature as a prognostic biomarker in thyroid cancer based on transcriptomics. Biochem Genet. 2025;63(2):1584–1604. doi:10.1007/s10528-024-10767-9
  9. Coca-Pelaz A, Shah JP, Hernandez-Prera JC, et al. Papillary thyroid cancer: Aggressive variants and impact on management. A narrative review. Adv Ther. 2020;37(7):3112–3128. doi:10.1007/s12325-020-01391-1
  10. Yi M, Li T, Niu M, et al. Exploiting innate immunity for cancer immunotherapy. Mol Cancer. 2023;22(1):187. doi:10.1186/s12943-023-01885-w
  11. Li C, Yu X, Han X, et al. Innate immune cells in tumor microenvironment: A new frontier in cancer immunotherapy. iScience. 2024;27(9):110750. doi:10.1016/j.isci.2024.110750
  12. Gonzalez H, Hagerling C, Werb Z. Roles of the immune system in cancer: From tumor initiation to metastatic progression. Genes Dev. 2018;32(19–20):1267–1284. doi:10.1101/gad.314617.118
  13. Szczerba BM, Castro-Giner F, Vetter M, et al. Neutrophils escort circulating tumour cells to enable cell cycle progression. Nature. 2019;566(7745):553–557. doi:10.1038/s41586-019-0915-y
  14. Wu CF, Andzinski L, Kasnitz N, et al. The lack of type I interferon induces neutrophil-mediated pre-metastatic niche formation in the mouse lung: Type I IFN suppresses metastasis formation. Int J Cancer. 2015;137(4):837–847. doi:10.1002/ijc.29444
  15. Bozan MB, Yazar FM, Kale İT, Yüzbaşıoğlu MF, Boran ÖF, Azak Bozan A. Delta neutrophil index and neutrophil-to-lymphocyte ratio in the differentiation of thyroid malignancy and nodular goiter. World J Surg. 2021;45(2):507–514. doi:10.1007/s00268-020-05822-6
  16. Galdiero MR, Varricchi G, Loffredo S, et al. Potential involvement of neutrophils in human thyroid cancer. PLoS One. 2018;13(6):e0199740. doi:10.1371/journal.pone.0199740
  17. Gao Q, Quan M, Zhang L, Ran Y, Zhong J, Wang B. Neutrophil-to-lymphocyte ratio as a prognostic indicator in thyroid cancer. Cancer Control. 2024;31:10732748241309048. doi:10.1177/10732748241309048
  18. Lee H, Kim TS, Gu JY, et al. Value of circulating neutrophil elastase for detecting recurrence of differentiated thyroid cancer. Endocr Connect. 2023;12(12):e230400. doi:10.1530/EC-23-0400
  19. Wang XS, Wu SL, Peng Z, Zhu HH. SLCO4A1 is a prognosis-associated biomarker involved in neutrophil-mediated immunity in thyroid cancer. Int J Gen Med. 2021;14:9615–9628. doi:10.2147/IJGM.S339921
  20. Feng G, Bajpai G, Ma P, et al. CCL17 aggravates myocardial injury by suppressing recruitment of regulatory T cells. Circulation. 2022;145(10):765–782. doi:10.1161/CIRCULATIONAHA.121.055888
  21. Zhang A, Xu Y, Xu H, et al. Lactate-induced M2 polarization of tumor-associated macrophages promotes the invasion of pituitary adenoma by secreting CCL17. Theranostics. 2021;11(8):3839–3852. doi:10.7150/thno.53749
  22. Zhu F, Li X, Chen S, Zeng Q, Zhao Y, Luo F. Tumor-associated macrophage or chemokine ligand CCL17 positively regulates the tumorigenesis of hepatocellular carcinoma. Med Oncol. 2016;33(2):17. doi:10.1007/s12032-016-0729-9
  23. Liu LB, Xie F, Chang KK, et al. Chemokine CCL17 induced by hypoxia promotes the proliferation of cervical cancer cell. Am J Cancer Res. 2015;5(10):3072–3084. PMID:26693060. PMCID:PMC4656731.
  24. Huo X, Zhou X, Peng P, et al. Identification of a six-gene signature for predicting the overall survival of cervical cancer patients. Onco Targets Ther. 2021;14:809–822. doi:10.2147/OTT.S276553
  25. Ye T, Zhang X, Dong Y, et al. Chemokine CCL17 affects local immune infiltration characteristics and early prognosis value of lung adenocarcinoma. Front Cell Dev Biol. 2022;10:816927. doi:10.3389/fcell.2022.816927
  26. Gu X, Chen B, Zhang S, Zhai X, Hu Y, Ye H. The expression of CCL17 and potential prognostic value on tumor immunity in thyroid carcinoma based on bioinformatics analysis. Sci Rep. 2024;14(1):31580. doi:10.1038/s41598-024-75750-1
  27. Chae YS, Kim H. NPC2 expression in thyroid tumors and its possible diagnostic utility. Int J Clin Exp Pathol. 2021;14(1):126–132. PMID:33532030. PMCID:PMC7847500.
  28. Liao YJ, Lin MW, Yen CH, et al. Characterization of Niemann–Pick type C2 protein expression in multiple cancers using a novel NPC2 monoclonal antibody. PLoS One. 2013;8(10):e77586. doi:10.1371/journal.pone.0077586
  29. Wang ZM, Ning ZL, Ma C, Liu TB, Tao B, Guo L. Low expression of lysosome-related genes KCNE1, NPC2, and SFTPD promote cancer cell proliferation and tumor associated M2 macrophage polarization in lung adenocarcinoma. Heliyon. 2024;10(6):e27575. doi:10.1016/j.heliyon.2024.e27575
  30. Robles J, Pintado-Berninches L, Boukich I, et al. A prognostic six-gene expression risk-score derived from proteomic profiling of the metastatic colorectal cancer secretome. J Pathol Clin Res. 2022;8(6):495–508. doi:10.1002/cjp2.294
  31. Asakawa JI, Kodaira M, Ishikawa N, et al. Two-dimensional complementary deoxyribonucleic acid electrophoresis revealing upregulated human epididymal protein-1 and downregulated CL-100 in thyroid papillary carcinoma. Endocrinology. 2002;143(11):4422–4428. doi:10.1210/en.2002-220550
  32. Sun L, Ma Y, Geng C, et al. DPP4, a potential tumor biomarker, and tumor therapeutic target: Review. Mol Biol Rep. 2025;52(1):126. doi:10.1007/s11033-025-10235-6
  33. Gao X, Le Y, Geng C, Jiang Z, Zhao G, Zhang P. DPP4 is a potential prognostic marker of thyroid carcinoma and a target for immunotherapy. Int J Endocrinol. 2022;2022:5181386. doi:10.1155/2022/5181386
  34. Zheng B, Niu Z, Si S, et al. Comprehensive analysis of new prognostic signature based on ferroptosis-related genes in clear cell renal cell carcinoma. Aging (Albany NY). 2021;13(15):19789–19804. doi:10.18632/aging.203390
  35. Yan L, Tian X, Ye C, et al. CD26 as a promising biomarker for predicting prognosis in patients with pancreatic tumors. Onco Targets Ther. 2020;13:12615–12623. doi:10.2147/OTT.S278736
  36. Ng L, Wong SKM, Huang Z, et al. CD26 induces colorectal cancer angiogenesis and metastasis through CAV1/MMP1 signaling. Int J Mol Sci. 2022;23(3):1181. doi:10.3390/ijms23031181
  37. Niazmand A, Nedaeinia R, Vatandoost N, et al. The impacts of dipeptidyl-peptidase 4 (DPP-4) inhibitors on common female malignancies: A systematic review. Gene. 2024;927:148659. doi:10.1016/j.gene.2024.148659
  38. Wen Y, Lin C, Hsiao C, et al. Genetic variants of dipeptidyl peptidase IV are linked to the clinicopathologic development of prostate cancer. J Cell Mol Med. 2023;27(17):2507–2516. doi:10.1111/jcmm.17845
  39. Ren X, Du H, Cheng W, et al. Construction of a ferroptosis-related eight gene signature for predicting the prognosis and immune infiltration of thyroid cancer. Front Endocrinol (Lausanne). 2022;13:997873. doi:10.3389/fendo.2022.997873
  40. Cheng SY, Wu ATH, Batiha GES, et al. Identification of DPP4/CTNNB1/MET as a theranostic signature of thyroid cancer and evaluation of the therapeutic potential of sitagliptin. Biology (Basel). 2022;11(2):324. doi:10.3390/biology11020324
  41. Ertek S, Cicero AF. State of the art paper hyperthyroidism and cardiovascular complications: A narrative review on the basis of pathophysiology. Arch Med Sci. 2013;5:944–952. doi:10.5114/aoms.2013.38685
  42. Hong HS, Mbah NE, Shan M, et al. OXPHOS promotes apoptotic resistance and cellular persistence in TH17 cells in the periphery and tumor microenvironment. Sci Immunol. 2022;7(77):eabm8182. doi:10.1126/sciimmunol.abm8182
  43. Feng Y, Yang Z, Wang J, Zhao H. Cuproptosis: Unveiling a new frontier in cancer biology and therapeutics. Cell Commun Signal. 2024;22(1):249. doi:10.1186/s12964-024-01625-7
  44. Tsvetkov P, Coy S, Petrova B, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science. 2022;375(6586):1254–1261. doi:10.1126/science.abf0529
  45. Hupfer A, Brichkina A, Koeniger A, et al. Matrix stiffness drives stromal autophagy and promotes formation of a protumorigenic niche. Proc Natl Acad Sci U S A. 2021;118(40):e2105367118. doi:10.1073/pnas.2105367118
  46. Cox TR. The matrix in cancer. Nat Rev Cancer. 2021;21(4):217–238. doi:10.1038/s41568-020-00329-7
  47. Riesco-Eizaguirre G, Rodríguez I, De La Vieja A, et al. The BRAFV600E oncogene induces transforming growth factor β secretion leading to sodium iodide symporter repression and increased malignancy in thyroid cancer. Cancer Res. 2009;69(21):8317–8325. doi:10.1158/0008-5472.CAN-09-1248
  48. Wang H, Zheng Q, Lu Z, et al. Role of the nervous system in cancers: A review. Cell Death Discov. 2021;7(1):76. doi:10.1038/s41420-021-00450-y
  49. Revilla G, Al Qtaish N, Caruana P, et al. Lenvatinib-loaded poly(lactic-co-glycolic acid) nanoparticles with epidermal growth factor receptor antibody conjugation as a preclinical approach to therapeutically improve thyroid cancer with aggressive behavior. Biomolecules. 2023;13(11):1647. doi:10.3390/biom13111647
  50. Li J, Hussain SA, Rayalu Daddam J, Sun M. Bergapten attenuates human papillary thyroid cancer cell proliferation by triggering apoptosis and the GSK-3β, P13K and AKT pathways. Adv Clin Exp Med. 2024;34(1):113–122. doi:10.17219/acem/183877
  51. Fu T, Dai LJ, Wu SY, et al. Spatial architecture of the immune microenvironment orchestrates tumor immunity and therapeutic response. J Hematol Oncol. 2021;14(1):98. doi:10.1186/s13045-021-01103-4
  52. Lv B, Wang Y, Ma D, et al. Immunotherapy: Reshape the tumor immune microenvironment. Front Immunol. 2022;13:844142. doi:10.3389/fimmu.2022.844142
  53. Wu J, Sun Y, Li J, et al. Analysis of prognostic alternative splicing reveals the landscape of immune microenvironment in thyroid cancer. Front Oncol. 2021;11:763886. doi:10.3389/fonc.2021.763886
  54. Menicali E, Guzzetti M, Morelli S, Moretti S, Puxeddu E. Immune landscape of thyroid cancers: New insights. Front Endocrinol (Lausanne). 2021;11:637826. doi:10.3389/fendo.2020.637826
  55. Lei X, Lei Y, Li JK, et al. Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett. 2020;470:126–133. doi:10.1016/j.canlet.2019.11.009
  56. Scheffel RS, Dora JM, Maia AL. BRAF mutations in thyroid cancer. Curr Opin Oncol. 2022;34(1):9–18. doi:10.1097/CCO.0000000000000797
  57. Wu P, Sun W, Zhang H. An immune-related prognostic signature for thyroid carcinoma to predict survival and response to immune checkpoint inhibitors. Cancer Immunol Immunother. 2022;71(3):747–759. doi:10.1007/s00262-021-03020-4
  58. Abou-Alfa GK, Mayer R, Venook AP, et al. Phase II multicenter, open-label study of oral ENMD-2076 for the treatment of patients with advanced fibrolamellar carcinoma. Oncologist. 2020;25(12):e1837–e1845. doi:10.1634/theoncologist.2020-0093
  59. Ping W, Hong S, Xun Y, Li C. Comprehensive bioinformatics analysis of toll-like receptors (TLRs) in pan-cancer. Biomed Res Int. 2022;2022(1):4436646. doi:10.1155/2022/4436646
  60. Li Z, Zhao J, Wu Y, et al. TRAF2 decrease promotes the TGF-β-mTORC1 signal in MAFLD-HCC through enhancing AXIN1-mediated Smad7 degradation. FASEB J. 2024;38(4):e23491. doi:10.1096/fj.202302307R
  61. Zhang Y, Zhang W, Wang Y, et al. Emerging nanotaxanes for cancer therapy. Biomaterials. 2021;272:120790. doi:10.1016/j.biomaterials.2021.120790
© Wroclaw Medical University