Abstract

   Cancer stem cells (CSCs) play a crucial role in the progression of
   tumors, resistance to immunotherapy, and the development of thyroid
   carcinoma (THCA). However, the precise mechanisms involved remain
   poorly understood. In this study, 26 mRNAs linked to cancer stemness
   were identified through an analysis of 13 publicly available THCA
   transcriptomic datasets alongside a CRISPR dataset for thyroid cancer
   cell lines. Notably, we revealed an inverse relationship between the
   stemness of cancer cells and the efficacy of immunotherapy. Utilizing
   multiple machine learning approaches, we created and confirmed a model
   for tumor stemness cells (TSC), which included genes associated with
   cancer stemness. This model successfully forecasted patient outcomes
   and responses to treatment with Cytotoxic T-Lymphocyte Associated
   Protein 4 (CTLA4) immune checkpoint inhibitors (ICIs) across different
   types of cancer. Additionally, chromosomal amplifications in regions 1q
   and 8p were identified as intrinsic contributors to the onset of
   thyroid cancer. Specifically, based on single cell dataset, the
   amplification of CKS1B on chromosome 1q was revealed to advance the
   progression of thyroid cancer cells by boosting their proliferation.
   Furthermore, experimental results revealed that the heightened
   expression of CKS1B considerably enhanced the proliferation and
   invasion capabilities of THCA, suggesting it could serve as a possible
   therapeutic target for this condition. In summary, our study sheds
   light on the relationship between the traits of cancer stem cells and
   the efficacy of immunotherapy, providing a novel model for predicting
   the prognosis and individual responses to immunotherapy in patients
   with THCA.

Supplementary Information

   The online version contains supplementary material available at
   10.1007/s12672-025-02813-8.

   Keywords: Cancer stemness, Immunotherapy, Thyroid cancer, Single cell,
   CKS1B

Introduction

   Thyroid cancer, which arises from the cells of the thyroid gland, is a
   prevalent malignancy and the most frequently diagnosed malignant tumor
   within the endocrine system, with an annual incidence increasing by an
   average of 4.5% [[40]1]. Despite advancements in diagnostic techniques
   and therapeutic strategies, the ability to rapidly detect and precisely
   predict patient prognosis remains limited [[41]2]. Focusing on the
   inhibitory interactions that occur between immune checkpoint receptors
   and their corresponding ligands, including the programmed death 1
   (PD-1) pathway, has transformed the approach to treating advanced
   cancers through the use of immune checkpoint inhibitors. The
   integration of immunotherapy has transformed the approach to cancer
   care and enhanced patient survival rates. However, due to the genetic
   diversity of tumors, current biomarkers are not adequate for accurately
   categorizing patients for optimal treatment outcomes [[42]3]. Certain
   patients might not react positively and may encounter adverse effects,
   highlighting the significance of recognizing strong markers and
   enhancing combinations of biomarkers. Precisely identifying patients
   who will benefit from immunotherapy could expand the pool of eligible
   candidates for this treatment, potentially lessening the impact of
   chemotherapy side effects in individuals with thyroid cancer.

   CSCs exhibit a pronounced capacity for self-renewal, display
   considerable heterogeneity along with metastatic potential, and
   contribute significantly to tumor initiation, progression, and
   dissemination [[43]4, [44]5]. The diversity of cellular components and
   the immune-suppressive milieu within the tumor serve as distinct yet
   interrelated factors driving tumor progression and fostering treatment
   resistance [[45]6]. Prior investigations have demonstrated extensive
   negative associations between cancer stem cell traits and the immune
   system's response to cancer therapy. This phenomenon was noted even
   though tumors with high stemness exhibited elevated mutation rates,
   cancer-specific antigen expression, and intratumor heterogeneity across
   21 different solid cancer types [[46]7]. Recent study revealed a
   negative correlation between tumor stemness and the therapeutic
   outcomes of ICIs in various cancers. Their research introduced a robust
   predictive model to investigate how cancer stemness influences
   immunotherapy responsiveness [[47]8]. Despite these advancements, the
   specific impact of tumor stemness on ICI efficacy in thyroid cancer
   remains underexplored.

   In this research, we effectively employed the robust mRNA-based
   stemness index (mRNAsi) [[48]9] to pinpoint 26 mRNAs associated with
   cancer stemness. Through comprehensive analysis of transcriptomic
   profiles, single-cell data, and CRISPR-based cell line studies, our
   findings demonstrated an inverse association between cancer stemness
   and the efficacy of immunotherapy. Subsequently, we created and
   validated a scoring model for TSC, designed to forecast prognosis and
   responses to immunotherapy in individuals diagnosed with thyroid
   cancer. Additionally, we discovered that the overexpression of CKS1B, a
   crucial gene associated with cancer stemness, can stimulate the growth
   of thyroid cancer cells. This gene could potentially act as a valuable
   biomarker for predicting the cellular origins of malignant thyroid
   cells, supported by experiments conducted with cell lines.

   Our research offers fresh insights into the connection between CSCs and
   immunotherapy in the context of THCA. Furthermore, it introduces a
   well-performed predicting model for categorizing THCA patients based on
   their immune responses and survival outcomes.

