Single-cell and bulk transcriptomic analyses identify B-cell senescence-associated biomarkers in papillary thyroid carcinoma

Introduction

Thyroid cancer (THCA) is the most prevalent malignant neoplasm of the endocrine system worldwide, with a steadily increasing incidence (1). Environmental factors, such as radiation exposure, imbalanced iodine intake, and genetic mutations (e.g., BRAFV600E), are primary predisposing elements (2). Papillary thyroid carcinoma (PTC) constitutes over 80% of THCA cases, typically characterized by indolent growth and a favorable prognosis. However, some instances display heightened invasiveness, lymph node metastasis, and a substantial risk of recurrence (3). Clinically, this heterogeneity presents a significant challenge: while the majority of PTC patients have an excellent outcome with standard treatment, a subset experiences aggressive disease progression and recurrence. Accurately distinguishing between these indolent and aggressive forms at the time of diagnosis remains a major unmet clinical need, as over-treatment of low-risk patients and under-treatment of high-risk patients are both common concerns. Current PTC management primarily comprises surgical excision, radioactive iodine (RAI) therapy, targeted pharmacotherapy (e.g., tyrosine kinase inhibitors), and emerging immunotherapies (4). Nonetheless, these modalities exhibit limited efficacy in high-risk or refractory patients, with challenges of drug resistance and relapse underscoring the imperative for comprehending the tumor microenvironment (TME) and immune mechanisms to devise more precise therapeutic strategies.

B cell senescence denotes a condition of enduring cell-cycle cessation in B cells, triggered by stressors like oxidative damage or tumor-associated cues. This senescent state is linked with the emission of the senescence-associated secretory phenotype (SASP), encompassing cytokines and chemokines (5). Senescent B cells manifest compromised antigen-presenting and antibody-generating capacities, culminating in a compromised or skewed immune response. Within the tumor milieu, senescent B cells have the potential to foster an immunosuppressive milieu, bolster tumor evasion via heightened expression of immune checkpoints or restructuring of the TME, and conceivably amplify inflammation-driven tumor advancement (6). Recent findings propose a correlation between B cell senescence and unfavorable outcomes in diverse cancers (7). Comprehensive multi-omics investigations have unveiled the involvement of B cell senescence in configuring the immune microenvironment. In PTC, the infiltration of B cells is associated with either quiescence or advancement. Analysis of single-cell transcriptomes indicates that B cells infiltrating tumors can dictate the quiescent trajectory of PTC (8). Moreover, variations in the expression of senescence-related genes (SRGs) in PTC are linked to alterations in the immune cell composition, with age-related factors further impacting genomic variability and the immune milieu (9). Nonetheless, the precise mechanisms through which B cell senescence impacts PTC necessitate further elucidation, particularly its interaction with tumor heterogeneity and immune evasion.

Single-cell RNA sequencing (scRNA-seq) is an advanced high-throughput sequencing method designed to identify gene expression patterns within individual cells. This technology enables the exploration of cellular heterogeneity across various cell types, offering novel perspectives on tumor transcriptomics at the single-cell RNA level (8). By overcoming the limitations of conventional transcriptome sequencing, scRNA-seq provides valuable insights into the intricate landscape of different tumors and their microenvironments. Moreover, it facilitates the study of cell differentiation and the investigation of intercellular communications within biological systems (10).

This study integrated single-cell and bulk transcriptome sequencing datasets from public databases, including The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO), and conducted quantitative reverse transcription polymerase chain reaction (qRT-PCR) validation on clinical samples to identify biomarkers associated with B cell senescence in PTC. Through single-cell analysis and differential expression analysis, the biological pathways related to these biomarkers were elucidated, while the status of immune cell infiltration was assessed. The findings of this study aim to offer new insights for targeted therapy and clinical diagnosis in patients with papillary THCA, as well as informing future basic research on this disease. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2828/rc).

Methods

Data source

