Introduction

Thyroid carcinoma (THCA), the prevalent endocrine malignancy globally, has witnessed a steady escalation in its incidence over recent years, manifesting as an annual growth rate of approximately 4%.^1–3 Despite the promising prognosis for most THCA patients following individualized treatment modalities such as surgery, radioactive iodine therapy, radiotherapy, systemic therapy, and personalized medical interventions, a subset of patients unfortunately succumbs to disease metastasis or recurrence.^4–6 Consequently, unraveling the intricate molecular regulatory mechanisms underlying THCA and pinpointing novel therapeutic targets represent crucial milestones in the pursuit of innovative treatment strategies for this malignancy.

The ubiquitin-proteasome system (UPS) orchestrates the degradation of misfolded and aggregated proteins as well as the turnover of cellular proteins, thereby maintaining the homeostatic balance and functionality of organisms.^7,8 This system encompasses ubiquitin (ub) activating enzymes (E1), ubiquitin conjugating enzymes (E2), ubiquitin ligases (E3), ubiquitin itself, and the 26S proteasome.⁹ It governs the degradation of approximately 80% of cellular proteins in eukaryotic cells and is integral to several fundamental cellular processes.¹⁰ Emerging research underscores the UPS’s involvement in a myriad of biological processes including the cell cycle, apoptosis, cell differentiation, autophagy, epigenetic regulation, angiogenesis, drug resistance, signal transduction, inflammation, and immune response.^11,12 Dysregulation of UPS components through functional mutations or abnormal expression plays a role in a range of diseases, such as neurodegenerative disorders and various immune-related conditions.⁹ Notably, growing evidence links the UPS to the onset and progression of malignancies, positioning it as a potential target for cancer treatment. For instance, one study showed that the expression of the E3 ubiquitin ligase MYCBP2 is correlated with inflammatory cell infiltration and the prognostic outcomes of immunotherapy in patients with thyroid carcinoma.^13,14 Since 2020, inhibitors of the UPS have been explored in clinical trials as potential anticancer agents.¹⁵ The combinatorial use of 26S proteasome inhibitors with other therapeutic agents has demonstrated encouraging outcomes in breast cancer treatment.¹⁶ Additionally, the proteasome inhibitor bortezomib has shown efficacy against anaplastic thyroid carcinoma. However, the functional roles and clinical implications of UPS-related genes in THCA remain to be comprehensively defined.¹⁷

In this study, we formulated a prognostic model for the central genes involved in the UPS. Following this, we stratified the THCA samples in the training set into two distinct subtypes based on UPS-PMS characteristics. Subsequent analyses identified KCNA1 as a pivotal prognostic gene within the UPS-PMS framework. Finally, we conducted comprehensive immunotherapy and drug sensitivity analyses to evaluate the therapeutic potential of targeting KCNA1.

Materials and Methods

Data Collection and Processing