Methods

Public data collection

   The dataset concerning thyroid cancer with survival data from The
   Cancer Genome Atlas (TCGA) was acquired from the UCSC Xena database
   [[49]10]. In addition, 12 GEO cohorts were studied to focus on thyroid
   cancer ([50]GSE104005, [51]GSE196264, [52]GSE29265, [53]GSE32662,
   [54]GSE35570, [55]GSE53157, [56]GSE5364, [57]GSE60542, [58]GSE65074,
   [59]GSE65144, [60]GSE76039, and [61]GSE82208) to examine mRNAs
   associated with cancer stemness. Furthermore, the gene expression
   profiles of a single-cell dataset for thyroid cancer with accession
   number [62]GSE184362 [[63]11] were collected from the GEO database to
   explore mRNAs linked to intrinsic stemness in cancer cells.

   To authenticate the precision of the model in forecasting the response
   to CTLA4 immunotherapy, we collected diverse datasets from cohorts
   treated with anti-CTLA4. These included the cohort by Riaz N [[64]12]
   ([65]GSE91061: melanoma patients administered with anti-CTLA4 and PD-1)
   and the cohort by Necchi [[66]13] (IMvigor210: individuals with
   advanced or metastatic urothelial carcinoma treated with Atezolizumab)
   obtained via the "IMvigor210CoreBiologies" R package. We amassed gene
   expression and clinical data from these cohorts.

   To explore the relationship between the stemness properties of cancer
   cells and their response to immunotherapy, a cohort of melanoma
   patients treated with ICIs was analyzed. This analysis incorporated
   single-cell RNA sequencing data and ICI response information,
   accessible from GEO under accession number [67]GSE115978 [[68]14].

Experimental data and material

   Human thyroid cancer cell line (KHM-5M) was obtained from the Cell Bank
   of the Committee for Conservation of Typical Cultures at the Chinese
   Academy of Sciences. These cell lines were cultured in Dulbecco's
   Modified Eagle Medium (DMEM) (HyClone, Logan, UT, USA), supplemented
   with 10% fetal bovine serum, 100 IU/mL penicillin, and 100 IU/mL
   streptomycin (Gibco, New York, USA).

   The siRNA targeting CSK1B was synthesized by Biotend Co., Ltd. Cells
   were transfected with 50 nM of siRNA for 24 h utilizing the jetPRIME
   transfection kit from Polyplus-transfection (Strasbourg, France).

   The Western blotting protocol was performed as follows: Cultured cells
   were washed with ice-cold PBS, and total protein lysates were extracted
   at 4 °C using RIPA lysis buffer (Beyotime, Shanghai, China)
   supplemented with 1% protease inhibitor cocktail (MedChemExpress, New
   Jersey, USA). After centrifugation at 12,000 g for 20 min at 4 °C, the
   supernatant was collected and mixed with loading buffer. The samples
   were then separated by 10% SDS-PAGE and transferred to a PVDF membrane.
   The membrane was blocked with 5% skim milk for 2 h at room temperature
   and incubated overnight at 4 °C with primary antibodies. After washing
   with Tris Buffered Saline, the membrane was incubated with secondary
   antibodies. Protein detection was performed using enhanced
   chemiluminescence reagents (Beyotime, Shanghai, China). The following
   antibodies were used for Western blot analysis: CKS1B (ab72639, Abcam,
   Waltham, Massachusetts, USA) and GAPDH (60004-1-Ig, Proteintech, Wuhan,
   China).

   To evaluate clonogenic potential, 1000 to 2000 cells were seeded into
   each well of a 6-well plate and cultured for 7 to 10 days. Colonies
   containing more than 50 cells were fixed with 0.2% crystal violet
   solution for 30 min. After washing three times with PBS, the colonies
   were photographed and quantified. Cells were evenly plated in a 6-well
   plate. Upon reaching confluence, a straight scratch was created using a
   10-μL pipette tip guided by a ruler. After PBS washing, the cells were
   maintained in serum-free medium. Incubation was carried out at 37 °C
   with 5% CO2, and images were captured at 0 and 24 h post-scratch. To
   assess cell migration, 3–5 × 10^4 cells were suspended in 200 μL of
   serum-free medium and added to the upper chamber of Transwell plates
   (BD Biosciences, Bedford, MA, USA). The lower chamber was filled with
   600 μL of medium containing 10% FBS. After 24 h of incubation at 37 °C,
   the cells were fixed with 4% paraformaldehyde for 30 min and stained
   with 0.25% crystal violet for an additional 30 min. Non-migrated cells
   on the upper side of the membrane were removed, and migrated cells on
   the lower side were imaged and counted.