The TCGA-THCA dataset was downloaded from UCSC-Xena database (https://xenabrowser.net/datapages/) (accessed on April 10, 2025) as a bulk transcriptome dataset for PTC. The TCGA-THCA dataset included 505 PTC tumor tissue samples with complete progression-free interval (PFI) information and 59 adjacent non-tumor control samples. Clinical data and survival information for PTC samples were simultaneously retrieved from the UCSC-Xena database. Meanwhile, the TCGA-THCA dataset was divided into the TCGA-THCA-training set (354 samples) and TCGA-THCA-validation set (151 samples) according to a 7:3 split method. Furthermore, the GSE193581 dataset (GPL24676) was likewise downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) for the purpose of obtaining scRNA-seq data from PTC patients. Only 7 PTC tumor tissue samples and 4 PTC-adjacent non-tumor control samples were selected from this dataset; the other 2 control samples, which were adjacent non-tumor samples of differentiated thyroid carcinoma, were excluded.

A total of 866 SRGs (available online: https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-1.xls) were retrieved from the Aging Atlas database (11).

Single-cell analysis

The GSE193581 dataset was subjected to quality control using the Seurat package (v 5.0.1) (12). The selection criteria were as follows: 500< nFeature RNA (number of genes measured per cell) <5,500, nCount RNA (total gene count per cell) <24,000, and percentmt (percentage of mitochondrial genes per cell) <20%. To filter for highly variable genes (HVGs), the expression value of each gene in each cell was divided by the total expression value of all genes in that cell, then multiplied by a scaling factor of 10,000. The result was logarithmically transformed to obtain the normalized data. The FindVariableFeatures function was then used to identify the top 2,000 HVGs. The results were visualized using scRNAtoolVis package (v 0.0.7) (13) highlighting the top 5 HVGs. Subsequently, HVGs were subjected to principal component analysis (PCA) using the RunPCA function. The appropriate principal components (PCs) were selected using the ElbowPlot function from the Seurat package (v 5.0.1) (P<0.05). The cells were clustered using the FindNeighbors and FindClusters functions. Dimensionality reduction was then performed using uniform manifold approximation and projection (UMAP) with a resolution of 0.4. Subsequently, the marker genes provided in the literature (14) were applied to annotate different cell types. The result of cell annotation was visualized using the UMAP function. The expression levels of these marker genes were visualized using a bubble plot. Finally, histograms were used to display the abundance of each cell type across all samples in the PTC group and the control group.

Candidate genes acquisition

In the GSE193581 dataset, the FindMarkers function from the Seurat package (v 5.0.1) was used to select for differentially expressed genes (DEGs) in B cells between PTC samples and control group samples. The resulting genes were named “scRNA-seq-DEGs” with thresholds set at |average log₂ fold change (FC)| >0.25 and P<0.05. The top 10 scRNA-seq-DEGs were then ranked based on their |average log₂FC| values.

The DESeq2 package (v 1.38.0) (15) was utilized to identify bulk-DEGs between the PTC and control groups in the TCGA-THCA dataset (PTC vs. control), with thresholds set at |log₂FC| >0.5 and P<0.05. The bulk-DEGs were sorted by log₂FC, and the top 10 upregulated and downregulated genes were visualized using ggplot2 package (v 3.4.4) (16). The ComplexHeatmap package (v 2.14.0) (17) was subsequently used to create a heatmap illustrating expression profiles of the top 10 bulk-DEGs from each regulatory category.

Subsequently, the VennDiagram package (v 1.7.3) (18) was used to identify the intersections between scRNA-seq-DEGs-UP and bulk-DEGs-UP, as well as between scRNA-seq-DEGs-DOWN and bulk-DEGs-DOWN. The resulting genes were named candidate genes 1 and candidate genes 2, respectively. Candidate genes 3 were obtained by integrating candidate genes 1 and candidate genes 2. Furthermore, the VennDiagram package (v 1.7.3) was used to intersect the candidate genes 3 and SRGs. The obtained genes were designated as candidate genes. Enrichment analysis of candidate genes was conducted using the clusterProfiler package (v 4.7.1.003) (19), including Gene Ontology (GO) (P<0.05) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichments (P<0.05). The top 5 terms for each category and the top 10 KEGG pathways were presented according to the p-value. Additionally, a protein-protein interaction (PPI) network was established based on candidate genes using the STRING database (confidence score >0.4). The result was visualized using Cytoscape (v 3.8.1) (20).

Acquisition of biomarkers and construction of the risk model

In the TCGA-THCA-training set, univariate Cox regression analysis was conducted using the survival package (v 3.5-3) (21) [hazard ratio (HR) ≠1, P<0.2] (22) to identify candidate biomarkers. A forest plot was generated using the forestplot package (v 3.1.1) (23) to visualize the results. Then, candidate biomarkers were selected based on the proportional hazards (PH) assumption test (P>0.05). Subsequently, in the TCGA-THCA-training set, the randomForestSRC package (v 3.2.2) (24) was used to construct a random survival forest (RSF) model based on the biomarkers. Risk score for each patient was calculated using this model, with the number of trees (ntree) set to 17 and the maximum number of variables considered for splitting at each node (mtry) set to 1. The 354 PTC samples in the TCGA-THCA-training set were categorized into a high- (HRG) and low-risk group (LRG) based on the optimal cutoff value of risk score. The survminer package (v 0.4.9) (25) was used to plot Kaplan-Meier (K-M) curve (log-rank test, P<0.05) for the HRG and LRG to determine the difference in PFI between the 2 groups. In addition, a risk curve and survival status map were plotted to demonstrate the distribution of patients with PTC. The receiver operating characteristic (ROC) curve was used to evaluate the predictive performance of the model for 1-, 3-, and 5-year PFI using the survivalROC package (v 1.0.3.1) (26) and the area under the curve (AUC) value (AUC >0.6) was computed. A heatmap was used to display the expression profiles of biomarkers. Finally, to validate the risk model, the same analytical approach was used to assess the model in 151 PTC samples from TCGA-THCA-validation set.

Nomogram construction and assessment

In the TCGA-THCA-training set, risk score was combined with clinical characteristics (age, T stage, N stage, stage, and gender). The prognostic factors were identified via univariate Cox regression analysis (HR ≠1 and P<0.05). Subsequently, the PH assumption testing (P>0.05) and multivariate Cox regression analysis (HR ≠1 and P<0.05) were performed to screen for prognostic factors. A nomogram model was constructed based on prognostic factors using the rms package (v 6.5-0) (27) to predict the 1-, 3-, and 5-year survival probabilities of PTC patients. To assess the model’s predictive accuracy, The ROC curve was generated and AUC values were calculated using the survivalROC package (v 1.0.3.1) to evaluate the accuracy of the model in predicting PFI. The rms package (v 6.5-0) was used to plot the calibration curve. Additionally, a decision curve analysis (DCA) plot was generated using the rmda package (v 1.6) (28).

Analysis of the association between risk score and clinical characteristics

In the TCGA-THCA-training set, the association between risk score and clinical characteristics (including age, T stage, N stage, gender, and stage) was analyzed using the Wilcoxon test (P<0.05).

Immune microenvironment analysis

The single-sample gene set enrichment analysis (ssGSEA) was applied to assess the infiltration scores of 28 immune cells (29) in the TCGA-THCA-training set. The Wilcoxon test (P<0.05) was used to compare the immune cell scores between the HRG and LRG. Finally, Spearman correlation analysis via the psych package (v 2.2.9) (30) was used to evaluate associations between biomarkers and differentially expressed immune cells, providing insights into potential immune regulatory networks [|correlation coefficient (cor)| >0.30, P<0.05].

Tumor mutation burden (TMB) and drug sensitivity analysis

In the TCGA-THCA-training set, somatic mutation data were analyzed for mutation frequencies in HRG and LRG using the maftools package (v 2.16.0) (31). The top 20 genes with the highest mutation frequencies were displayed using a waterfall plot.

Additionally, to assess the drug sensitivity of HRG and LRG from the TCGA-THCA-training set, the half maximal inhibitory concentration (IC50) values for 138 conventional drugs, obtained from the GDSC database, were calculated using the pRRophetic package (v 1.2) (21). Spearman correlation analysis via the psych package (v 2.2.9) was used to evaluate associations between risk score and IC50 values (|cor| >0.30, P<0.05). The top 5 drugs with the strongest positive and negative correlations were separately displayed. Subsequently, the Wilcoxon test (P<0.05) was used to compare the differences in IC50 values of drugs significantly correlated with risk score between the HRG and LRG, and the top 10 drugs were presented.

Gene set enrichment analysis (GSEA) and expression analysis

The reference gene set was obtained from the “c2.cp.kegg.v7.4.symbols.gmt” file in the MSigDB. GSEA was performed using the clusterProfiler package (v 4.7.1.003) to identify pathways enriched for each biomarker. The correlation coefficients between each biomarker and all genes in the TCGA-THCA dataset were computed using the psych package (v 2.2.9). These correlation coefficients were ranked to generate a list of genes associated with each biomarker, with thresholds set at P<0.05. The top 5 signaling pathways were displayed according to ascending P values.

Additionally, in the TCGA-THCA dataset, the Wilcoxon test (P<0.05) was used to compare the expression differences of biomarkers between the PTC group and the control group.

Cell communication analysis

Based on all samples from GSE193581 dataset, cell-cell communication analysis among all cell types was performed using the CellChat package (v 1.6.1) (32). Potential ligand-receptor interactions were calculated to identify the signaling interactions between cell types.

Pseudo-time analysis and cell enrichment analysis

Since this study was related to B cells, B cells were defined as key cells. The dimensionality reduction and clustering methods used for key cells were the same as those in the “Single-cell analysis” section. The resolution was set to 0.1. These annotations were then visualized via UMAP. To explore the differentiation status and trajectory of key cell subclusters, pseudo-time analysis was carried out using the Monocle package (v 2.24.0) (33). Simultaneous analysis was conducted on the expression of biomarkers at each stage of cell differentiation.

To further investigate the association between our prognostic genes and canonical senescence markers, we analyzed the expression of CDKN1A (p21), CDKN2A (p16), and GLB1 in B-cell subclusters. Spearman correlation coefficients were calculated between each prognostic gene and each senescence marker across all B cells using the psych package (v 2.2.9) (P<0.05). For significant gene pairs, we identified co-expressing cells. The proportion of co-expressing cells was calculated relative to the total B-cell population. Co-expression patterns were visualized on the pseudotime trajectory, with co-expressing cells highlighted in red.

Additionally, the ReactomeGSA package (v 1.12.0) (34) was applied to the key cell subclusters in the GSE193581 dataset to identify enriched functional pathways. Using the lowest enrichment score as the benchmark, the enrichment score of each pathway was calculated, and the difference from this benchmark was determined. Pathways were then sorted in descending order of the difference, and the top 15 pathways with the largest differences were displayed.

Clinical samples

This study randomly obtained primary PTC tissues and matched noncancerous tissues from 40 patients who had previously undergone surgery at Fujian Medical University Union Hospital between August 2025 and October 2025. All tissues were snap-frozen in liquid nitrogen. None of these patients had received radiotherapy or chemotherapy before surgery, and they all signed the informed consent form prior to the operation. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Institutional Review Board (Ethics Committee) of Fujian Medical University Union Hospital (approval No. 2025KY823) and individual consent for this analysis was waived due to the retrospective nature.

qRT-PCR

The total RNA of samples was extracted using FastPure^® Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, Nanjing, China). The NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA) was used to quantify RNA. The cDNA was obtained by reverse transcription using a Vazyme HiScript III RT SuperMix for qPCR reagent kit. The qRT-PCR was performed using a 7500 real-time PCR system (Applied Biosystems, Waltham, USA). Relative mRNA expression levels of the indicated genes were normalized to β-actin as an internal control, with each sample analyzed in triplicate. The primer sequences were presented in https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-2.xls.

Statistical analysis

R (v 4.2.2) was used to conduct the bioinformatics analysis. Differences between the two groups were assessed using the Wilcoxon rank-sum test (P<0.05).

Results

Annotation of 8 cell types in the PTC