We downloaded 512 thyroid carcinoma (THCA) tissues and 59 normal tissues from The Cancer Genome Atlas (TCGA) database¹⁸ (https://portal.gdc.cancer.gov), with file types in HTSEQ-FPKM. Among these, 507 thyroid carcinoma samples had clinical information. To ensure data completeness, we excluded samples lacking clinical data, resulting in 503 samples for clinical factor analysis. Additionally, the GSE29265 dataset was employed to validate the expression of genes within the prognostic model. All differential expression analyses were performed using the limma package.¹⁹ Based on previously published literature, we collected 797 ubiquitin proteasome system related genes, as shown in Supplemental Table 1.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA²⁰ is a robust bioinformatics framework that effectively aggregates highly correlated genes into modules and evaluates their associations with phenotypes. In this study, we utilized the R package (version 4.2.3) to analyze differentially expressed genes from TCGA-THCA. We determined power parameters using a soft threshold function, generated a cluster map through hierarchical clustering techniques, and constructed a module-trait association network. Modules were selected for further investigation based on Pearson’s correlation coefficients.

Functional Enrichment Analyses

Gene Ontology (GO)²¹ analysis is a commonly employed method for large-scale functional enrichment studies, encompassing biological processes, molecular functions, and cellular components. The Kyoto Encyclopedia of Genes and Genomes (KEGG)²² is a curated database resource that integrates diverse biological data across systematic, genomic, chemical, and health domains. Differentially expressed genes (DEGs) were annotated using GO and subjected to KEGG pathway enrichment analysis via the R package clusterProfiler. Statistical significance was set at a false discovery rate of < 0.05.

Gene expression data were analyzed using GSEA software (version 4.10) with the gene set database c2.cp.kegg.v7.4.symbols.gmt. Enrichment scores, p-values, and FDR-q values were computed. Criteria for inclusion in the enrichment evaluation required an FDR < 0.25 and p < 0.05.

Construction of UPS-Based Prognostic Model and Nomogram

We first identified DEGs across different clusters and selected intersecting DEGs for Lasso regression analysis. The glmnet R package facilitated this process, with lambda values adjusted to reduce the number of genes. A 10-fold cross-validation was used to build the model, with the optimal model selected based on lambda. The final risk score model was defined as Risk Score = Σ(exp_i * beta_i), where beta represents the multivariable coefficient, exp signifies expression levels, and i denotes the gene. A nomogram was constructed using the R package “rms”, incorporating clinical pathological features and core UPS-related genes.

Immune Characteristics Analyses

The single-sample Gene Set Enrichment Analysis (ssGSEA)²³ method was used to assess the enrichment scores for 19 immune cell subtypes. This approach estimates immune cell abundance in individual samples based on overall marker genes, representing the relative proportions of each subtype as enrichment scores. Differences in immune cell proportions between groups were evaluated using the Wilcoxon rank-sum test. Additionally, the ESTIMATE algorithm was employed to assess stromal and immune infiltration. Statistical significance was determined with p-values less than 0.05.

Identification of Different UPS-Related Subtypes of THCA

Cluster analysis was conducted utilizing the R package “ConsensusClusterPlus”, a robust tool for consensus-based clustering. Through this analysis, we identified two distinct subtypes of the data. The optimal number of clusters (k) was identified using the cumulative distribution function (CDF), which allowed for an objective selection of the most appropriate k value.

Constructing Potential ceRNA Networks for KCNA1

To identify miRNAs potentially interacting with KCNA1, we employed TargetScanHuman 8.0²⁴ (http://www.targetscan.org) and DIANA-microT²⁵ (http://diana.imis.athena-innovation.gr/DianaTools/index). Subsequently, miRNet 2.0²⁶ (https://www.mirnet.ca/miRNet/home) and starBase 3.0²⁷ (https://rnasysu.com/encori/) were used to predict and analyze the lncRNAs targeted by these miRNAs, facilitating the construction of a potential ceRNA network.

Drug Sensitivity Analyses

Connectivity Map (CMap)²⁸ provides an open resource for exploring the relationships between diseases and potential targets using transcriptomic data. By comparing DEGs between low and high expression groups of KCNA1, we submitted these DEGs to CMap to identify potential therapeutic targets and drugs for THCA. The heatmap displays the top 50 compounds with the most negative predicted scores. Additionally, the R package “oncoPredict” was used to estimate the half-maximal inhibitory concentration (IC50) values of anti-cancer drugs for each tumor sample in TCGA-THCA, with sensitivity assessed through ridge regression analysis.

Molecular Docking Analysis

HDOCK²⁹ (http://hdock.phys.hust.edu.cn/) is a widely utilized computing resource for molecular docking of DNA/RNA and protein-protein interactions. We utilized HDOCK for molecular docking modeling between KCNA1 and the drug gemcitabine. Subsequently, PyMOL (version 2.3.6)³⁰ was used to remove water molecules and dissociate the original ligand from the core target protein, with the results being visualized.

Western Blotting

Cell lysates were prepared with RIPA buffer (Sigma-Aldrich) containing protease and phosphatase inhibitors. Cells were harvested, lysed in the buffer for 30 minutes, and then centrifuged for 15 minutes. Proteins (30 µg) were mixed with loading buffer and denatured at 95°C for 5 minutes. SDS-PAGE was performed on a 10% polyacrylamide gel, followed by transfer to a nitrocellulose membrane (Bio-Rad) using a wet transfer system at 100V for 1 hour. The membrane was blocked with 5% non-fat milk in TBST (Tris-buffered saline with 0.1% Tween 20) at room temperature for 1 hour. Primary antibodies (rabbit anti-KCNA1/Kv1.1 monoclonal antibody (1:1000; ab32433, Abcam) and anti-Tubulin monoclonal antibody (1:1000 dilution; Proteintech, 11224-1-AP)) were diluted in 5% BSA in TBST and incubated overnight at 4°C. After TBST washes, the membrane was incubated with HRP-conjugated secondary antibodies at room temperature for 1 hour. Protein bands were visualized with the ECL detection system and imaged using the Chemidoc XRS+ system (Bio-Rad). Band intensities were quantified using ImageJ software (NIH).

Cell Counting Kit-8 (CCK-8) Assay

The human thyroid carcinoma cell lines (Bcpap) were procured from the Shanghai Institute of Cell Biology (Shanghai, China). Cells were cultured in complete DMEM (with 10% FBS and 1% penicillin-streptomycin) at 37°C in a 5% CO₂ incubator. For the CCK-8 assay, cells were seeded at 5×10³ cells per well in 96-well plates and allowed to adhere overnight. Following treatment with 100 μg of gemcitabine for 24, 48, or 72 hours, 10 µL of CCK-8 solution (Dojindo) was added to each well, and cells were incubated at 37°C for 2 hours. Absorbance was measured at 450 nm using a microplate reader (BioTek).

Half-Maximal Inhibitory Concentration (IC50)

To determine the IC50, the Bcpap cells were exposed to serially diluted gemcitabine. The viability data were analyzed using nonlinear regression to fit a dose-response curve, employing GraphPad Prism software. The IC50 value was derived from the curve as the concentration at which cell viability was reduced to 50% of the control.

Statistical Analyses

Bioinformatics, statistical analyses, and visualizations were carried out utilizing R software version 4.2.1. The Wilcoxon test was employed for comparing two groups, with statistical significance set at a P-value below 0.05. Pearson’s correlation coefficient was utilized for correlation assessments. Survival analysis was conducted with the survival and survminer packages in R, and Kaplan-Meier curves were generated for overall survival (OS). The Log rank test was used to assess differences in survival, with P-values <0.05 considered significant. Data were presented as mean ± standard deviation and analyzed using GraphPad Prism 8.0 software.

Results

Integrative Analysis Unveiled a Prognostic Model of Core Genes Associated with the UPS in THCA

To pinpoint critical genes in thyroid carcinoma (THCA), we analyzed differential expression data from TCGA-THCA, identifying 5854 genes with significant expression changes (p < 0.05) (Figure S1A). Among these, 3449 genes were upregulated and 2345 were downregulated. Furthermore, considering the vital role that Ubiquitin-proteasome system (UPS) plays in the development of THCA, we used a venn diagram and found 60 intersecting genes (26 upregulated and 34 downregulated) (Figure 1A and B). To further delineate the key genes, we conducted a weighted gene co-expression network analysis (WGCNA) using a soft threshold power of 7 to establish a scale-free network (Figure S1B). This analysis identified 16 gene co-expression modules through hierarchical clustering (Figure S1C). Following a comprehensive evaluation and screening process, four core modules were selected based on their significant association with THCA through gene expression patterns (Figure 1C). Subsequently, we employed a Venn diagram and successfully identified 28 DEGs by analyzing the intersection of the four key gene modules with the UPS (Figure 1D), which were entered into the LASSO-Cox regression model for further screening (Figure 1E and F). As a result, we successfully identified six candidate genes: potassium voltage-gated channel modifier subfamily S member 3 (KCNS3), deltex E3 ubiquitin ligase 4 [Homo sapiens (DTX4), zinc finger and BTB domain containing 21 (ZBTB21), potassium channel tetramerization domain containing 16 (KCTD16), potassium voltage-gated channel subfamily A member 1 (KCNA1), and F-box protein 15 (FBXO15), showing some degree of correlation with each other (Figure 1G). Concurrently, we determined the penalty coefficients for key genes using LASSO analysis and established a risk index by applying these coefficients to gene expression levels. The resulting formula for the ubiquitin-proteasome system prognostic model score (UPS-PMS) is as follows: UPS-PMS = (0.073 × KCNS3 exp) + (−0.008 × DTX4 exp) + (−0.221 × KCTD16 exp) + (0.463 × ZBTB21 exp) + (1.992 × KCNA1 exp) + (0.554 × FBXO15 exp). To explore the relationship between the six core genes and clinicopathologic features, including age, gender, T, M, N and pathological stage, we employed a heatmap analysis (Figure S1D). Our analysis revealed a notable similarity in the expression profiles of KCNS3 and DTX4, which exhibited higher expression levels in samples with more advanced clinical features, particularly in the T, M, and N staging parameters. In contrast, genes KCTD16, ZBTB21, KCNA1, and FBXO15 demonstrated a greater degree of expression similarity, with higher expression levels observed in samples with early clinical features according to T, M, and N staging. We identified six core genes associated with UPS in thyroid carcinoma and developed a prognostic model that predicts patient prognosis with high accuracy.

Figure 1 Development of a prognostic model for UPS core genes in thyroid cancer. (A) Cross-over genes were found using Venn diagram. (B) A volcano plot identified 26 up-regulated and 34 down-regulated genes. (C and D) Weighted gene co-expression network analysis (WGCNA) was employed to pinpoint core modules and identify pivotal genes. (E and F) The prognostic model was formulated via LASSO-Cox regression. (G) Correlation analysis elucidated the interrelationships among the genes within the model (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

The Development of a Novel Prognostic Assessment Tool UPS-PMS for THCA

Initially, samples were categorized into low and high UPS-PMS groups. To assess differences between these groups, we visualized patient survival times and statuses alongside the expression levels of the six core genes. This analysis revealed a notably higher incidence of mortality in the high UPS-PMS group. Additionally, gene expression analyses revealed a statistically significant upregulation of ZBTB21, KCTD16, KCNA1, and FBXO15 in the high UPS-PMS group compared to the low UPS-PMS group. Conversely, the expression levels of KCNS3 and DTX4 were significantly reduced in the high UPS-PMS group (Figure 2A). Subsequent survival analysis demonstrated that patients in the high UPS-PMS group had a significantly poorer prognosis compared to those in the low UPS-PMS group (Figure 2B). Receiver operating characteristic (ROC) curve analysis was employed to assess the sensitivity and predictive specificity of the risk scores. The area under the curve (AUC) values, specifically 0.764, 0.913, and 0.831 (Figure 2C), established UPS-PMS as a notable prognostic indicator for patients with THCA. Following that, we conducted the univariate and multivariate Cox regression analyses which demonstrated that UPS-PMS was an independent prognostic factor for THCA (Figure 2D and E). Furthermore, among the numerous clinicopathologic characteristics, gender and stage exhibited significant differences between the UPS-PMS high and low groups as shown in Table 1. The results were integrated into the nomogram, which also included other clinicopathological characteristics, namely age, gender, tumor stage, and T, M, and N classification. The integration of these variables facilitates accurate risk stratification and prognostic assessment of THCA, in line with personalized medicine by integrating comprehensive clinical and biomarker panels (Figure 2F). In the calibration plot, the nomogram-predicted result is in close alignment with the reference line, which suggests that the prediction nomogram is an optimal predictive model (Figure 2G). In conclusion, we successfully developed and validated a novel ubiquitin-proteasome system-based prognostic assessment tool for thyroid carcinoma, which provides accurate risk assessment and prognostic prediction for personalized medicine.

Table 1 Difference in Clinicopathological Characteristics Between High and Low UPS-PMS Groups

Figure 2 Development of a thyroid cancer prognosis assessment tool. (A) Survival time and status of patients in the high-UPS-PMS and low-UPS-PMS groups, and visualization of the expression of six core genes. (B) Survival analysis and (C) ROC curve analysis were performed between high-UPS-PMS group and low-UPS-PMS group. (D and E) Univariate and multivariate Cox regression analysis of UPS-PMS. (F) Prediction nomogram and (G) calibration curve of the model.

Association Analysis of UPS-PMS with Immune and Cancer Progression Pathways in THCA Patients

To elucidate the biological functions of DEGs, we utilized volcano plot to visualize gene expression variations between the high and low UPS-PMS groups. This analysis identified 2371 genes with significant upregulation and 3482 genes with downregulation (Figure 3A). Top 100 most upregulated and downregulated genes were shown in Figure 3B. To elucidate the pathways enriched in THCA patients with high versus low UPS-PMS, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses with pathway annotation shown in Supplementary Tables 2 and 3. Downregulated genes were predominantly associated with immune-related pathways, including antigen processing and presentation, lymphocyte-mediated immunity, MHC protein complex, T cell receptor binding, and Th1/Th2 cell differentiation. In contrast, upregulated genes were primarily linked to ubiquitination and cancer progression pathways, such as ubiquitin-dependent ERAD, positive regulation of protein ubiquitination, ubiquitin ligase binding, vascular endothelial growth factor receptor binding, and TGF-beta signaling (Figure 3C and D). Additionally, Gene Set Enrichment Analysis (GSEA) demonstrated that UPS-PMS was negatively correlated with the innate immune system, infectious disease response, IL-18 signaling, translation, and natural killer cell-mediated cytotoxicity and so on (Figure 3E and F). Taken together, these pathways involve immune responses, cell-cell interactions, and protein degradation and are central to cancer progression and immune regulation.

Figure 3 Pathway enrichment analysis of differential genes in UPS-PMS. (A) A volcano plot visualized the differential expression patterns between the high and low UPS-PMS groups. (B) The expression profiles of the top 100 most significantly differentially expressed genes were depicted. (C) Gene ontology (GO), (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) and (E and F) gene set enrichment analysis (GSEA) were harnessed to uncover the primary enriched pathways associated with the differential genes.

Associations Between Immune Cell Heterogeneity and UPS-PMS in THCA

To investigate the relationship between risk groups and immune cell profiles, we utilized the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm to assess immune cell enrichment within the tumor microenvironment. The heatmap (Figure 4A) illustrated the distribution of various immune cells across samples and their association with UPS-PMS. The correlation analysis revealed a negative correlation between UPS-PMS and T cells (r = −0.516, p = 1.43e-08), helper T cells (r = −0.245, p < 0.001), and Th2 cells (r = −0.583, p < 0.001), while showing a positive correlation with NK cells (r = 0.530, p < 0.001), DC cells (r = 0.367, p < 0.001), Th1 cells (r = 0.430, p < 0.001), and macrophages (r = 0.436, p < 0.001) (Figure 4B). Additionally, ESTIMATE algorithm analysis demonstrated that the Immune Score and ESTIMATE Score were negatively correlated with UPS-PMS, with the low UPS-PMS group exhibiting higher scores. The Stromal Score, however, did not show a significant correlation with UPS-PMS (Figure 4C). Next, we examined immune cell infiltration in the samples, revealing distinct levels of immune infiltration between high and low UPS-PMS groups (Figure 4D). Furthermore, it was observed that immune cells in the high UPS-PMS group exhibited a greater degree of interconnectivity compared to the low-UPS-PMS group. This phenomenon may be indicative of enhanced synergism between immune cells in the context of high UPS-PMS tumor development (Figure 5A). Moreover, the expression of six immune checkpoints was consistently elevated in the low UPS-PMS group compared to the high UPS-PMS group (p < 0.05). Subsequent studies revealed a negative correlation between immune checkpoint expression and UPS-PMS (Figure 5B–G). By examining the immune microenvironment, the poor prognosis of patients with high UPS-PMS thyroid carcinoma may be associated with immunosuppression. This suggests the potential for new immunotherapy strategies.

Figure 4 The heterogeneity of immune cells and their correlations with UPS-PMS. (A) Immune infiltration levels of 18 immune cells in high-UPS-PMS group vs low-UPS-PMS group. The higher or lower the Immune infiltration, the darker the color (red is high infiltration level and blue is low infiltration level), (B) Lollipop plot illustrating the correlation of UPS-PMS with various immune cells. Stem length denotes correlation magnitude, with greater green saturation indicating stronger associations. Circle size at stem termini represents P-value significance, where smaller circles denote higher statistical relevance, (C) The scatterplot depicts the relationship between UPS-PMS and ImmuneScore, StromalScore, and ESTIMATEScore, with a straight line fitting the overall correlation. The box plots visualize the differences in the above three scores between the high and low UPS-PMS groups, (D) Box plots delineate immune infiltration disparities among immune cell types between high- and low-UPS-PMS groups (*p < 0.05, ***p < 0.001).

Figure 5 The level of immune cell connectivity and checkpoint expression suggests the potential for immunotherapy. (A) Heatmap of Mantel correlations between UPS-PM core gene expression and the infiltration scores of 18 immune cells in high- and low-UPS-PMS group, (B-G) Scatter plots reveal the associations between the expression of six immune checkpoints and UPS-PMS levels, while box plots highlight the expression disparities of these checkpoints across high and low UPS-PMS groups (*p < 0.05, **p < 0.01, ***p < 0.001).

Identification and Clinical Implications of UPS-PM-Related Subtypes in THCA

Using consensus clustering based on the expression levels of six key UPS-PM genes in THCA samples, we classified the 503 samples into two distinct subtypes (cluster_1 and cluster_2) with the optimal K-value of 2 (Figure 6A–C). Cluster_1, comprising 252 patients, exhibited higher UPS-PMS levels compared to cluster_2, which included 251 patients (Figure 6D). Survival analysis indicated that patients in cluster_2 had significantly better outcomes than those in cluster_1 (Figure 6E). By combining subtypes with UPS-PMS co-analysis, we discovered that cluster_1, which exhibited a high UPS-PMS, demonstrated the poorest survival outcome. Conversely, cluster_2, which displayed a low UPS-PMS, exhibited the most favorable survival outcome (Figure 6F). To further explore the differences between the two isoforms, we investigated clinical pathologies and found characteristic N, M, stage showed evident difference (Figure 6G–I) while gender and T had no significant variations (Table 2). Our findings indicate that cluster_1 exhibited significantly elevated UPS-PMS and were associated with a higher prevalence of terminal N and M stages and overall staging. In contrast, cluster_2 demonstrated lower UPS-PMS and a more favorable prognosis, exhibiting earlier clinicopathologic features compared to cluster_1 (Figure 6J).

Table 2 Difference in Clinicopathological Characteristics Between Cluster_1 and Cluster_2

Figure 6 The recognition of two different UPS-PM-related subtypes of THCA. (A)The cumulative distribution curve of each k value from 2 to 10 in the consensus cluster analysis, (B) The area under the curve of each k value, (C) The heat map of the two different UPS-PM-related subtypes of THCA, (D) The UPS-PMS value of cluster_1 and cluster_2, (E) The survival probability analysis of cluster_1 and cluster_2, (F) Survival analysis of THCA patients with the combined subtype and UPS-PMS, (G–I) Bar scale diagram showing differences in distinct subtypes between (G) distinct stage, (H) metastasis and (I) nodes metastasis, (J) The heatmap demonstrated the different clinicopathological features of each sample, survival, UPS-PMS and subtypes (***p < 0.001).

Identifying Biological Difference in UPS-PM-Related THCA Subtypes

To further elucidate the biological differences between the UPS-PM-related subtypes, we conducted a differential expression analysis comparing cluster_2 to cluster_1. This analysis revealed 2680 DEGs from the TCGA-THCA dataset, with logFC > 2 and p < 0.05 (Figure S2A). Moreover, the expression profiles of the top 100 UPS-PMS-related DEGs between the two THCA subtypes were displayed in a heatmap, highlighting the unique genetic heterogeneity of each subtype (Figure S2B). To elucidate the pathway enrichment of up- and down-regulated genes, we performed GO and KEGG enrichment analyses with pathway annotation shown in Supplementary Tables 4 and 5. The upregulated genes were predominantly associated with terms related to immune processes, including antigen processing and presentation, leukocyte cell-cell adhesion, cytokine production regulation, IκB kinase/NF-κB signaling, and adaptive immune response. These terms are linked to immune cell activation, proliferation, and differentiation. Conversely, the downregulated genes were primarily involved in pathways such as oxidative stress response, ubiquitin-dependent ERAD, cellular oxidative stress response, intrinsic apoptotic signaling, and oxidative phosphorylation—pathways associated with disease onset and progression (Figure 7). In summary, our genetic analysis reveals that the distinct biological characteristics between the two UPS-PMS-related subtypes are closely associated with immune processes.

Figure 7 GO & KEGG analyses of the up- and down-regulated genes. (A and B) The number of genes enriched in each KEGG pathway term is shown as the circle size, and the categories are shown by different colors.

Cellular Infiltration and Immune Checkpoint Differences in THCA Subtypes

Given the evident link between pathway enrichment and immune functions, we employed the ssGSEA algorithm to assess immune cell infiltration in THCA samples. This analysis provided deeper insights into the immune profiles and differences between the two subtypes (Figure 8A). A notable discrepancy was observed in the ESTIMATE algorithm analysis, with the Cluster_2 group exhibiting a higher score than the Cluster_1 in the IMMUNE score category. In contrast, Cluster_1 exhibited a significantly higher Stromal score compared to Cluster_2 (Figure 8B). We examined the immune landscape of Cluster_1 and Cluster_2 by analyzing the correlation between the expression levels of six key ubiquitin-proteasome system pathway molecules (UPSPMs) and immune cell infiltration. Our results reveal a stronger correlation between the expression profiles of six UPSPMs and immune cell types in Cluster_1 compared to Cluster_2. This enhanced correlation in Cluster_1 indicates a more complex interplay between the immune landscape and the ubiquitin-proteasome system pathway, underscoring its significance in the disease’s immunopathogenesis (Figure 8C). We employed the ssGSEA method to quantify the enrichment scores of various immune cell subgroups and explore their correlations with the two clusters. The analysis revealed that NK cells, macrophages, and Th1 cells were more abundant in Cluster_2, while T cells, helper T cells, and Th2 cells were more prevalent in Cluster_1 (Figure 8E). Additionally, given the pivotal role of immune checkpoints in immunotherapy, we compared their expression between the clusters. Cluster_2 exhibited higher levels of six key immune checkpoints, including TIGIT and PDCD1LG2, relative to Cluster_1 (Figure 8D). In summary, our analysis using ssGSEA and ESTIMATE algorithms uncovered notable differences in immune cell infiltration and immune checkpoint expression between two thyroid carcinoma subtypes. These findings underscore the importance of the immune microenvironment in relation to UPS pathway molecules and their role in the immunopathology of these subtypes.

Figure 8 Cellular infiltration and immune checkpoint differences in THCA subtypes. (A) Heatmap illustrating the immune infiltration levels of 19 distinct immune cell types across thyroid cancer samples, (B) Box plots delineate differences in ImmuneScore, StromalScore, and ESTIMATEScore between cluster_1 and cluster_2, (C) Heatmap of Mantel correlations between UPS-PM core gene expression and the infiltration scores of 18 immune cells in cluster_1 and cluster_2, (D) Comparative box plot analysis of immune checkpoint expression across clusters cluster_1 and cluster_2, (E) Box plots delineate immune infiltration disparities among 19 immune cell types between cluster_1 and cluster_2 (*p < 0.05, **p < 0.01, ***p < 0.001).

Identifying a Core Prognostic Marker in THCA

To find the core genes that best characterize the model’s prognosis, the following analyses were conducted. Initially, the GSE29265 dataset was introduced for the purpose of analyzing the expression of the six genes that constitute the model in cancerous and non-cancerous tissues. The results showed that the expression of KCNA1, ZBTB21, KCTD16, FBXO15 was observed to be elevated in cancerous tissues in comparison to their expression levels in normal tissues. However, the opposite was true for KCNS3 and DTX4 (Figure 9A). In order to compare the influence of each candidate gene on survival in UPSPM prognostic signature, the overall survival was further evaluated. The results demonstrated that only two candidate genes, KCNA1 and ZBTB21, exhibited a poor prognosis when expressed at high levels (Figures 9B and S3E). To assess the prognostic value of the candidate genes, we performed a time-dependent ROC curve analysis for one-year, three-year, and five-year overall survival, calculating the corresponding AUC values. The analysis revealed that the AUCs for KCNA1 were 0.788, 0.794, and 0.642 at one, three, and five years, respectively. In contrast, ZBTB21 exhibited AUCs of 0.356, 0.851, and 0.704 for the same time points (Figures 9C and S3J). Additionally, the survival status of other genes is either as protective factors in the development of thyroid cancer or not significantly meaningful (Figure S3A–D, F-I). Considering the better integrated prediction performance of KCNA1, KCNA1 was identified as the leading gene in the prognostic characterization of UPSPM.

Figure 9 Identification of KCNA1 as the Leading Gene in the UPS-PM Signature. (A) The expression level of KCNS3, KCNA1, DTX4, ZBTB21, FBXO15 and KCTD16 in the GSE29265 dataset, (B) The survival analysis of KCNA1, (C) The area under the curve of KCNA1, (D) The miRNA intersection between TargetScan and DIANA-microT database, (E-I) The correlation analysis between KCNA1 expression and (E) hsa-miR-92a-3p, (F) hsa-miR-92b-3p, (G) hsa-miR-32-5p, (H) hsa-miR-125a-5p and (I) hsa-miR-125b-5p expression, (J) The site prediction between KCNA1 and hsa-miR-92a-3p and hsa-miR-92b-3p and hsa-miR-125b-5p, (K) The lncRNA intersection between miRNet and STARBASE database, (L) The construction of ceRNA regulatory networks (*p < 0.05, **p < 0.01, ***p < 0.001).

To elucidate the molecular mechanisms and regulatory factors driving the upregulation of KCNA1 expression, we constructed a competing endogenous RNA (ceRNA) regulatory network for KCNA1 in THCA. Utilizing predictions from TargetScan and the DIANA-microT database, we identified a total of 17 miRNAs potentially involved in this regulatory network (Figure 9D). Scatter plots demonstrated the correlation between KCNA1 expression and the target miRNAs (Figure 9E–I). Notably, KCNA1 expression exhibited significant negative associations with several miRNAs: hsa-miR-92a-3p (r = −0.203, p < 0.001), hsa-miR-92b-3p (r = −0.233, p < 0.001), hsa-miR-32-5p (r = −0.151, p < 0.001), hsa-miR-125a-5p (r = −0.255, p < 0.001), and hsa-miR-125b-5p (r = −0.316, p < 0.001). Subsequently, we utilized the TargetScan database to identify potential binding sites for the target miRNAs within KCNA1 (Figure 9J). To further investigate, we employed the miRNet and starBase databases, which predicted 13 long non-coding RNAs (lncRNAs) potentially interacting with these miRNAs (Figure 9K). To identify lncRNAs negatively correlated with the miRNAs, we further screened and validated these lncRNAs in THCA using the starBase database. Based on these findings, we constructed four ceRNA regulatory networks implicated in THCA (Figure 9L). In conclusion, in order to deeply investigate the complex relationship between the genetic elements affecting the expression of UPSPM, we chose KCNA1 gene as the main research object. And a three-layer ceRNA regulatory network centered on KCNA1 was constructed.

Drug Sensitivity Analysis and Immunotherapy Potential of KCNA1 Expression in THCA

In light of the growing significance of immunotherapy in the management of the disease, we undertook an investigation into the therapeutic potential of immunotherapy. The findings showed that the group with higher KCNA1 gene expression had higher tumor immune dysfunction and exclusion (TIDE) scores and lower immunotherapy response rates (Figure 10A and B). To explore potential therapeutic interventions targeting KCNA1, we conducted a Connectivity Map (CMap) analysis. We developed a KCNA1-associated gene signature, which included 168 upregulated and 162 downregulated genes. Heatmap clustering of this signature highlighted the enrichment of the top 50 compounds (Figure 10D). Besides, we also conducted another drug sensitivity exploration with the Cancer Therapeutics Response Portal (CTRP), obtaining 218 potential drugs (Figure 10C). Subsequently, the intersection of the two resulted in the identification of five drugs (gemcitabine, TW-37, PAC-1, LY-2183240 and isoliquiritigenin) (Figure 10E). A review of the literature indicates that gemcitabine is the only drug extensively used in clinical practice for THCA. Besides, our findings revealed that the half maximum inhibitory concentration (IC50) for gemcitabine was significantly elevated in the cohort exhibiting high KCNA1 gene expression (Figure 10F). To determine the effects of gemcitabine on KCNA1, we performed a comprehensive set of experiments. Firstly, the CCK8 assay demonstrated that gemcitabine is capable of inhibiting the viability of cancer cells (Figure 10G). Subsequently, the WB assay revealed that gemcitabine can suppress the protein expression of KCNA1 (Figure 10H). The IC50 assay demonstrated that KCNA1 knockdown increased cellular sensitivity to gemcitabine (Figure 10I). To investigate the three-dimensional structure of KCNA1, we utilized AlphaFold for tertiary structure prediction and consulted the PDB database. Subsequently, we employed the H-DOCK server to predict potential binding sites between KCNA1 and gemcitabine, as illustrated in Figure 10J.

Figure 10 Drug Sensitivity analysis and Immunotherapy Potential of KCNA1 expression in THCA. (A) The TIDE score between low KCNA1 expression group and high KCNA1 expression group, (B) The response ratio between low KCNA1 expression group and high KCNA1 expression group, (C) Correlation analysis of KCNA1 expression and chemotherapeutic agents in the CTPR database, (D) Reflection of DEGs between KCNA1 high and low expression groups in the Cmap database to drugs in different cell lines, (E) The chemotherapeutic agents intersection between CTPR and Cmap database, (F) CCK8 assay for inhibition of Bcpap cell proliferation in administered gemcitabine and control groups, (G) Protein expression of KCNA1 after treatment with gemcitabine and control groups, (H) The IC50 changes of gemcitabine after KCNA1 knockdown, (I) The IC50 profile of gemcitabine in the KCNA1 highland expression group, (J) The docking site prediction between KCNA1 and gemcitabine (**p < 0.01, ***p < 0.001, ****p < 0.0001).

Discussion

The ubiquitin-proteasome system (UPS) plays a crucial role in protein quality control and the regulation of cellular homeostasis. Increasing evidence highlights the UPS as a focal point and therapeutic target in cancer treatment.³¹ However, previous findings have often been confined to the role of individual ubiquitination regulators in prognostic efficacy or chemotherapy resistance, lacking a comprehensive and integrated investigation of UPS regulators.³² In our research, we integrated known regulatory factors within the UPS, identified differentially expressed UPS regulators, and established a UPS-based model to predict outcomes and immune microenvironment in thyroid carcinoma (THCA) patients.

Utilizing the expression profiles of UPS regulators in THCA and machine learning techniques, we developed a six-gene prognostic risk signature comprising KCNS3, DTX4, ZBTB21, KCTD16, FBXO15, and KCNA1. This signature demonstrates strong predictive capability regarding the prognosis of THCA patients. KCNS3, a member of the voltage-gated potassium channel modifier subfamily S, has been implicated in the progression of esophageal cancer.³³ Previous studies confirmed that Deltex E3 ubiquitin ligase 4 (DTX4) acts as a regulator of both the Notch signaling pathway and the interferon (IFN)-I signaling pathway in innate immunity, promoting thyroid carcinoma progression via stearoyl-CoA desaturase 1.^34–36 ZBTB21, a zinc-finger protein, may function as a DNA-binding motif in transcriptional regulation, closely associated with cancer development.^37,38 KCTD16, containing 16 potassium channel oligomerization domains, predominantly operates upstream or within the regulatory scope of G protein-coupled receptor signaling pathways.³⁹ FBXO15, an F-box protein, suppresses breast cancer progression through the EGFR/STAT3/miR-92a-3p axis.⁴⁰ KCNA1, a member of the voltage-gated potassium channel subfamily A, promotes tumor progression by regulating mitochondrial damage and ferroptosis while limiting tumorigenesis through potassium channel-dependent senescence.^41–43 Collectively, these six significant UPS regulators play direct or potential roles in cancer progression. Additionally, based on the aforementioned six-gene signature, we identified two distinct UPS-associated subtypes of THCA. Notably, these subtypes effectively differentiate the prognostic outcomes and varying clinical-pathological features of THCA patients. In summary, our research characterizes the relationship and features between THCA, its subtypes, and the UPS signature, potentially providing novel insights into personalized treatment strategies for THCA.

During tumorigenesis and evolution, dysregulation of the UPS manifests in various biological functions, including signal transduction, immune infiltration, and other critical biological processes.⁴⁴ Further pathway enrichment analyses revealed that the UPS signature is significantly positively correlated with pathways related to oxidative stress response, epithelial-mesenchymal transition, vascular endothelial growth factor receptor binding, positive regulation of protein ubiquitination, and ubiquitin ligase binding. Conversely, significant negative correlations were observed with pathways involving UPS antigen processing and presentation, B cell-mediated immune responses, MHC protein complexes, Th17 cell differentiation, and Th1 and Th2 cell differentiation. These results suggest that the UPS signature regulates the progression of THCA through interactions with cancer-related pathways and immune-related pathways.

The UPS plays a pivotal role in the tumor immune microenvironment by regulating critical aspects of both innate and adaptive immunity, including signaling pathways, protein homeostasis, and antigen processing.⁴⁵ Through pattern recognition receptors (PRRs), the UPS modulates intracellular signaling triggered by danger signals, thus influencing innate immune responses.⁴⁶ Additionally, the UPS affects adaptive immunity by altering levels and functions of T cell signaling intermediates or influencing transcription factors that drive CD4⁺ T cell identity and function.⁴⁷ Notably, emerging evidence underscores the UPS’s significant regulatory role in the differentiation of Th1 and Th2 cells. The E3 ubiquitin ligase CUL4B inhibits Th1 and Th2 differentiation by ubiquitinating H2AK119.⁴⁸ Our study reveals a marked negative correlation between UPS signature and myeloid-derived immune cells, including NK cells, dendritic cells, and macrophages, suggesting potential immunosuppressive features of UPS signature. Interestingly, UPS signature is positively correlated with T cells, helper T cells, and Th2 cells, but negatively correlated with Th1 cells, indicating that UPS signature may reshape the tumor microenvironment by modulating the Th1/Th2 balance, thereby promoting immune suppression and escape. Furthermore, high-risk UPS signatures often display a “hot” immune state with highly active interactions and cross-talk among immune cells. In THCA subtypes, the worse prognosis of the C1 subtype was associated with an immune infiltration profile similar to that of high UPS signature states. In summary, the UPS signature may contribute to the progression of THCA by creating an immunosuppressive microenvironment through immune infiltration.

As research progresses, the UPS has been recognized as a key regulator of resistance to various therapeutic modalities, including chemotherapy, targeted therapy, immunotherapy, and radiotherapy.⁴⁵ Targeting the UPS presents an effective strategy to counteract resistance across different cancer types, thereby enhancing therapeutic approaches and improving patient outcomes.⁴⁹ Our initial screening identified KCNA1 as having significant clinical value. Further investigation into KCNA1’s impact on immunotherapy and chemotherapy sensitivity revealed that high KCNA1 expression was associated with higher tumor immune dysfunction and exclusion (TIDE) scores and lower response rates to immunotherapy compared to low KCNA1 expression, confirming its detrimental effect on immunotherapy. Additionally, bioinformatics and experimental analyses demonstrated that silencing KCNA1 expression effectively reduced the IC50 of gemcitabine, alleviating gemcitabine resistance in THCA patients. Thus, targeting KCNA1 offers a novel therapeutic avenue for THCA treatment.

It is widely acknowledged that E3 ubiquitin ligases, the proteasome, and deubiquitinating enzymes have emerged as novel targets for anticancer therapeutics.^31,50 Nonetheless, the development of UPS inhibitors in clinical trials remains limited.³¹ The underlying mechanisms through which these components regulate cellular processes are intricate, with much of the previous research focusing primarily on apoptosis and autophagy.^51,52 Presently, there is an urgent need for exploration into the mechanisms governing UPS immunological functions, as well as the development of anticancer drugs that target these processes.⁸ Our recent findings unveil a novel interplay between UPS prognostic risk signatures and thyroid immunity, offering a glimmer of hope for the development of targeted UPS inhibitors. Notably, these UPS signatures may reshape the immunological microenvironment by modulating the Th1/Th2 ratio, thereby fostering an immunosuppressive milieu. In essence, our study not only underscores the untapped potential of UPS as a therapeutic target but also provides a compelling rationale for future investigations into the interplay between UPS and immune regulation in THCA.

However, this study has limitations. Due to the lack of survival data for THCA patients in the GEO database, we only validated the UPS signature at the expression level without prognostic validation. Additionally, future research should include functional experiments to elucidate the mechanisms of KCNA1. Lastly, the reliability of this prognostic feature should be improved through iterative long-term clinical application.

In conclusion, the UPS signature accurately predicts the prognosis of THCA patients and can serve as a decision-making tool for guiding their treatment. Furthermore, KCNA1, a key gene within the UPS signature, holds clinical potential in THCA therapy.

Conclusion

Research on the ubiquitin-proteasome system (UPS) in thyroid carcinoma remains limited. Our study elucidated the UPS landscape in thyroid carcinoma by developing a UPS-mediated model to predict patient survival accurately and to delve into the tumor immune microenvironment. We subsequently identified two UPS-associated subtypes by evaluating ubiquitin-proteasome system prognostic model score (UPS-PMS) and their relationship with various clinical parameters and immune responses, which provide insights into patient prognosis and immune profiles. Additionally, we validated the clinical relevance of the prognostically significant E3 ubiquitin ligase KCNA1 in immune therapy and chemotherapy sensitivity through bioinformatics and experimental validation. Collectively, our work introduces an accurate and innovative model aimed at enhancing precise diagnosis and personalized treatment for thyroid carcinoma, deepening the understanding of the relationship between thyroid carcinoma and UPS, and potentially setting the stage for future research breakthroughs.