Statistical analysis

   CRISPR screening of thyroid carcinoma cells was conducted using the
   DepMap portal ([69]https://depmap.org/portal/download). The CERES
   algorithm was employed to calculate dependency scores for approximately
   17,000 genes. Genes identified as essential exhibited a CERES score
   lower than −1 in 75% of the 21 THCA cell lines analyzed [[70]15].

   A novel pipeline was developed to establish a predictive model for
   cancer stem cells, as depicted in Fig. [71]1. Firstly, the calculation
   of mRNAsi was conducted on 12 GEO datasets as well as the TCGA-THCA
   dataset based on the pipeline described in previous research [[72]16].
   Utilizing a predefined reference gene signature containing specific
   stem cell-associated genes, each with assigned weight coefficients
   [[73]16], we extracted the common genes between this signature and our
   datasets. Stemness scores were computed by calculating the Spearman
   correlation between each sample's gene expression profile and the
   reference stem cell gene signature, which quantifies the similarity of
   gene expression patterns to the canonical stem cell gene set. To
   facilitate interpretability and comparative analysis, we linearly
   transformed the correlation scores to a standardized 0–1 range, where
   higher scores indicate greater stem cell-like characteristics.
   Subsequently, the connection between overall mRNA and mRNAsi was
   evaluated, and mRNAs showing noteworthy positive correlations in a
   minimum of seven out of thirteen sets (Cor > 0.2 and P < 0.05) were
   recognized as mRNAsi-related mRNAs, leading to the identification of
   135 mRNAsi-linked mRNAs.

Fig. 1.

   [74]Fig. 1
   [75]Open in a new tab

   Identification of key mRNAs linked to cancer stemness and their
   association with immunotherapy response. A Thyroid cancer cohorts used
   to identify mRNAs correlated with cancer stemness. B Schematic
   representation of the CRISPR screening process for mRNAs involved in
   thyroid cancer progression. C Venn diagram depicting the overlap
   between mRNAs associated with cancer stemness and those linked to
   thyroid cancer prognosis. D Box plot comparing the normalized
   expression levels of 26 stemness-related mRNAs in normal versus tumor
   thyroid tissues. E KEGG pathway enrichment analysis of mRNAs related to
   cancer stemness. F tSNE plot classifying malignant cells based on
   response phenotypes. Red points represent malignant cells from NR
   patients, and dark blue points represent malignant cells from TN
   patients. G tSNE plot visualizing AUCell scores for cancer stem
   cell-related gene sets in malignant cells, with red representing high
   stemness and blue indicating low stemness. H Box plot illustrating the
   distribution of AUCell scores across response phenotypes (NR vs. TN) in
   the SKCM cohort (Wilcoxon test; *P < 0.05; **P < 0.01; ***P < 0.001)

   As a result of the scarcity of clinical survival data on thyroid
   cancer, the TCGA-THCA dataset was split into two halves, Cohort1 and
   Cohort2, while all samples were merged into Cohort3. We utilized
   Cohort1 as the training set to develop machine learning model, then
   validated the model's accuracy using both Cohort2 and Cohort3.
   Afterwards, the effectiveness of the stemness model for cancer was
   created and assessed. To develop the cancer stemness model, various
   machine learning algorithms were utilized, such as random forest (RSF)
   [[76]17], elastic net (Enet) [[77]18], gradient boosting (GBM)
   [[78]19], survival-SVM [[79]20], ridge regression [[80]21], Stepcox
   [[81]22], plsRcox [[82]23], CoxBoost [[83]24], lasso [[84]25], and
   SuperPC [[85]26]. Specifically, through the input of the expression
   matrix of 26 selected cancer stemness related genes, the corresponding
   machine learning algorithm was used to predict the survival time and
   survival state of patients, and finally the score for each patient was
   calculated.

   To evaluate the effectiveness of PD-1/CTLA4 immunotherapy, we initially
   computed tumor immune dysfunction and exclusion (TIDE) scores using
   expression data from thyroid cancer patients. These expression profiles
   were subsequently analyzed on the TIDE database platform [[86]27]
   ([87]http://tide.dfci.harvard.edu/) to predict patient outcomes by
   uploading the expression matrix into [88]http://tide.dfci.harvard.edu/
   and ‘cancer type’ section selected ‘other’ and ‘Previous immunotherapy’
   section selected ‘no’. We prepared input files for SubMap analysis
   through a meticulous data processing workflow. Initially, we utilized
   raw expression data from the original dataset, applying rigorous
   filtering to remove low-expressed genes (defined as those with values
   less than 1 in over 90% of samples). To standardize the gene expression
   data, we identified the intersection of gene lists between different
   datasets and performed log2 transformation by adding 1 to each value
   (log2(x + 1)) to normalize the expression scale and mitigate potential
   numerical issues. Sample classification was conducted by categorizing
   samples into two distinct clusters: low-TSC and high-TSC. Each cluster
   was assigned a numeric rank to facilitate computational analysis, with
   TSC-low designated as rank 1 and TSC-high as rank 2. R function,
   generateInputFileForSubMap(), was used to systematically generate two
   critical files required for SubMap analysis: SubMap.gct and SubMap.cls.
   The SubMap.gct file encapsulates the log2-transformed gene expression
   matrix, preserving the molecular details of each sample, while the
   SubMap.cls file contains the corresponding sample classification
   information with their respective numeric ranks. This approach ensures
   that the files meet the stringent input requirements of the SubMap
   algorithm, enabling precise comparative analysis of gene expression
   profiles across different immunological clusters. Subsequently, the
   submap algorithm on the GenePattern portal was utilized to compare
   response probabilities between low- and high-TSC groups.

   In single cell analysis, the expression profiles for all samples were
   combined through the Merge() function in R (version 3.6.1). Next,
   library size normalization took place using the NormalizeData()
   function within Seurat [[89]28], resulting in normalized counts that
   were then log-transformed. The 2000 most highly variable genes were
   identified using the FindVariableFeatures() function. Principal
   components (PCs) were calculated from the expression profiles of these
   top 2000 variable genes, focusing on the top 30 components. For cell
   clustering, the FindNeighbors() and FindClusters() functions in Seurat
   were employed. When suitable, visualization was carried out using the
   RunUMAP() functions. A 2D Uniform Manifold Approximation and Projection
   (UMAP) algorithm facilitated the visualization of cells using the
   RunUMAP() and DimPlot() functions. Marker genes for each cluster were
   identified through the FindAllMarkers() function in Seurat. For any
   specific cluster, the FindAllMarkers() function revealed positive
   markers by comparing them with all other clusters.

   Clinical characteristics between high and low TSC score groups were
   compared using the Wilcoxon rank-sum test. KEGG enrichment is carried
   out through the enrichKEGG function in the clusterProfiler package
   [[90]29]. Additionally, GSEA enrichment of HALLMARK pathways was
   performed using the GSEA function from the clusterProfiler package. The
   HALLMARK dataset was obtained from the MisgDB database
   ([91]https://www.gsea-msigdb.org/gsea/index.jsp). The logFC values of
   input genes were calculated using the FindMarkers function based on the
   comparison between CKS1B+ malignant and CKS1B- malignant cells.
   Differences in immunotherapy response rates between the two groups were
   analyzed using the Chi-square test. Pearson’s correlation coefficient
   was calculated to examine the association between mRNA expression and
   mRNAsi. The relationship between TSC, CKS1B, and patient survival was
   evaluated using Kaplan-Meier survival analysis, with statistical
   significance assessed by the log-rank test. Time-dependent receiver
   operating characteristic (ROC) curves, generated using the 'pROC'
   [[92]30] R package, were employed to evaluate the prognostic and
   immunotherapy predictive value of TSC. The C-index was calculated based
   on the model-predicted scores, utilizing the coxph function to evaluate
   patients' survival time and survival status. Furthermore, to evaluate
   whether the TSC model could serve as an independent prognostic factor,
   we utilized the coxph function to perform multivariate Cox regression
   analysis. We calculated the impact of various clinical features on
   patient prognosis. Specifically, features with P < 0.05 were considered
   to have an independent predictive capability for patient survival,
   unaffected by the influence of other clinical features. Patients were
   stratified into groups based on the optimal cutoff value determined by
   the 'survminer' R package. A significance level of P or adjP < 0.05 was
   considered statistically significant.

Results

Cancer stemness is negatively associated with immunotherapy benefits

   A comprehensive analysis was conducted to investigate the association
   between cancer stemness and immunotherapy efficacy. Thirteen
   transcriptomic datasets for THCA were retrieved from the GEO and TCGA
   databases to identify mRNAs linked to cancer stemness, with mRNAsi
   values calculated for each patient (Fig. [93]1A). Using Pearson’s
   correlation coefficient, mRNAs exhibiting a significant positive
   correlation with mRNAsi in most cohorts (over 50%; 7 out of 13; Cor >
   0.2 and P < 0.05) were selected. A total of 135 mRNAs were then singled
   out as being related to cancer stemness, signifying the stemness level
   within the tumor cells. Furthermore, a thorough examination of
   loss-of-function screens based on CRISPR technology was conducted
   globally utilizing the DepMap database to pinpoint significant
   candidate genes linked to THCA malignancy. We finally identified 647
   genes critical for the survival of 21 THCA cell lines (CERES score <
   − 1 in 75% of malignant THCA cells). Among these, 26 mRNAs were
   selected based on the overlap between THCA mRNAsi-associated mRNAs and
   those identified through CRISPR screening (Fig. [94]1C).

   Consistent with above results, the majority of the identified genes
   displayed significant overexpression in tumor tissues compared to
   adjacent tissues, providing evidence of their role in tumor progression
   (Fig. [95]1D). Additionally, the KEGG enrichment demonstrated that
   these genes are mainly engaged in biological processes like the cell
   cycle and DNA repair, which further reinforces their connection to the
   stemness of tumor stem cells (Fig. [96]1E). All selected genes were
   prioritized for further investigation. To assess the impact of cancer
   stemness-related genes on immunotherapy efficacy, we analyzed a
   published single-cell dataset from melanoma patients treated with
   immune checkpoint inhibitors (ICIs), focusing on the association
   between cancer stemness and ICI response. After excluding samples
   lacking malignant cells, we concentrated on 23 patients, including 10
   non-responders (NR) and 13 treatment-naïve (TN) individuals. Using
   AUCell [[97]31], based on the 26 selected genes related to cancer
   stemness, we calculated the scores of cancer stemness-related genes in
   each cell using the AUCell_calcAUC method. We observed that malignant
   cells from the NR subgroup displayed significantly higher stemness
   levels (P < 0.001, Fig. [98]1H), indicating an inverse correlation
   between cancer stemness and ICI treatment outcomes.

   Taken together, our results indicated that these 26 mRNAs above
   selected can accurately characterized cancer stemness, and the stemness
   of cancer is negatively correlated with immunotherapy benefits.

Establishment of cancer stemness score based on THCA RNA-seq dataset

   To advance the development of a prediction model aimed at assessing the
   response to immunotherapy based on mRNAs associated with cancer
   stemness, a variety of machine learning algorithms were employed. These
   predictors were used to calculate scores for each sample across three
   cohorts. The evaluation of performance involved determining the average
   C-index for every algorithm, which showed that the majority of
   predictors yielded a high average C-index (Fig. [99]2A). Among the
   algorithms, the RSF approach demonstrated the greatest accuracy,
   achieving an average C-index of 0.993 (Fig. [100]2A), and was
   designated as the TSC score. Subsequently, patients were stratified
   according to the optimal cutoff determined by the 'survminer' R
   package, revealing that individuals with higher TSC scores had
   significantly reduced overall survival (OS) compared to those with
   lower TSC scores in cohorts 1–3 (Figs. [101]2B, S1A, B). Additionally,
   our results demonstrated that TSC is a reliable prognostic indicator,
   as evidenced by area under the curve (AUC) values. The AUC analysis
   highlighted strong predictive accuracy, with an average AUC exceeding
   0.98 across the three cohorts analyzed (Fig. [102]2C). To enhance the
   clinical applicability of our model, we investigated whether TSC could
   function independently of other clinical variables. A comprehensive
   multi-variable cox analysis indicated that TSC could act as an
   autonomous prognostic factor for assessing THCA prognosis (Fig.
   [103]2D). Recent progress in next-generation sequencing and large-scale
   data mining has enabled a comprehensive exploration and refinement of
   gene expression-based prognostic markers. To rigorously evaluate the
   performance of TSC compared to other markers, we systematically
   reviewed 79 models developed over the past decade (Table S1). Notably,
   TSC outperformed all other models in predicting survival outcomes
   across three independent THCA cohorts, achieving an average AUC greater
   than 0.98 in these cohorts (Fig. [104]2E).

Fig. 2.

   [105]Fig. 2
   [106]Open in a new tab

   Construction of a cancer stemness model using diverse machine learning
   approaches. A C-index for each algorithm across three independent
   cohorts, where the color of each grid represents the magnitude of the
   C-index obtained by different algorithms/algorithm combinations, and
   the bar graph on the right represents the average C-index across the
   three cohorts. B Kaplan-Meier survival curves comparing overall
   survival in TCGA-THCA patients stratified by high and low TSC scores in
   cohort3. C AUC values of the TSC model for predicting clinical status
   in three THCA cohorts. D Multivariate Cox regression analysis
   evaluating the independent prognostic significance of the TSC score in
   TCGA-THCA samples. E ROC values demonstrating the predictive
   performance of the TSC model compared to 79 other models for clinical
   status in three THCA cohorts

   In summary, the results demonstrated that the TSC model is a reliable
   prognostic model for THCA, exhibiting substantially greater predictive
   accuracy than other existing models.

Tumor stemness cell score as a potential predictor for CTLA4 immunotherapy
outcomes

   Given the inverse correlation between cancer stemness and immunotherapy
   efficacy, we investigated the potential of TSC as a predictive
   biomarker for immunotherapy outcomes in THCA patients. Using the submap
   algorithm on the GenePattern platform, we assessed the likelihood of
   immunotherapy response in patients stratified by high or low TSC
   levels. Notably, the low TSC group showed a significantly improved
   response to CTLA4 immunotherapy (P = 0.072, Bonferroni corrected P <
   0.01, Fig. [107]3A), underscoring the utility of the TSC model in
   predicting CTLA4 immunotherapy benefits.

Fig. 3.

   [108]Fig. 3
   [109]Open in a new tab

   Evaluation of the cancer stemness model as a predictor of immunotherapy
   efficacy. A Predicted response rates to immune checkpoint inhibitors
   (CTLA4 and PD1) across TSC subgroups (R: Response, NR: No Response). B,
   C Kaplan-Meier survival curves depicting overall survival in patients
   receiving CTLA4 immune checkpoint inhibitor (ICI) therapy in metastatic
   urothelial carcinoma (IMvigor210; B) and melanoma ([110]GSE91061; C). D
   Immunotherapy response rates stratified by TSC groups, with responses
   classified as complete (CR) or partial (PR), and non-responses as
   progressive (PD) or stable disease (SD). E, F Receiver operating
   characteristic (ROC) curves illustrating the TSC model's ability to
   predict response status in CTLA4 ICI-treated cohorts

   To further confirm the ability of TSC to predict the therapeutic
   effectiveness of CTLA4 ICIs, we gathered a variety of datasets
   associated with CTLA4 ICIs. Our findings consistently revealed that
   patients exhibiting lower TSC scores experienced notably improved
   overall survival rates after undergoing CTLA4 immunotherapy when
   compared to those presenting higher TSC scores (Fig. [111]3B and C).
   This indicates that elevated TSC scores might hinder the advantages of
   CTLA4 immunotherapy. Patients with elevated TSC scores demonstrated
   less favorable responses to the therapy, with merely 48% exhibiting a
   positive reaction (Fig. [112]3D). Significantly, our analysis
   demonstrated that TSC reliably predicts response to CTLA4-ICI therapy,
   as supported by high area under the curve (AUC) values. The AUC
   analysis revealed strong predictive performance, with values of 0.996
   at 1 year, 0.989 at 2 years, and 0.989 at 3 years in the IMvigor210
   cohort, and 0.900 at 1 year, 0.975 at 2 years, and 0.854 at 3 years in
   the SKCM (melanoma) cohort (Fig. [113]3E and F).

   To summarize, the TSC model serves as a dependable predictor for
   projecting the outlook of THCA patients and acts as an efficient
   biomarker for anticipating patient response to immunotherapy.

The amplification of CKS1B is associated with the proliferation of malignant
cells in THCA

   To investigate the clonal composition and cellular heterogeneity of
   thyroid cancer, we collected single-cell RNA dataset of thyroid
   carcinoma patients. Cells expressing fewer than 200 genes or exhibiting
   mitochondrial gene content exceeding 20% were filtered out. Following
   quality control, the retained cells were categorized into eight major
   cell types using established biomarkers. Figure [114]4A and B
   illustrate distinct clusters, encompassing B cells, endothelial cells,
   epithelial cells, fibroblasts, myeloid cells, plasma cells, T & NK
   cells, and dendritic cells.

Fig. 4.

   [115]Fig. 4
   [116]Open in a new tab

   Exploring core endogenous heterogeneity driven cancer stemness
   associated gene influencing prognosis of thyroid cancer. A UMAP plot
   depicting cell types annotated in thyroid cancer dataset with the
   accession number of [117]GSE184362. B Heatmap showing expression of
   each marker in each cell type. C Heatmap showing copy number variation
   of reference cells and epithelial cells. D Venn plot selected the
   cancer stemness associated mRNAs driven from thyroid malignant cell
   endogenous heterogeneity. E Box plot illustrating the normalized
   expression of CKS1B between normal and tumor tissues in thyroid cancer.
   F CKS1B expression in THCA tissues from human protein atlas. G Box
   plots comparing CKS1B expression between early (Stage I/II or T1/T2)
   and advanced (Stage III/IV or T3/T4) THCA groups. H GSEA enrichment of
   HALLMARK pathways between CKS1B- and CKS1B+ malignant cells

   Subsequently, we employed the inferCNV computational approach to
   examine both copy number variations (CNV) and clonal diversity within
   thyroid malignancy cells originating from epithelial cells (ECs). In
   the analysis, a total of 29,499 ECs from THCA tissues were examined,
   revealing that 21,052 of these cells displayed elevated CNV scores.
   This substantial number of cells exhibiting high CNV is indicative of
   malignancy, highlighting a significant characteristic of tumorigenesis
   in these instances (Fig. [118]4C). Moreover, particular amplifications
   were noted in chromosomal 8p and 1q, which emerged as significant
   driving variations associated with the pathogenesis of thyroid cancer.
   Notably, the amplification of the CKS1B gene located on chromosomal 1q
   was found to be linked with genes related to stemness (Fig. [119]4D).
   Previous research has established that the CKS1B gene has a pivotal
   role in promoting cell proliferation and invasion through the
   activation of the STAT3/AKT signaling pathway, particularly in cases of
   papillary thyroid carcinoma [[120]32]. This connection underscores the
   potential of CKS1B as a critical factor in the progression of thyroid
   malignancies.

   CKS1B expression was consistently found to be markedly increased in
   tumors when compared to the adjacent para-cancerous tissues in TCGA
   (Fig. [121]4E). Additionally, aligning with our bioinformatics
   findings, the protein levels of CKS1B were significantly higher in
   thyroid cancer tissues in contrast to normal thyroid samples (Fig.
   [122]4F). Moreover, our observations indicated that elevated CKS1B
   expression was particularly associated with the progression of tumor
   size and stage in TCGA (Fig. [123]4G).

   In conclusion, the upregulation of CKS1B may propel the advancement of
   malignant thyroid cells and could function as a valuable biomarker for
   predicting the cellular origins of THCA cancerous cells.

CKS1B enhances thyroid malignant cells’ proliferation, migration, and
invasion

   To investigate the impact of CKS1B on thyroid cancer cells, we
   categorized malignant tumor cells with an expression value above 1 as
   CKS1B+ malignant tumor cells, while those with an expression value of 0
   were classified as CKS1B- malignant tumor cells based on CKS1B
   expression levels. Our HALLMARK enrichment analysis revealed that CKS1B
   is capable of activating cell cycle-related pathways such as the G2M
   checkpoint and E2F targets (Fig. [124]4H). To further validate the
   oncogenic role of CKS1B in thyroid cancer, we generated CKS1B knockdown
   model in KHM-5M cell line. Their efficacy was confirmed through western
   blot analyses at the protein level (Fig. [125]5A). Notably, CKS1B
   knockdown led to a marked inhibition of clonogenic capacity in KHM-5M
   cells (Fig. [126]5B), suggesting a role for CKS1B in promoting KHM-5M
   cell growth potential. Furthermore, in wound healing and Transwell
   assays, CKS1B knockdown significantly impaired the migratory capacity
   of KHM-5M cells (Fig. [127]5C), indicating an important role for CKS1B
   in promoting the metastatic potential of thyroid cancer.

Fig. 5.

   [128]Fig. 5
   [129]Open in a new tab

   CKS1B enhances thyroid cancer proliferation and migration in vitro. A
   Validation of CKS1B knockout in KHM-5M cells via Western blot analysis.
   B Clonal formation ability of CKS1B knockout KHM-5M cells. C Evaluation
   of migratory abilities in CKS1B knockout KHM-5M cells using wound
   healing and Transwell assays. Scale bar, 100 μm

Discussion

   Immunotherapy, particularly immune checkpoint inhibitors, is now widely
   accepted as a second-line treatment for advanced cancers and is
   increasingly being approved for first-line use. Variability in tumors
   and their stem-like properties play crucial roles in how tumors evade
   the immune system and respond to immunotherapy. While significant
   research has been conducted on the relationship between cancer
   stem-like characteristics and anti-tumor immunity, direct evidence
   connecting tumor stemness to immune checkpoint inhibitor responses in
   thyroid cancer is presently lacking. Immunotherapy for thyroid cancer
   faces various challenges, with the main hurdle being the identification
   of suitable candidates for ICI treatment.

   In this research study, we employed an established stemness index
   [[130]9] to detect 28 cancer stemness-linked mRNAs in 13 thyroid cancer
   datasets. These gene sets were identified to participate in various
   important pathways including cell cycle, DNA replication, nucleotide
   excision repair, and the p53 signaling route. These pathways are
   crucial for regulating cancer proliferation [[131]33–[132]35].
   Furthermore, our investigation revealed an inverse relationship between
   stemness and the outcomes of immune checkpoint inhibitors, supported by
   data from an ICI-treated single-cell cohort of SKCM [[133]14].
   Moreover, several studies have indicated that patients with heightened
   cancer stemness experience deteriorated prognosis and heightened
   resistance to immunotherapy [[134]27, [135]36]. Subsequently, to assess
   the prognostic effect of tumor stemness genes in thyroid cancer
   patients, we utilized diverse machine learning algorithms to establish
   a prognostic model for survival and response to immunotherapy. The
   performance of this model was authenticated in three datasets, with the
   Random Survival Forest model being identified as the favored TSC model
   due to its improved reliability and precision. Analysis using the
   Kaplan-Meier method suggested that individuals with elevated TSC scores
   in THCA showed markedly poorer survival rates than those with lower TSC
   scores. Additionally, the AUC curve demonstrated an outstanding ability
   to predict outcomes, with the average AUC surpassing 0.98 across the
   three datasets. Additionally, our TSC score was validated as an
   autonomous predictor in multivariable Cox regression analysis. The TSC
   model showcased superior predictive performance for survival in
   contrast to conventional clinical and pathological characteristics,
   surpassing seventy-nine previously proposed models.

   The predictive value of cancer stemness for immunotherapy in various
   cancer types has been well-established, but its potential role in
   thyroid cancer remains largely unexplored [[136]8]. To assess the
   effectiveness of the cancer stemness index in forecasting the response
   to immunotherapy, we utilized the TIDE and Submap web server tools to
   analyze the immunotherapy response in thyroid cancer [[137]27]. The
   effectiveness of the TSC model in anticipating the reaction to CTLA4
   immune checkpoint inhibitors, particularly in cases of thyroid cancer,
   is emphasized by our findings. These findings suggest that the TSC
   model may have broad utility in predicting CTLA4 ICI treatment
   responses across different cancer types. Then, a thorough analysis
   showed that TSC exhibited a high level of accuracy in forecasting CTLA4
   ICIs responses across different datasets utilizing bulk RNA-Seq data,
   achieving an average AUC greater than 0.9.

   The presence of stemness is positively associated with increased
   intratumoral heterogeneity, potentially providing a protective
   mechanism for antigenic clones against elimination by the immune system
   [[138]7]. CNV in cancerous cells plays a crucial role in shaping tumor
   heterogeneity. Initial identification of inherent CNVs in malignant
   thyroid cells at the single-cell level indicated amplification of
   chromosomal 8p and 1q. Subsequent analysis linking mRNAsi and CNV
   profiles unveiled CKS1B as an endogenous mRNA linked to
   amplification-driven cancer stemness. As a protein-encoding gene, CKS1B
   associates with catalytic elements of cyclin-dependent kinases and
   enhances retinoblastoma cell proliferation, migration, invasion, and
   neovascularization via the MEK/ERK signaling cascade [[139]37]. Studies
   have earlier established CKS1B's utility as an immunotherapy response
   marker in pancreatic cancer patients [[140]38]. In line with these
   findings, our study confirms a notable upregulation of CKS1B in thyroid
   tumor tissue, correlating with tumor advancement. On a mechanistic
   level, overexpression of CKS1B triggers the G2M checkpoint and E2F
   target-related pathways, bolstering proliferation, invasion, and
   metastasis of thyroid cancer cells, presenting CKS1B as a potential
   therapeutic target for THCA.

   To summarize, based on single-cell dataset analysis from melanoma
   patients, our study suggests a potential inverse correlation between
   cancer stemness and immunotherapy response, which we applied to explore
   implications in thyroid cancer. Notably, we have developed a strong
   cancer stemness signature that can assess prognosis and potential
   benefits of immunotherapy, providing a useful tool for improving
   decision-making and monitoring strategies for individual patients with
   THCA. While we have shown the impressive predictive precision of the
   TSC in evaluating thyroid cancer prognosis and immunotherapy benefits,
   it is crucial to recognize some limitations in our research. The TSC
   model's ability to forecast immunotherapy outcomes for thyroid cancer
   is dependent on predictions from the submap and TIDE algorithm, and
   further validation of the TSC model using real thyroid cancer CTLA4
   immunotherapy cohorts is required. Moreover, conducting studies with
   larger cohort data is vital to validate the performance of the TSC
   model in thyroid cancer.

Supplementary Information

   [141]12672_2025_2813_MOESM1_ESM.tif^ (4.3MB, tif)

   Fig. S1 Validation of a cancer stemness model.A, B Kaplan-Meier
   survival curves comparing overall survival in TCGA-THCA patients
   stratified by high and low TSC scores in cohort1 (A) and cohort2 (B)

   Additional file1 (TIF 4417 kb)
   [142]Additional file2 (DOCX 92 kb)^ (92.6KB, docx)
   [143]Additional file3 (XLS 99 kb)^ (99KB, xls)

Acknowledgments