In the GSE193581 dataset, there were 29,832 cells and 23,537 genes before quality control. After quality control, the dataset contained 26,201 cells, with the number of genes remaining unchanged at 23,537 (Figure 1A,1B). Next, 2,000 HVGs were identified, and the top 10 HVGs were annotated, including TFF3 (Figure S1A). The PCA was performed, and the top 20 PCs were selected for further analysis (Figure S1B). The cells were divided into 18 clusters (Figure 1C). The clusters were then classified into 8 cell types: natural killer (NK) cells, tumor cells, T cells, B cells, myeloid cells, transitional follicular cells (TFCs), fibroblasts, and endothelial cells (Figure 1D). The expression levels of marker genes were presented using a bubble plot (Figure 1E). Additionally, in the control group, T cells accounted for the highest proportion at 68.68%, while in the PTC group, tumor cells accounted for the highest proportion at 33.32% (Figure 1F).

Figure 1 Annotation of cell types through scRNA-seq analysis. Distribution of nFeature_RNA, nCount_RNA, and percent.mt in samples before (A) and after (B) quality control. (C) Cells were divided into 18 cell clusters via the UMAP clustering method. (D) Eight cell types were identified after annotation. (E) Expression levels of marker genes in each annotated cell type. (F) Abundance of each cell type in the annotated samples. NK, natural killer; PTC, papillary thyroid carcinoma; TFC, transitional follicular cell; UMAP, uniform manifold approximation and projection.

Acquisition of 26 candidate genes

In the GSE193581 dataset, 3,267 scRNA-seq-DEGs were identified, with 935 genes upregulated and 2,332 genes downregulated in the PTC group (Figure 2A).

Figure 2 Acquisition of candidate genes. (A) Volcano plot of DEGs in the GSE193581 dataset. (B,C) Volcano plot and heatmap of DEGs in the TCGA-THCA dataset. (D) Candidate genes 1 was obtained by taking the intersection of scRNA-seq-DEGs-UP and bulk-DEGs-UP. (E) Candidate genes 2 was obtained by taking the intersection of scRNA-seq-DEGs-DOWN and bulk-DEGs-DOWN. (F) Candidate genes 3 were obtained by merging candidate genes 1 and 2, and 26 candidate genes were obtained by taking the intersection of candidate genes 3 and SRGs. (G,H) GO and KEGG enrichment analyses were performed on the candidate genes. (I) Construction of the PPI network. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; PPI, protein-protein interaction; PTC, papillary thyroid carcinoma; SRG, senescence-related gene; TCGA, The Cancer Genome Atlas; THCA, thyroid cancer.

In the TCGA-THCA dataset, 6,492 bulk-DEGs were identified, with 3,488 genes upregulated and 3,004 genes downregulated in the PTC group (Figure 2B,2C). By intersecting scRNA-seq-DEGs-UP and bulk-DEGs-UP, 106 candidate genes 1 were identified (Figure 2D). By intersecting scRNA-seq-DEGs-DOWN and bulk-DEGs-DOWN, 242 candidate genes 2 were identified (Figure 2E). A total of 348 candidate genes 3 were obtained by integrating candidate genes 1 and candidate genes 2. Furthermore, 26 candidate genes were obtained by integrating candidate genes 3 and SRGs (Figure 2F). Enrichment analysis demonstrated that the candidate genes were significantly associated with 548 GO terms and 74 KEGG pathways. For instance, GO terms included “neuron death” and “endoplasmic reticulum chaperone complex” (Figure 2G, available online: https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-3.xls). The KEGG pathways primarily included “apoptosis” and “human T-cell leukemia virus 1 infection” (Figure 2H, available online: https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-4.xls). In the PPI network, excluding the isolated proteins, the remaining 20 interacting proteins were identified, such as glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and EGR2 (Figure 2I).

Development of a risk model with high predictive accuracy

A total of 4 biomarkers (FOS, HSPA5, EGR2, and ABI3) were identified through univariate Cox regression analysis and the PH assumption test (Figure 3A, available online: https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-5.xls). Subsequently, the risk score for each patient was calculated. Based on the 4 biomarkers, a risk model was constructed. In the TCGA-THCA training set, the optimal cutoff value for risk score was −2.435901. Then, the 354 PTC samples in the TCGA-THCA-training set were classified into 145 HRG and 209 LRG samples. The risk curve and survival status map showed that the HRG experienced more deaths and had a higher mortality rate than the LRG (Figure 3B). The HRG had a shorter PFI with a statistically significant difference, indicating poorer progression-free survival outcomes (P<0.001) (Figure 3C). The ROC curve was used to evaluate the predictive performance of the model for 1-, 3-, and 5-year PFI. The AUC for 1-year PFI was 0.91, for 3-year PFI was 0.89, and for 5-year PFI was 0.84. These results indicated that the risk model exhibited moderate discriminative ability in predicting PFI at the 3 time points, with higher AUC values suggesting better predictive performance (Figure 3D). Expression level analysis showed that all 4 biomarkers were highly expressed in the LRG (Figure 3E).

Figure 3 Construction of the risk model. (A) Univariate Cox regression analysis and the PH assumption test. (B) Risk curve diagram and survival status distribution of the TCGA-THCA training set. (C) K-M curves drawn according to HRG and LRG in the TCGA-THCA training set. (D) ROC curve of the risk model in the TCGA-THCA training set. (E) Expression levels of prognostic genes in the HRG and LRG in the TCGA-THCA training set. AUC, area under the curve; CI, confidence interval; HRG, high-risk group; K-M, Kaplan-Meier; LRG, low-risk group; PFI, progression-free interval; PH, proportional hazards; TCGA, The Cancer Genome Atlas; THCA, thyroid cancer.

Additionally, with an optimal risk score cutoff value of 2.203105, the 151 PTC samples in the TCGA-THCA-validation set were divided into 66 HRG and 85 LRG samples, and analyzed in the same way, producing results consistent with those from the TCGA-THCA-training set (Figure 4A-4D). These results emphasized the reliability of the risk model’s predictions, showcasing its potential for guiding personalized treatment in PTC.

Figure 4 Validation of the risk model. (A) Risk curve and survival status distribution of the TCGA-THCA-validation set. (B) K-M curves drawn according to HRG and LRG in the TCGA-THCA-validation set. (C) ROC curve of the risk model in the TCGA-THCA-validation set. (D) Expression levels of prognostic genes in HRG and LRG in the TCGA-THCA-validation set. AUC, area under the curve; HRG, high-risk group; K-M, Kaplan-Meier; LRG, low-risk group; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; THCA, thyroid cancer.

Establishment of a good nomogram model based on prognostic factors

After selection via Cox regression analysis and PH assumption test, risk score and stage were selected as prognostic factors (Figure 5A,5B, available online: https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-6.xls). Based on these results, a nomogram model was developed. Each prognostic factor was assigned a single-item score (points), and the sum of these scores yielded the total score (total points). The total points were used to predict 1-, 3-, and 5-year survival probabilities, where higher scores were associated with higher mortality rates (Figure 5C). The ROC curve was used to evaluate the predictive performance of the nomogram model for 1-, 3-, and 5-year PFI. The AUC for 1-year PFI was 0.91, for 3-year PFI was 0.92, and for 5-year PFI was 0.89 (Figure 5D). The calibration curves for the 1-, 3-, and 5-year outcomes closely aligned with the reference curves, confirming the nomogram’s accuracy in predictions (Figure 5E). The DCA result indicated that the model exhibited favorable overall predictive performance (Figure 5F). These results suggested that prognostic factors could effectively predict the prognosis of PTC.

Figure 5 Independent prognostic analysis of clinical factors and risk scores. (A) Univariate Cox regression analysis and PH hypothesis test of clinical factors and risk scores. (B) Multivariate Cox regression analysis was used to identify independent prognostic factors. (C) Nomogram model for predicting the survival rate. (D) ROC curve of the nomogram model. (E) Calibration curve of the nomogram. (F) DCA curve of the nomogram model. **, P<0.01; ***, P<0.001. AUC, area under the curve; CI, confidence interval; DCA, decision curve analysis; N, node; PFI, progression-free interval; PH, proportional hazards; ROC, receiver operating characteristic; T, tumor.

Association between risk score and clinical characteristics

Risk scores showed significant differences across age, T stage, and stage. Patients in the age group over 60 had higher risk scores. Moreover, compared with the stage 1 group, risk scores were higher in the stage 2 and stage 3 groups (P<0.05) (Figure 6A). These results indicated that the risk scores were associated with the clinical characteristics.

Figure 6 Correlation analysis and immune infiltration analysis. (A) Correlation analysis between risk scores and clinical characteristics. (B) Analysis of immune cell infiltration in the HRG and LRG. (C) Correlation analysis between prognostic genes and differential immune cells. NS, no significance; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. HRG, high-risk group; LRG, low-risk group; N, node; T, tumor.

Significant association between risk score and immune microenvironment

The scores of 12 immune cell types, including monocytes, showed significant differences (P<0.05). The scores of most differentially expressed immune cell types were higher in the HRG, including central memory CD4 T cells and central memory CD8 T cells (Figure 6B). Spearman correlation analysis revealed that EGR2 showed the strongest positive correlation with eosinophils (cor =0.58, P<0.05), and FOS also showed the strongest positive correlation with eosinophils (cor =0.59, P<0.05), while HSPA5 showed the strongest negative correlation with immature dendritic cells (cor =−0.53, P<0.05) (Figure 6C). These results indicated that immune cells interacted with biomarkers and were involved in the progression of PTC.

High mutation of BRAF

In HRG and LRG, BRAF and NRAS exhibited the highest mutation frequencies, with missense mutations being the most common (Figure 7A,7B). BRAF hypermutation generally indicated enhanced tumor proliferative activity in PTC, which had implications for the treatment of PTC.

Figure 7 TMB, drug sensitivity and GSEA enrichment analysis. Gene mutation profiles of thyroid cancer samples in the HRG (A) and LRG (B). The waterfall plots show the top 20 genes with the highest mutation frequencies. (C) Correlation analysis between the risk score and the IC50 values of drugs. The drugs marked in the figure are the top 10 drugs with positive and negative correlations. (D) Differences in the IC50 values of the top 10 drugs between the HRG and LRG. (E-H) GSEA enrichment analysis of prognostic genes. The top 5 KEGG pathways associated with each gene are marked in the figure. ****, P<0.0001. GSEA, gene set enrichment analysis; HRG, high-risk group; IC50, the half maximal inhibitory concentration; KEGG, Kyoto Encyclopedia of Genes and Genomes; LRG, low-risk group; TMB, tumor mutation burden.

Additionally, the IC50 values of 17 drugs were statistically significantly correlated with risk scores, among which 12 showed a negative correlation (e.g., AZD6244) and 5 showed a positive correlation (e.g., sorafenib) (Figure 7C). Compared with the LRG, drugs such as LFM.A13 had lower IC50 values in the HRG, indicating that the HRG was sensitive to this drug and suggesting that patients in the HRG might have a better therapeutic response to this drug (Figure 7D).

The diversity of signaling pathways

GSEA revealed that FOS was enriched in 19 pathways, including “TGF-β signaling pathway”. HSPA5 was enriched in 30 pathways, including “peroxisome pathway”. EGR2 was enriched in 27 pathways, including “WNT signaling pathway”. ABI3 was enriched in 31 pathways, including “calcium signaling pathway” (Figure 7E-7H, available online: https://cdn.amegroups.cn/static/public/tcr-2025-1-2828-7.xls). These signaling pathways played a key role in the development of PTC.

Strong communication relationships between cell types

Compared with the control group, there was a communication from NK cells to B cells in PTC. Meanwhile, in the PTC group, tumor cells had more communication with other cell types (Figure 8A,8B). In addition, compared with the control group, the weights of communication from B cells to T cells and myeloid cells in the PTC group increased (Figure 8C,8D). In the control group, relatively high communication probabilities were observed between TFCs and myeloid cells, as well as between TFCs and B cells, with the ligand-receptor pair being MIF-(CD74 + CXCR4) in both cases, while in the PTC group, relatively high communication probabilities were observed between tumor cells and TFCs, with the ligand-receptor pair being MDK-SDC2 (Figure 8E,8F). The communication networks of these cell types provided insights for exploring the therapeutic mechanisms of PTC.

Figure 8 Cell-cell communication analysis. The number of interactions between key cells in the control group (A) and the PTC group (B). The intensity of interactions between key cells in the control group (C) and the PTC group (D). The nodes and colors in the figure represent different cell types, the size of the node circles indicates the number of cells of that type, and the thickness of the lines represents the number or intensity of communication. The bubble chart shows the ligand-receptor pairs between key cells in the control group (E) and the PTC group (F). The size of the bubbles in the figure represents the P value, and the color represents the average expression of the ligand-receptor pairs. NK, natural killer; PTC, papillary thyroid carcinoma; TFC, transitional follicular cell.

Analysis of B cell differentiation trajectory

A secondary dimensionality reduction and clustering analysis was performed on B cells. The top 15 PCs were selected for further analysis (Figure S2). B cells were divided into 3 clusters (Figure 9A). Pseudo-time analysis showed that the differentiation trajectories of B cells progressed from dark-colored to light-colored cells (Figure 9B). The cell development stages were divided into 5 periods, corresponding to different colors, with stage 1 representing the early differentiation phase (Figure 9C). Cluster 1 was predominantly found in the early stages (Figure 9D). Additionally, the expression level of ABI3 gradually increased from the early stage to the late and final stages. The expression level of EGR2 decreased from the early to middle stages, remained stable in the middle stage, and then increased sharply in the final stage. The expression level of FOS gradually increased slowly from the early to middle stages, decreased in the late stage, and rose again in the final stage. The expression level of HSPA5 gradually decreased from the early to middle stages, then increased again in the late and final stages. Meanwhile, the biomarkers were highly expressed in clusters 0 and 1 (Figure 9E,9F). The differentiation status of these cells was crucial in the development of disease.

Figure 9 Pseudo-time trajectory analysis and functional enrichment analysis of B cells. (A) B cells were divided into three cell subclusters via the UMAP clustering method. (B) Pseudo-time trajectory of B-cell differentiation. Dark blue indicates the early stage of differentiation, whereas light blue indicates the late stage of differentiation. (C) Different states of B-cell differentiation. (D) Differentiation trajectory of B-cell subclusters. (E,F) Expression trends of prognostic genes along the pseudo-time trajectory of B-cell differentiation. (G) Functional enrichment analysis of B-cell subclusters. The figure shows the expression heatmap of the top 15 differential pathways. (H) Heatmap showing Spearman correlations between prognostic genes (FOS, HSPA5, EGR2, ABI3) and canonical senescence markers (CDKN1A, CDKN2A, GLB1). (I-L) Co-expression patterns of significant gene pairs along the B-cell pseudotime trajectory: (I) FOS and CDKN1A, (J) HSPA5 and CDKN1A, (K) EGR2 and CDKN1A, (L) ABI3 and GLB1. Red circles indicate cells co-expressing both genes above their respective median expression levels. Percentages represent the proportion of co-expressing cells relative to total B cells. ns, no significance; **, P<0.01; ***, P<0.001. UMAP, uniform manifold approximation and projection.

Cell functional enrichment analysis showed that the “TWIK-related acid-sensitive K+ channel (TASK)” pathway was associated with key cell subclusters such as cluster 1 and cluster 2 (Figure 9G). These results helped to understand the biological functions of cells, thus providing references for the treatment of PTC.

To provide direct evidence linking our prognostic genes to functional senescence programs, we analyzed their correlation with canonical senescence markers (CDKN1A/p21, CDKN2A/p16, GLB1) in B cells. Spearman correlation analysis revealed four significant positive correlations among the 12 possible combinations: FOS with CDKN1A (r=0.20, P<0.001), HSPA5 with CDKN1A (r=0.11, P<0.001), EGR2 with CDKN1A (r=0.11, P<0.001), and ABI3 with GLB1 (r=0.07, P<0.01) (Figure 9H). The remaining eight combinations showed no significant correlation. We further quantified the co-expression patterns of these significant gene pairs along the B-cell pseudotime trajectory. Cells co-expressing both genes above their respective median expression levels were identified. The highest co-expression proportion was observed for FOS and CDKN1A (19.6% of total B cells), followed by HSPA5 and CDKN1A (16.9%), EGR2 and CDKN1A (3.1%), and ABI3 and GLB1 (3.0%) (Figure 9I-9L). Notably, co-expressing cells were distributed across multiple differentiation states, suggesting that senescence-associated programs may be activated at various stages of B-cell development within the tumor microenvironment. These findings provide transcriptional evidence that B cells expressing our identified biomarkers also co-express canonical senescence markers, particularly CDKN1A. The FOS-CDKN1A axis, with the highest correlation and co-expression proportion, may represent a key regulatory module in B-cell senescence within the PTC microenvironment.

Prognostic gene expression analysis

To elucidate the differential expression of the four prognostic genes in THCA, we conducted an expression analysis using the TCGA-THCA-training set. The results showed that ABI3 expression was significantly upregulated in THCA compared to normal samples, whereas FOS, HSPA5, and EGR2 expressions were downregulated in THCA relative to normal samples (P<0.001) (Figure 10A).

Figure 10 Expression of the prognostic gene. (A) Inter-group expression differences of prognostic genes in the TCGA-THCA-training set. (B) The expression differences of prognostic genes in clinical samples by qRT-PCR analysis. ***, P<0.001; ****, P<0.0001. PTC, papillary thyroid carcinoma; qRT-PCR, quantitative reverse transcription polymerase chain reaction; TCGA, The Cancer Genome Atlas; THCA, thyroid cancer.

Subsequently, we further verified their differential expression in clinical tissue samples through qRT-PCR experiments. The results were basically consistent with the bioinformatics analysis, confirming that the expression of ABI3 in PTC tissues was significantly higher than that in adjacent normal tissues, while FOS, HSPA5, and EGR2 showed a downward trend (P<0.001) (Figure 10B).

Discussion

PTC is the most common subtype of THCA, and its clinical outcomes are highly heterogeneous, ranging from indolent growth to invasive metastasis (35). Emerging evidence indicates that B-cell senescence significantly influences the tumor immune microenvironment (TIME), potentially facilitating immune evasion and hastening disease progression (36,37). This study integrates single-cell and bulk transcriptomic data to identify four biomarkers linked to B-cell senescence in PTC: FOS, HSPA5, EGR2, and ABI3. Using these biomarkers, we developed and validated a robust risk model and a nomogram model to effectively predict the PFI for PTC patients. Furthermore, we characterized differences in immune infiltration, drug sensitivity profiles, signal pathway enrichment, and B-cell differentiation trajectories. These findings offer novel insights into the interplay between B-cell senescence and PTC.

The study identified four distinct biomarkers with interconnected roles in tumorigenesis and development. The FOS gene, a pivotal member of the Fos family, encodes the FOS protein, a core component of the transcription factor AP-1 (38). This protein plays a crucial role in regulating cell proliferation, differentiation, and apoptosis. Studies have linked FOS to THCA stemness, with low expression correlating with poor prognosis (39). Meanwhile, studies have shown that the FOS/AP-1 signaling pathway plays a critical role in cellular senescence (40), suggesting a potential mechanistic link between this pathway and B cell aging that warrants further investigation. Heat shock protein family A (Hsp70) Member 5 (HSPA5), a member of the HSP70 family, acts as a calcium-binding protein in the endoplasmic reticulum (ER). It inhibits the activation of pro-apoptotic components under stress conditions and facilitates the transfer of folded proteins to maintain cellular protein synthesis, thereby preserving calcium homeostasis and ER stability in tumor cells (41). While HSPA5 overexpression is linked to cancer progression in various types of cancer (42), its down-regulation in PTC is associated with a poor prognosis, indicating disease-specific effects (43). Some studies have reported that HSPA5 is a key regulator in the unfolded protein response (UPR) (44). When ER stress occurs, it promotes cellular senescence by upregulating the UPR signaling pathway (45). Therefore, the downregulation of HSPA5 expression in PTC may disrupt the ER homeostasis of B cells and contribute to the acquisition of senescence-associated phenotypes. EGR2, a zinc finger transcription factor from the EGR gene family, modulates downstream gene expression by binding to specific DNA sequences (46). Reduced EGR2 levels in PTC tissues and cell lines have been observed, and its overexpression mitigates the carcinogenic impact of miR-224-5p in PTC (47). Furthermore, decreased EGR2 expression in PTC tumor tissue is linked to a poor prognosis. EGR2 overexpression inhibits IHH-4 and BCPAP cell proliferation, while EGR2 knockdown has the opposite effect (48). It is known that EGR2 plays a crucial role in the development and function of lymphocytes (49). Meanwhile, our research has observed that the dynamic expression pattern of EGR2 along the pseudotime trajectory of B cells suggests that it may be a key regulator of B-cell fate and might induce senescence in the PTC microenvironment. Targeting EGR2 may represent an effective novel therapeutic strategy for preventing tumor progression. ABI3, a novel SH3 domain molecule in the ABI protein family, is closely associated with tumor development. ABI3 expression correlates with immune cells and immunomodulatory genes in various cancers, impacting immune infiltration in tumor cells, patient prognosis, and serving as potential targets for immunomodulation (50). Studies have shown that FAT10 enhances cancer cell migration by stabilizing phosphorylated ABI3, potentially activating WRC (51). Although there are relatively few direct studies on ABI3 during the aging process, its involvement in cytoskeleton remodeling and immune function may affect the polarity of B cells and the interactions within the TME, thereby indirectly promoting the formation of the aging microenvironment. Subsequently, 40 clinical samples from PTC patients were collected from local hospitals, and the differential expression of these four biomarkers was confirmed using qRT-PCR, consistent with previous research.

The TIME comprises diverse cell types beyond cancer cells, such as macrophages, neutrophils, dendritic cells, myeloid-derived suppressor cells, and lymphocytes. Cancer cells can manipulate these cells to support cancer progression and evade the immune response (52). Our study conducted immune infiltration analysis, revealing significant discrepancies in the levels of 12 immune cell types between HRG and LRG. Particularly, the HRG displayed a more intricate immune cell infiltration pattern, potentially linked to immune dysfunction induced by B cell senescence (53). Elevated regulatory T (Treg) cells levels have been associated with poor prognosis in various cancers like breast cancer, ovarian cancer, lymphoma, and THCA aggressiveness (54-57). Studies have indicated that glioblastoma (GBM) patients lacking central memory CD4 T cells and NK cells exhibit improved overall survival (OS) (58). Conversely, oral squamous cell carcinoma (OSCC) patients with low central memory CD4 T cell infiltration had significantly shorter OS compared to those with abundant infiltration (P=0.0085) (59), suggesting central memory CD4 T cells as an independent prognostic biomarker for OSCC. Interestingly, eosinophils showed markedly higher infiltration in the LRG, aligning with previous research indicating their potential contribution to anti-tumor immunity and positive correlation with prognosis (60). Recent studies have reported reduced eosinophil levels in PTC tissues compared to controls (61). Eosinophils have been shown to exert an anti-tumor effect by modulating B cells to produce antibodies, albeit at a relatively low proportion (62). Spearman analysis revealed a negative correlation between the expression of most differential immune cells and prognostic genes. Collectively, these findings suggest that B cell SRGs could influence the PTC immune microenvironment by modulating specific immune cell subset functions.

Notably, the results of TMB analysis in this study showed that BRAF had the highest mutation rate among genes in both HRG and LRG. The BRAFV600E mutation is a pivotal driver in the pathogenesis and progression of PTC, facilitating tumor proliferation, invasion, and dedifferentiation through sustained activation of the MAPK pathway. Clinically, the BRAFV600E mutation serves as a molecular indicator for distinguishing between benign and malignant thyroid nodules, an independent prognostic factor, and a crucial target for precision therapy in advanced THCA (63). Combining ultrasound imaging characteristics with the BRAFV600E mutation has shown promise in predicting central lymph node metastasis risk preoperatively in papillary thyroid microcarcinoma patients (64), enhancing PTC prognosis and reducing recurrence likelihood. Regarding pharmacological sensitivity, LFM-A13 functions as a Bruton’s tyrosine kinase (BTK) inhibitor. Studies have demonstrated its ability to impede NF-κB-DNA binding, diminish cell proliferation, survival, and migration of cells, disrupt integrin-mediated adhesion to fibronectin, attenuate cellular responsiveness to tissue chemokines (CXCL12, CXCL13, and CCL19), and induce apoptosis (65). BTK is a key signaling molecule in B cells and is crucial for their development, activation, and survival (66). The predicted sensitivity of the HRG to the BTK inhibitor LFM-A13 in this study is particularly striking. We speculate that the B cell senescence signature identified in this study may mark a type of B cells that highly depend on the B cell receptor signaling pathway for survival in the tumor microenvironment. Therefore, patients with this signature may gain clinical benefits from BTK inhibitor treatment. While existing research has highlighted the anti-proliferative and pro-apoptotic properties of this inhibitor in cancers like breast and prostate cancer (67-69), its specific role in PTC remains ambiguous and necessitates further exploration. Future studies should prioritize experimental validation to determine the efficacy of LFM-A13, particularly in PTC models that can recapitulate the aging-associated tumor microenvironment enriched in B cells as observed in the HRG.

GSEA revealed the involvement of all four biomarkers in pivotal tumor-related signaling pathways, including TGF-β, Wnt, and calcium signaling pathways. The TGF-β pathway, known for its diverse biological functions, is intricately linked to the initiation and progression of malignant tumors. Activation of the TGF-β pathway has been associated with the downregulation of miR-124, leading to the promotion of epithelial-mesenchymal transition (EMT) and metastasis in non-small cell lung cancer (NSCLC) (70). The Wnt/β-catenin signaling pathway, a highly conserved pathway implicated in various cancers, has been extensively studied. Anaplastic THCA has been shown to exhibit elevated Wnt ligand expression, with tenascin-C amplifying Wnt signaling and tumor growth (71). Recent investigations have highlighted the role of KSRP in enhancing the nuclear localization and transcriptional activity of β-catenin by suppressing Wnt inhibitors, thereby fostering the progression and stemness of follicular THCA (72). Consequently, targeted interventions aimed at the Wnt/β-catenin pathway present a promising therapeutic avenue for managing THCA. Calcium ions serve as crucial regulators in cellular processes such as cell cycle regulation, apoptosis, and proliferation. Perturbations in the calcium signaling pathway can disrupt normal cellular functions, leading to uncontrolled cell proliferation and impaired programmed cell death, both of which are hallmark features of tumorigenesis (73). A Mendelian randomization study provided evidence of a causal relationship between elevated serum calcium levels and an increased susceptibility to pancreatic cancer (74).

At the single-cell level, our analysis of cell communication and pseudotemporal patterns revealed the pivotal involvement of B cells in PTC. Notably, we observed robust and extensive intercellular communication between B cells and myeloid cells within the PTC microenvironment, suggesting a potential impact on immune surveillance mechanisms. Moreover, our investigation unveiled the dynamic differentiation trajectory of B cells, showcasing stage-specific marker expression. Recent studies have reported a reduction in B cell numbers in PTC patients with lymph node metastasis. While excessive activation of B cells can trigger autoimmune disorders by fostering their proliferation, it concurrently plays a crucial role in bolstering anti-tumor immune responses (75). Furthermore, researchers have corroborated a substantial inverse relationship between B cell infiltration in PTC specimens and the occurrence of lymph node metastasis using immunohistochemistry (IHC) (76). Collectively, these findings underscore the potential influence of B cell senescence on immune evasion and lymph node metastasis in PTC.

To further confirm the association between our gene signature and functional B-cell senescence, we examined the co-expression of four prognostic genes with classic senescence markers along the B-cell differentiation trajectory. The study found that there was a significant correlation, especially between three genes (FOS, HSPA5, EGR2) and CDKN1A (p21, a key mediator of cell-cycle arrest). This not only provides transcriptional evidence but also suggests that B cells expressing these prognostic genes may indeed be in a senescent state. The FOS-CDKN1A axis showed the strongest correlation and the highest co-expression ratio (19.6%), indicating that FOS may play a particularly important role in regulating p21-mediated cell-cycle arrest in senescent B cells (77). Similarly, through the correlation between ABI3 and GLB1, our signature was linked to another senescence marker, namely the increased lysosomal content measured by β-galactosidase activity (78). These findings are consistent with the current understanding that senescence is not a single, homogeneous state but a complex program involving multiple interconnected pathways. The co-expression of our prognostic genes with unique senescence markers in different B-cell subsets indicates the existence of heterogeneity within senescent B cells, which may have different functional impacts on tumor immune regulation. However, we acknowledge that transcriptional co-expression, although supportive, does not constitute definitive evidence of functional senescence. Future studies should prospectively sort B cells from fresh PTC samples by flow cytometry and perform functional assays to conclusively determine the presence of functionally senescent B cells in the PTC microenvironment.

Nevertheless, this study is subject to several limitations. Firstly, while we conducted bioinformatics analysis and initial validation through qRT-PCR, our findings are primarily correlative. The inferred mechanistic role of FOS, HSPA5, EGR2, and ABI3 in B-cell senescence within the PTC microenvironment requires direct functional validation. Future studies employing in vitro models of B-cell senescence (e.g., using human primary B cells or cell lines treated with senescence-inducing agents) with gain- and loss-of-function experiments for these genes are essential to establish causality and elucidate the underlying molecular pathways. Secondly, the findings from drug sensitivity analysis, such as the potential of LFM-A13, are based on computational algorithms and should be interpreted with caution. These findings serve as a valuable resource for generating hypotheses but require prospective validation through in vitro and in vivo experiments before any clinical translation can be considered. Lastly, despite integrating single-cell and bulk transcriptome data, the heterogeneity among patients may obscure certain context-specific effects. Future research could incorporate other multi-omics analyses to delve deeper into the investigation.

Conclusions

In conclusion, this study conducted a systematic investigation into the transcriptional programs associated with of B-cell senescence in PTC, pinpointing four pivotal biomarkers: FOS, HSPA5, EGR2, and ABI3. Through the development of a risk model and integration of immune microenvironment analysis, a foundation was established for tailoring immunotherapeutic approaches for PTC.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2828/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2828/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2828/prf

Funding: The study was supported by the Clinical Research Center for Precision Management of Thyroid Cancer of Fujian Province (Grant No. 2022Y2006).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2828/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Institutional Review Board (Ethics Committee) of Fujian Medical University Union Hospital (approval No. 2025KY823) and individual consent for this analysis was waived due to the retrospective nature.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Pizzato M, Li M, Vignat J, et al. The epidemiological landscape of thyroid cancer worldwide: GLOBOCAN estimates for incidence and mortality rates in 2020. Lancet Diabetes Endocrinol 2022;10:264-72. [Crossref] [PubMed]
2. Banik GL, Shindo ML, Kraimer KL, et al. Prevalence and Risk Factors for Multifocality in Pediatric Thyroid Cancer. JAMA Otolaryngol Head Neck Surg 2021;147:1100-6. [Crossref] [PubMed]
3. Zhang Q, Li J, Shen H, et al. Screening and validation of lymph node metastasis risk-factor genes in papillary thyroid carcinoma. Front Endocrinol (Lausanne) 2022;13:991906. [Crossref] [PubMed]
4. Mohanty A, Afkhami M, Reyes A, et al. Exploring markers of immunoresponsiveness in papillary thyroid carcinoma and future treatment strategies. J Immunother Cancer 2024;12:e008505. [Crossref] [PubMed]
5. Laouris P, Muñoz-Espín D. Current Methodologies to Assess Cellular Senescence in Cancer. Methods Mol Biol 2025;2906:21-44. [Crossref] [PubMed]
6. Rodriguez JE, Naigeon M, Goldschmidt V, et al. Immunosenescence, inflammaging, and cancer immunotherapy efficacy. Expert Rev Anticancer Ther 2022;22:915-26. [Crossref] [PubMed]
7. Fettucciari K, Fruganti A, Stracci F, et al. Clostridioides difficile Toxin B Induced Senescence: A New Pathologic Player for Colorectal Cancer? Int J Mol Sci 2023;24:8155. [Crossref] [PubMed]
8. Li C, Wang P, Dong Z, et al. Single-cell transcriptomics analysis reveals that the tumor-infiltrating B cells determine the indolent fate of papillary thyroid carcinoma. J Exp Clin Cancer Res 2025;44:91. [Crossref] [PubMed]
9. Zhang Y, Chen Q, Niu L, et al. Multi-omics data analysis reveals the complex roles of age in differentiated thyroid cancer. Heliyon 2024;10:e33595. [Crossref] [PubMed]
10. Tan JK, Awuah WA, Roy S, et al. Exploring the advances of single-cell RNA sequencing in thyroid cancer: a narrative review. Med Oncol 2023;41:27. [Crossref] [PubMed]
11. Zhou X, Wu Y, Qin L, et al. Investigation of differentially expressed genes related to cellular senescence between high-risk and non-high-risk groups in neuroblastoma. Front Cell Dev Biol 2024;12:1421673. [Crossref] [PubMed]
12. Nagy C, Maitra M, Tanti A, et al. Single-nucleus transcriptomics of the prefrontal cortex in major depressive disorder implicates oligodendrocyte precursor cells and excitatory neurons. Nat Neurosci 2020;23:771-81. [Crossref] [PubMed]
13. Alani M, Altarturih H, Pars S, et al. A Roadmap for Selecting and Utilizing Optimal Features in scRNA Sequencing Data Analysis for Stem Cell Research: A Comprehensive Review. Int J Stem Cells 2024;17:347-62. [Crossref] [PubMed]
14. Lu L, Wang JR, Henderson YC, et al. Anaplastic transformation in thyroid cancer revealed by single-cell transcriptomics. J Clin Invest 2023;133:e169653. [Crossref] [PubMed]
15. Huang A, Sun Z, Hong H, et al. Novel hypoxia- and lactate metabolism-related molecular subtyping and prognostic signature for colorectal cancer. J Transl Med 2024;22:587. [Crossref] [PubMed]
16. Wang J, Wu N, Feng X, et al. PROS1 shapes the immune-suppressive tumor microenvironment and predicts poor prognosis in glioma. Front Immunol 2022;13:1052692. [Crossref] [PubMed]
17. Zhang X, Chao P, Zhang L, et al. Single-cell RNA and transcriptome sequencing profiles identify immune-associated key genes in the development of diabetic kidney disease. Front Immunol 2023;14:1030198. [Crossref] [PubMed]
18. Luo X, Xiang T, Huang H, et al. Identification of significant immune-related genes for epilepsy via bioinformatics analysis. Ann Transl Med 2021;9:1161. [Crossref] [PubMed]
19. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 2011;12:35. [Crossref] [PubMed]
20. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
21. Wang Y, Bin T, Tang J, et al. Construction of an acute myeloid leukemia prognostic model based on m6A-related efferocytosis-related genes. Front Immunol 2023;14:1268090. [Crossref] [PubMed]
22. Heng L, Bian H, Zhao C, et al. Identification of prognostic genes associated with mitochondria and macrophage polarization in prostate adenocarcinoma based on transcriptome and Mendelian randomization analysis. Discov Oncol 2025;17:5. [Crossref] [PubMed]
23. Liao X, Wang W, Yu B, et al. Thrombospondin-2 acts as a bridge between tumor extracellular matrix and immune infiltration in pancreatic and stomach adenocarcinomas: an integrative pan-cancer analysis. Cancer Cell Int 2022;22:213. [Crossref] [PubMed]
24. Zhao P, Zhen H, Zhao H, et al. Identification of hub genes and potential molecular mechanisms related to radiotherapy sensitivity in rectal cancer based on multiple datasets. J Transl Med 2023;21:176. [Crossref] [PubMed]
25. Shi Y, Wang Y, Dong H, et al. Crosstalk of ferroptosis regulators and tumor immunity in pancreatic adenocarcinoma: novel perspective to mRNA vaccines and personalized immunotherapy. Apoptosis 2023;28:1423-35. [Crossref] [PubMed]
26. Zhang S, Sun L, Cai D, et al. Development and Validation of PET/CT-Based Nomogram for Preoperative Prediction of Lymph Node Status in Esophageal Squamous Cell Carcinoma. Ann Surg Oncol 2023;30:7452-60. [Crossref] [PubMed]
27. Sachs MC. plotROC: A Tool for Plotting ROC Curves. J Stat Softw 2017;79:2. [Crossref] [PubMed]
28. Chen Y, He J, Pan X. Prediction of risk factors for preoperative deep vein thrombosis in patients with pelvic fracture. Front Surg 2025;12:1585460. [Crossref] [PubMed]
29. Shi J, Wu P, Sheng L, et al. Ferroptosis-related gene signature predicts the prognosis of papillary thyroid carcinoma. Cancer Cell Int 2021;21:669. [Crossref] [PubMed]
30. Zheng J, Zhang T, Guo W, et al. Integrative Analysis of Multi-Omics Identified the Prognostic Biomarkers in Acute Myelogenous Leukemia. Front Oncol 2020;10:591937. [Crossref] [PubMed]
31. Xu Q, Chen S, Hu Y, et al. Landscape of Immune Microenvironment Under Immune Cell Infiltration Pattern in Breast Cancer. Front Immunol 2021;12:711433. [Crossref] [PubMed]
32. Fang Z, Li J, Cao F, et al. Integration of scRNA-Seq and Bulk RNA-Seq Reveals Molecular Characterization of the Immune Microenvironment in Acute Pancreatitis. Biomolecules 2022;13:78. [Crossref] [PubMed]
33. Ma J, Wu Y, Ma L, et al. A blueprint for tumor-infiltrating B cells across human cancers. Science 2024;384:eadj4857. [Crossref] [PubMed]
34. Wan Z, Wang Y, Li A, et al. Single-cell transcription analysis reveals the tumor origin and heterogeneity of human bilateral renal clear cell carcinoma. Open Life Sci 2023;18:20220569. [Crossref] [PubMed]
35. Masone S, Velotti N, Savastano S, et al. Morbid Obesity and Thyroid Cancer Rate. A Review of Literature. J Clin Med 2021;10:1894. [Crossref] [PubMed]
36. Zheng H, Jiang W, Zhu S, et al. Prognostic significance of B cell senescence-associated genes as risk markers in prostate adenocarcinoma. Transl Cancer Res 2024;13:5771-83. [Crossref] [PubMed]
37. Zhou R, Zhou J, Muhuitijiang B, et al. Construction and experimental validation of a B cell senescence-related gene signature to evaluate prognosis and immunotherapeutic sensitivity in bladder cancer. Funct Integr Genomics 2022;23:3. [Crossref] [PubMed]
38. Cordier F, Creytens D. New kids on the block: FOS and FOSB gene. J Clin Pathol 2023;76:721-6. [Crossref] [PubMed]
39. Zhang T, Yan C, Ye Z, et al. The Identification of Three Key Genes Related to Stemness in Thyroid Carcinoma through Comprehensive Analysis. Comb Chem High Throughput Screen 2021;24:423-32. [Crossref] [PubMed]
40. Jing Y, Zheng W, Zhou Z, et al. Recent research advances of c-fos in regulating cell senescence. Arch Biochem Biophys 2025;769:110423. [Crossref] [PubMed]
41. Zhang C, Liu Q, Zhou Y, et al. HSPA5 Could Be a Prognostic Biomarker Correlated with Immune Infiltration in Breast Cancer. Dis Markers 2022;2022:7177192. [Crossref] [PubMed]
42. Li T, Fu J, Cheng J, et al. New progresses on cell surface protein HSPA5/BiP/GRP78 in cancers and COVID-19. Front Immunol 2023;14:1166680. [Crossref] [PubMed]
43. Dong W, Du D, Huang H. HSPA5 is a prognostic biomarker correlated with immune infiltrates in thyroid carcinoma. Endokrynol Pol 2022;73:680-9. [Crossref] [PubMed]
44. Ni B, Li S, Zhao L, et al. Novel glycoprotein SBSPON suppressed bladder cancer through the AKT signal pathway by inhibiting HSPA5 membrane translocation. Int J Biol Sci 2025;21:4586-603. [Crossref] [PubMed]
45. Bian Y, Wei J, Zhao C, et al. Natural Polyphenols Targeting Senescence: A Novel Prevention and Therapy Strategy for Cancer. Int J Mol Sci 2020;21:684. [Crossref] [PubMed]
46. Ashraf MU, Kang MH, Lim YT, et al. Egr2-dependent Mo-DCs regulate immunosuppressive NK cells, promoting neutrophil-mediated protective immunity against acute Listeria infection. Mol Cells 2025;48:100287. [Crossref] [PubMed]
47. Zang CS, Huang HT, Qiu J, et al. MiR-224-5p targets EGR2 to promote the development of papillary thyroid carcinoma. Eur Rev Med Pharmacol Sci 2020;24:4890-900. [PubMed]
48. Guo H, Zhang L. EGR1/2 Inhibits Papillary Thyroid Carcinoma Cell Growth by Suppressing the Expression of PTEN and BAX. Biochem Genet 2021;59:1544-57. [Crossref] [PubMed]
49. Gregory KJ, Morin SM, Schneider SS. Regulation of early growth response 2 expression by secreted frizzled related protein 1. BMC Cancer 2017;17:473. [Crossref] [PubMed]
50. Sun A, Cai F, Xiong Q, et al. Comprehensive pan-cancer investigation: unraveling the oncogenic, prognostic, and immunological significance of Abelson interactor family member 3 gene in human malignancies. Front Mol Biosci 2023;10:1277830. [Crossref] [PubMed]
51. Um H, Jeong H, Lee B, et al. FAT10 Induces cancer cell migration by stabilizing phosphorylated ABI3/NESH. Anim Cells Syst (Seoul) 2023;27:53-60. [Crossref] [PubMed]
52. Fang Z, Jiang C, Li S. The Potential Regulatory Roles of Circular RNAs in Tumor Immunology and Immunotherapy. Front Immunol 2020;11:617583. [Crossref] [PubMed]
53. Soma T, Nagata M. Immunosenescence, Inflammaging, and Lung Senescence in Asthma in the Elderly. Biomolecules 2022;12:1456. [Crossref] [PubMed]
54. Zhang H, Felthaus O, Eigenberger A, et al. Treg Cell Therapeutic Strategies for Breast Cancer: Holistic to Local Aspects. Cells 2024;13:1526. [Crossref] [PubMed]
55. Shao B, Sun B, Xiao Z. Transcriptome Analysis Unravels CD4(+) T-Cell and Treg-Cell Differentiation in Ovarian Cancer. Biomolecules 2025;15:1241. [Crossref] [PubMed]
56. Stirm K, Leary P, Wüst D, et al. Treg-selective IL-2 starvation synergizes with CD40 activation to sustain durable responses in lymphoma models. J Immunother Cancer 2023;11:e006263. [Crossref] [PubMed]
57. Fang X, Huang X, Lu J, et al. Causal role of immune cells in thyroid cancer: a bidirectional Mendelian randomization study. Front Immunol 2024;15:1425873. [Crossref] [PubMed]
58. Wu S, Yang W, Zhang H, et al. The Prognostic Landscape of Tumor-Infiltrating Immune Cells and Immune Checkpoints in Glioblastoma. Technol Cancer Res Treat 2019;18:1533033819869949. [Crossref] [PubMed]
59. Wu J, Zhang T, Xiong H, et al. Tumor-Infiltrating CD4(+) Central Memory T Cells Correlated with Favorable Prognosis in Oral Squamous Cell Carcinoma. J Inflamm Res 2022;15:141-52. [Crossref] [PubMed]
60. Zhai WY, Duan FF, Chen S, et al. A Novel Inflammatory-Related Gene Signature Based Model for Risk Stratification and Prognosis Prediction in Lung Adenocarcinoma. Front Genet 2021;12:798131. [Crossref] [PubMed]
61. Li M, Zhang J, Zhang Z, et al. Identification of Transcriptional Pattern Related to Immune Cell Infiltration With Gene Co-Expression Network in Papillary Thyroid Cancer. Front Endocrinol (Lausanne) 2022;13:721569. [Crossref] [PubMed]
62. Xie Z, Li X, He Y, et al. Immune Cell Confrontation in the Papillary Thyroid Carcinoma Microenvironment. Front Endocrinol (Lausanne) 2020;11:570604. [Crossref] [PubMed]
63. Scheffel RS, Dora JM, Maia AL. BRAF mutations in thyroid cancer. Curr Opin Oncol 2022;34:9-18. [Crossref] [PubMed]
64. Zhu D, Wu X, Zhang L, et al. Predictive Value of Ultrasound Imaging Characteristics and a BRAF V600E Nomogram for Central Lymph Node Metastasis Risk in Papillary Thyroid Microcarcinoma. Altern Ther Health Med 2023;29:139-43. [PubMed]
65. Aalipour A, Advani RH. Bruton's tyrosine kinase inhibitors and their clinical potential in the treatment of B-cell malignancies: focus on ibrutinib. Ther Adv Hematol 2014;5:121-33. [Crossref] [PubMed]
66. Katarzyna PB, Wiktor S, Ewa D, et al. Current treatment of systemic lupus erythematosus: a clinician's perspective. Rheumatol Int 2023;43:1395-407. [Crossref] [PubMed]
67. Guo W, Liu R, Bhardwaj G, et al. Targeting Btk/Etk of prostate cancer cells by a novel dual inhibitor. Cell Death Dis 2014;5:e1409. [Crossref] [PubMed]
68. Rozkiewicz D, Hermanowicz JM, Tankiewicz-Kwedlo A, et al. The intensification of anticancer activity of LFM-A13 by erythropoietin as a possible option for inhibition of breast cancer. J Enzyme Inhib Med Chem 2020;35:1697-711. [Crossref] [PubMed]
69. Sahin K, Tuzcu M, Yabas M, et al. LFM-A13, a potent inhibitor of polo-like kinase, inhibits breast carcinogenesis by suppressing proliferation activity and inducing apoptosis in breast tumors of mice. Invest New Drugs 2018;36:388-95. [Crossref] [PubMed]
70. Zu L, Xue Y, Wang J, et al. The feedback loop between miR-124 and TGF-β pathway plays a significant role in non-small cell lung cancer metastasis. Carcinogenesis 2016;37:333-43. [Crossref] [PubMed]
71. Hartmann HA, Loberg MA, Xu GJ, et al. Tenascin-C Potentiates Wnt Signaling in Thyroid Cancer. Endocrinology 2025;166:bqaf030. [Crossref] [PubMed]
72. Pan KF, Chou HL, Wang WL, et al. KSRP-mediated Wnt/β-catenin activation promotes follicular thyroid cancer progression and stemness. Br J Cancer 2025;133:1111-21. [Crossref] [PubMed]
73. Zheng S, Wang X, Zhao D, et al. Calcium homeostasis and cancer: insights from endoplasmic reticulum-centered organelle communications. Trends Cell Biol 2023;33:312-23. [Crossref] [PubMed]
74. Liu C, Wu H, Cai D, et al. Unraveling the causal relationship between serum minerals and pancreatic cancer: a Mendelian randomization study. Discov Oncol 2024;15:788. [Crossref] [PubMed]
75. Yu Y, Ouyang W, Huang Y, et al. Artificial intelligence-based multi-modal multi-tasks analysis reveals tumor molecular heterogeneity, predicts preoperative lymph node metastasis and prognosis in papillary thyroid carcinoma: a retrospective study. Int J Surg 2025;111:839-56. [Crossref] [PubMed]
76. Yang Z, Yin L, Zeng Y, et al. Diagnostic and prognostic value of tumor-infiltrating B cells in lymph node metastases of papillary thyroid carcinoma. Virchows Arch 2021;479:947-59. [Crossref] [PubMed]
77. Manousakis E, Miralles CM, Esquerda MG, et al. CDKN1A/p21 in Breast Cancer: Part of the Problem, or Part of the Solution? Int J Mol Sci 2023;24:17488. [Crossref] [PubMed]
78. Pedersen JJ, Duno M, Wibrand F, et al. β-Galactosidase deficiency in the GLB1 spectrum of lysosomal storage disease can present with severe muscle weakness and atrophy. JIMD Rep 2022;63:540-5. [Crossref] [PubMed]