Mutation Signatures of the APOBEC Family in pan-cancer and Association between AMES and Molecular Characteristics in THCA

To investigate the genomic characteristics of 11 APOBEC family genes across various cancer types, we generated bubble plots to visualize their expression profiles across 28 solid tumors using data from The Cancer Genome Atlas (TCGA) (Figure S1A). Among these genes, APOBEC3C exhibited the highest expression levels. Furthermore, we observed that the expression of certain APOBEC genes was significantly elevated in specific cancer types, such as APOBEC3A in head and neck squamous cell carcinoma (HNSC). Pearson’s correlation analysis revealed significant positive correlations among most APOBEC family members, especially within the APOBEC3 subgroup (Figure S1B). Additionally, we compared gene expression profiles between primary tumors and matched normal tissues across 20 cancer types in TCGA, each featuring at least 36 normal tissue samples. This analysis revealed widespread dysregulation of APOBEC gene expression in tumors compared to normal tissues. For example, APOBEC3H and APOBEC2 were markedly upregulated in THCA (Figure S1C-D).

Next, we evaluated the tumor mutational burden (TMB), total TCW mutations, and the APOBEC mutation enrichment score (AMES) in each tumor sample. Among 33 different cancer types, THCA demonstrated particularly notable AMES values (Figure S1E-G). Using thresholds of 1.5 and 2.5, we classified 9,964 pan-cancer samples into three groups: AMES-Low (21.68%), AMES-Medium (19.65%), and AMES-High (58.67%) (Figure S1H). Random forest analysis identified APOBEC2 as a key contributor to TCW mutations and AMES across pan-cancer samples (Figure S1I-J).

Given the importance of TMB as a biomarker for cancer immunotherapy, we evaluated the relationship between TMB, AMEs, and APOBEC gene expression across cancer types. A heatmap summarized these correlations, with dot size reflecting significance. In THCA, a moderate positive correlation between AMES and TMB was observed (Fig. 1A). Further stratified analyses showed that as AMES increased, TMB) and intratumoral heterogeneity (ITH) also rose significantly, while tumor neoantigen burden (TNB) decreased (Fig. 1B). This suggests that while higher AMES levels are associated with greater genetic heterogeneity, they do not necessarily enhance tumor immunogenicity. Additionally, in THCA, AMES was significantly positively correlated with non-synonymous mutations (R = 0.43, p = 2.1e − 07) and synonymous mutations (R = 0.37, p = 1.5e − 09) (Fig. 1C).

Fig. 1

Mutation signatures of the APOBEC family in pan-cancer and association between AMES and molecular characteristics in THCA. A Evaluating the relationship among TMB, AMES, and APOBECs in pan-cancers. B Exhibiting the alterations in TMB, TNB, and ITH scores along with rising AMES levels. C Illustrating the correlation between AMES and non-synonymous and synonymous mutations in Thyroid Carcinoma samples. D Identifying driver genes in the high, medium, and low AMES groups by means of the OncodriveCLUST algorithm. E Demonstrating the change in NRAS mutation frequency with ascending AMES levels. F In contrast to the other two cohorts, the forest plot discloses significant modifier genes identified in the AMES-H group. G The distribution of non-silencing mutations in the three AMES categories is depicted through a ternary plot, with the 10 most prevalent non-silencing mutations emphasized in red. H-J Presenting the co-occurring and mutually exclusive mutations witnessed in the high, medium, and low AMES groups. K The mutation rates of seven DNA damage response (DDR) pathways-specifically, HRR, BER, NER, NHEJ, MMR, FA, and TLS—are summarized across various AMES cohorts. L Signifying that compared to DDR wild-type samples, the AMES of DDR mutant samples is significantly elevated. M An oncoplot diagram depicts the frequency distribution of the ten most commonly mutated DDR genes across different AMES groups. N Describing the DDR pathways related to the 10 most commonly mutated DDR genes

OncodriveCLUST analysis identified NRAS as a recurrent driver gene in the AMES-High and AMES-Moderate groups, but not in the AMES-Low group (Fig. 1D). NRAS mutation frequency progressively declined with decreasing AMES levels (Fig. 1E). Moreover, distinct mutation patterns were found in the AMES-High group. For instance, BRAF emerged as the most frequently mutated gene (Fig. 1F). A ternary plot showed the distribution of all non-silent mutations across the three AMES groups, highlighting the top ten frequent mutations (Fig. 1G). Additionally, the frequency of co-occurring and mutually exclusive mutations increased with higher AMES levels, further linking AMES to enhanced somatic mutational activity in THCA (Figs. 1H - J).

In THCA, DDR gene mutations are prevalent and have been associated with improved outcomes following ICB therapy [32]. We analyzed the mutation frequencies of seven DDR pathways (HRR, BER, NER, NHEJ, MMR, FA, TLS) across various AMES groups. The homologous recombination repair (HRR) repair pathway exhibited the highest mutation frequency, followed by the base excision repair (BER) pathway (Fig. 1K). In the AMES-high group, DDR pathway mutations were significantly more frequent, and samples of DDR mutations had notably higher AMES scores compared to DDR wild-type samples (Fig. 1L). We mapped the distribution of the top 10 frequently mutated DDR genes across AMES groups (Fig. 1M), including PRKDC and DCLRE1C involved in the NHEJ pathway, BRCA1 and BRCA2 involved in the FA and HRR pathways, PARP4 in the BER pathway, SLX4 in the HRR pathway, and CUL3, ELOA2, ERCC3, and ERCC6 involved in the NER pathway (Fig. 1N).

APOBEC Mutagenesis and its Correlation with AMES in Thyroid cancer Mutation Patterns

To better understand the mutational features of thyroid cancer, two whole exome sequencing (WES) cohorts were analyzed. Using Non-Negative Matrix Factorization (NMF) with an optimal factorization parameter of k = 5, We identified five unique mutational signatures in the TCGA-THCA cohort. Notably, Signature 1, characterized by cytidine deaminase activity (C > T), was attributed to APOBEC family mutagenesis (Fig. 2A). The relative occurrence rate of each mutation signature in the TCGA cohort showed a predominance of APOBEC Signature 1 (Fig. 2B). Furthermore, the abundance of the APOBEC signature was strongly positively correlated with AMES (Fig. 2C), and the number of TCW mutations significantly increased with rising AMES levels (Fig. 2D). Consistent results were observed in the independent DFCI/MSKCC-THCA cohort (Figs. 2E-H). In summary, these findings demonstrate that APOBEC-driven mutagenesis plays a major role in shaping the mutation landscape of THCA and is strongly positively correlated with AMES.

Fig. 2

APOBEC Mutagenesis and Its Correlation with AMES in Thyroid Cancer Mutation Patterns. A In the TCGA-THCA cohort, five distinct mutation signatures were discerned by employing an optimal factorization k value (k = 5) through the Non-negative Matrix Factorization (NMF) algorithm. B A pie chart represents the relative frequency of each mutation signature within the TCGA cohort. C Correlation analysis was conducted between the abundance of APOBEC signatures and AMES. D Changes in TCW mutations along with increasing AMES levels. E Five mutation signatures were identified in the DFCI/MSKCC - THCA cohort through using the optimal factorization k value (k = 5) in the NMF algorithm. F The relative abundance of each mutational feature within the MSKCC cohort is displayed in the form of a pie chart. G The correlation of the abundance of APOBEC signatures with AMES in the DFCI/MSKCC - THCA cohort. H Changes in TCW mutations in the DFCI/MSKCC - THCA cohort with rising AMES levels

In-depth Analysis of APOBEC Family Expression Patterns in THCA and their Effect on Malignant Tumor Cell Phenotypic Remodeling

To deeply investigate the expression patterns of the APOBEC family in THCA and its microenvironment, we analyzed two single-cell RNA sequencing (scRNA-seq) datasets (GSE191288 and GSE232237). The GSE191288 dataset includes three pairs of bilateral PTC samples and one non-tumorous thyroid sample. The distribution of these samples is shown in Fig. 3A. Using t-SNE dimensionality reduction, we visualized the distribution and disparities among various cell types (Fig. 3B).

Fig. 3

In-depth analysis of APOBEC family expression pattern in THCA and its effect on the phenotypic remodeling of malignant tumor cells. A Sample distribution in the single-cell dataset GSE191288. B The distribution state of each cell subpopulation in the single-cell GSE191288 was clarified through t-SNE dimensionality reduction. C Displays the expression status of APOBEC2 in each cell subpopulation. D Presents the expression situation of APOBEC3H in each cell subpopulation. E The heatmap reveals the expression situation of marker genes in each cell population. F With stromal cells as a reference, epithelial cells were separated from the total number of epithelial cells via extrapolated large-scale CNV. G The epithelial cell subpopulations were further delineated by utilizing t-SNE dimensionality reduction. H Different colors were employed to represent the expression level of APOBEC2 in the epithelial cell subpopulations. I The box plot clearly shows the normalized expression level of APOBEC2 in each cell subpopulation. J Illustrates the pseudo-time trajectory of epithelial cell evolution, and the arrow indicates the assumed developmental direction. K Reflects the dynamic changes of APOBEC2 in the pseudo-time trajectory of epithelial cells. L The ssGSEA and random forest algorithms jointly indicated that the RAS pathway is the most significant cancer biomarker among the five subgroups. M Explains the correlation between APOBEC2 and the RAS score in all malignant epithelial cells. N From normal thyroid tissue to adjacent tissues to THCA, APOBEC2 exhibits a progressively and markedly increasing trend. O The cumulative proportion curve reveals that the elevated APOBEC2 cohort, depicted by the red curve, is consistently situated to the right of the diminished APOBEC2 cohort, represented by the blue curve. This finding highlights the significant impact of APOBEC2 in catalyzing TCW mutations within THCA

We then assessed the expression patterns of 11 APOBEC family members across different cell types. In THCA, APOBEC3H was underexpressed in epithelial cells, while APOBEC2 exhibited specific expression patterns (Fig. 3C-D). The heatmap illustrates the expression level of marker genes in each cell subgroup (Fig. 3E), and the expression distribution of each marker across each cell subgroup, with t-SNE plots detailing the distribution of each marker (Figure S2). By inferring large-scale copy number variations (CNVs) and using stromal cells as a reference, we identified epithelial cell_c1 as malignant due to increased chromosomal instability (Fig. 3F). Further, t-SNE analysis subdivided malignant epithelial cells into six clusters (C_0, C_1, C_2, C_3, C_4, C_5) (Fig. 3G), with APOBEC2 expression notably visualized across these subgroups (Fig. 3H). Boxplot analysis revealed that APOBEC2 expression was almost absent in normal epithelial cells but significantly elevated in clusters C_0 and C_5 (Fig. 3I).

Next, we performed pseudo-time trajectory analysis to map the developmental pathways of epithelial cells, with arrows indicating potential progression (Fig. 3J). Dynamic changes in APOBEC2 expression across this trajectory are shown in Fig. 3K.

To explore biological heterogeneity among malignant subgroups, we integrated ssGSEA and a random forest algorithm to identify key cancer markers. The RAS signaling pathway emerged as the most influential (Fig. 3L). A significant positive correlation was observed between APOBEC2 expression and the RAS score across malignant epithelial cells (R = 0.48, p = 3.8e − 07; Fig. 3M). Additionally, APOBEC2 expression progressively increased from normal thyroid tissue to adjacent tissues, and even further in THCA samples (p = < 2e − 16; Fig. 3N). The cumulative distribution plot shows that TCW mutations are more prevalent in the high APOBEC2 group, suggesting a key role in the mutagenesis of THCA (Fig. 3O).

Exploring the Relationship Between Tumor-derived APOBEC Family Signaling and Immune Cells in THCA

The GSE232237 dataset encompasses single-cell RNA sequencing (scRNA-seq) data of 7 cases of PTC and 5 cases of ATC. Figure 4A summarizes the distribution of different cell types. Using t-SNE dimensionality reduction, we mapped the distribution of each sample across different cell subgroups (Fig. 4B). The bubble chart in Fig. 4C displays marker gene expression across cell subgroups, while Figure S3 further illustrates marker distribution via t-SNE.

Fig. 4

Exploring the relationship between tumor-derived APOBEC family signaling and immune cells in THCA. A t-SNE dimensionality reduction reveals the distribution of each cell subpopulation within the single-cell dataset GSE232237. B t-SNE further demonstrates the distribution of samples within the single-cell dataset GSE232237. C The bubble plot exhibits the expression level of the marker gene of each cell subpopulation. D The histogram presents the expression of APOBEC2 in each cell subpopulation. E The histogram depicts the signalling of the APOBEC family to tumour cells. F t-SNE dimensionality reduction clustering presents the distribution of the subpopulations of CD8+ T cells. G Illustrates the pseudo-time trajectory of CD8+ T cell development, and the arrow indicates the hypothesised development direction. H alterations in the differentiation score triggered by the APOBEC family during the development of CD8+ T cells. I The relevant heat map elucidates the relationship between the APOBEC family and the infiltration of diverse immune cells in a multitude of THCA samples. J The cumulative proportion curve and (K) the box plot jointly signify that the M2 infiltration abundance in THCA samples with a high APOBEC signal group is notably enhanced

We specifically assessed the expression of APOBEC family members across cell types. Bar graph analysis indicated elevated normalized APOBEC expression levels in tumor cells compared to other cell types (Fig. 4D). Moreover, the APOBEC family signal from tumor cells appeared to primarily act on CD8+ T cells (Fig. 4E), suggesting a potential link between high APOBEC expression and immune modulation in the tumor microenvironment.

Subsequent t-SNE clustering divided CD8+ T cells into five subpopulations (Fig. 4F). Pseudo-time trajectory analysis indicated two main developmental branches forming CD8+ Tex_c1 and CD8+ Tex_c2 subclusters (Fig. 4G). Dynamic expression profiling revealed significant enrichment of tumor-induced APOBEC signals in the terminal stages of CD8+ T cell differentiation (Fig. 4H).

We also generated a heatmap illustrating the correlation between APOBEC family expression and immune cell infiltration (based on CIBERSORT analysis of TCGA data; Fig. 4I). A strong positive correlation was identified between APOBEC signatures and M2 macrophages (r > 0.6, p < 0.00001). Cumulative proportion curves (Fig. 4J) and boxplots (Fig. 4K) further confirmed significantly increased M2 macrophage infiltration in THCA samples with high APOBEC expression. Collectively, these findings from both single-cell RNA and bulk RNA sequencing data imply that the APOBEC family may influence CD8+ T cell differentiation and immune suppression in THCA.

Biological Characteristics and Prognostic Implications of AMES in THCA

To explore the biological characteristics associated with AMES in THCA, we employed the fgsea algorithm across all GO_BP gene sets. Samples with high AMES scores exhibited significantly enhanced activity in several immune-suppressive pathways, most notably the “IMMUNE RESPONSE-INHIBITORY SIGNALING PATHWAY” (NES = 2.13, p = 1.7e − 3; Fig. 5A).

Fig. 5

Exploring the biological features and prognostic relevance of AMES in THCA. A The fgsea algorithm was employed on all GO_BP pathways in distinct AMES groups, and 10 pathways with significant variations in the high AMES group were identified. B For multiple samples, the association network between AMES and various immune cell types was depicted. C Changes in immune infiltration scores upon an increase in AMES levels. D Alterations in the abundance of CD8+ T cells as AMES levels increase. E As illustrated in the heat map, to assess the predictive capacity of AMES for cancer immunotherapy, the expression profiles of T cell inflammatory genes (18 genes) related to the ICB response were introduced. F Comparison of Trags cells between the high AMES group and the low AMES group. G - I GSEA analysis was executed on three related yet independent genes, and the IFN-γ response in the high AMES group was significantly diminished. J Comparison of IFNG mRNA expression (log 2 normalization) between the high AMES group and the low AMES group. K Significant decrements in the expression levels of the typical immune checkpoint markers namely PD-1, PD-L1, TIGIT, and ctla − 4 were observed in the high AMES group. L Comparison of Foxp3+T scores between the high AMES group and the low AMES group

A network analysis was conducted to examine the correlations between AMES and immune cell types (Fig. 5B). AMES showed negative correlations with activated dendritic cells (DCs), eosinophils, M0, CD8+ T cells, and CD4+ Naive cells, but a strong positive correlation with neutrophils. As AMES levels rose, the immune infiltration score (Fig. 5C) and CD8+ T cell abundance (Fig. 5D) significantly declined.

Further analysis of immune differences between low and high AMES groups revealed a reduction in T cell inflammatory gene expression (related to ICB response) in the high AMES group (Fig. 5E). Conversely, regulatory T cells (Tregs) were significantly enriched (Fig. 5F). GSEA of three distinct IFN-γ-related pathways indicated substantial reductions in IFN-γ responses in high AMES samples (Figs. 5G - I), with corresponding reductions in IFN-g mRNA expression (log2 normalization) (Fig. 5J).

Given the reliance of immune checkpoint blockade (ICB) therapies on suppressing key checkpoints, we assessed the expression of PD-1, PD-L1, TIGIT, and CTLA-4, all of which were markedly decreased in the high AMES group (Fig. 5K). Interestingly, the Foxp3+ T score, reflecting T reg-mediated suppression, increased significantly (Fig. 5L). Collectively, these findings suggest that AMES is associated with immune suppression in THCA.

Finally, we assessed the prognostic relevance of AMES in THCA. Kaplan-Meier survival analysis revealed that THCA patients with higher AMES scores had significantly poorer overall survival (HR = 5.38, CI = 1.95–14.81, p = 0.001; Figure S4).

Establishing an APOBEC Mutation Model Related To the Prognosis of THCA Patients

Given the predictive value of AMES for survival outcomes, we aimed to develop an OS risk assessment model based on the gene expression profiles associated with AMES.

First, we identified 327 DEGs through differential gene analysis (with an FDR q-value < 0.05). Next, univariate Cox regression analysis of these DEGs yielded 91 candidate genes significantly associated with OS (p < 0.05). After applying LASSO regularization (Fig. 6A), we retained 27 genes with nonzero Cox coefficients (Fig. 6B), enabling the calculation of a prognostic APOBEC mutagenesis-related risk score (AMrs) for each THCA patient (detailed in the Methods section), and AMrs address the inherent limitations of DNA-level AMEs assessments in predicting functional consequences (Figure S5). A spine plot revealed distinct differences in the expression of cancer markers between low- and high- AMrs groups (Fig. 6C). In TCGA-THCA training cohort, patients with elevated AMrs exhibited significantly worse survival outcomes (Fig. 6D). The prognostic value of AMrs was validated in three independent THCA cohorts: CPTAC (HR = 4.63, CI = 3.55–6.12, p = 0.008), QCMG (HR = 2.01, CI = 1.75–5.4, p = 0.002), and UTSW (HR = 3.1, CI = 1.51–8.6, p = 0.036) (Fig. 6E-G). These results indicate that AMrs is a robust and optimal tool for predicting prognosis in THCA patients.

Fig. 6

Establishing an APOBEC Mutation Model Related to the Prognosis of THCA Patients. A-B After the application of LASSO regularization, a total of 27 genes retained their Cox coefficients, facilitating the calculation of the prognostic risk score related to APOBEC mutagenesis (AMrs) for each patient with THCA. C A ridgeline plot reveals significant differences in the efficacy of various cancer markers when comparing samples classified as low AMrs with those identified as high AMrs. D The survival outcomes for both the low and high AMrs subgroups were evaluated using the training dataset (TCGA-THCA). E - G Survival prognosis of the high and low AMrs subgroups in the three independent THCA validation sets. H To identify potential therapeutic targets and drugs for patients with elevated AMr levels, a comprehensive cumulative assessment of 1927 compounds was conducted utilizing three independent drug response databases (particularly GDSC, CTRP, and PRISM). I Signal pathways and therapeutic targets of 13 candidate compounds in GDSC. J & L The projected AUC values for CTRP and PRISM compounds were collected for each sample from The Cancer Genome Atlas (TCGA), and a Spearman’s correlation analysis between the AUC values and AMr was performed. The eight compounds which manifested the most pronounced negative correlation coefficients in CTRP and PRISM are depicted in the dashed graph. K & M All estimated AUC values of these compounds were compared between the high AMr group and the low AMr group

We further explored drug sensitivity differences between high- and low-AMrs groups to uncover potential therapeutic targets and compounds. A total of 1927 compounds were screened across three distinct drug response databases: GDSC, CTRP, and PRISM. Starting with GDSC data, we estimated IC50 values for 288 compounds per TCGA sample and then carried out Spearman correlation analysis between AMrs and the IC50 values. Using a selection threshold of a negative correlation coefficient and p ˂ 0.05, we successfully identified 15 candidate compounds. Among them, the RSK inhibitor SL0101 had the strongest negative correlation with AMrs (Fig. 6H). A detailed analysis of signaling pathways and therapeutic targets for these compounds is shown in Fig. 6I. Similarly, we analyzed AUC values for compounds in the CTRP and PRISM databases. Eight compounds with the strongest negative correlations with AMrs were identified in each dataset (CTRP: PI-103, cucurbitacin I, SL0101, vincristine, GSK461364, SB-743921, BI −2536, lipomycin B; PRISM: SL0101, vincristine, tilmicosin, rigosertib, verubulin, PHA-848125, NVP -AUY922 and vindesine) (Figs. 6J & L). Their AUC values were significantly lower in high-AMrs patients (Figs. 6K & M). These findings highlight several compounds with potential therapeutic relevance for THCA patients at high risk.

Role of SL0101 in Modulating Malignant and the Immune Microenvironment of THCA

In the FTC-238 thyroid cancer cell line, SL0101 treatment significantly decreased the mRNA and protein levels of metastasis-related genes (Snail1, Vimentin, N-cadherin, and ZEB1) (Fig. 7A-B). Additionally, CCK8 assays demonstrated that SL0101 suppresses tumor cell proliferation (Fig. 7C). To evaluate SL0101’s effect on the immune microenvironment, we co-cultured FTC-238 cells with PBMCs isolated from mouse spleens. Flow cytometry revealed that SL0101 treatment decreased the proportion of Tregs (CD4 + CD25 + Foxp3 + T cells) and increased CD8 + T cell proportions (Fig. 7D - F). In vivo, SL0101 significantly reduced tumor burden in C57BL/6J mice compared to controls followed by treatment with SL0101 or a control vehicle, and tumors were harvested after 28 days (Fig. 7G-I). Immunohistochemical analysis of harvested tumors confirmed the decrease in Ki67 expression (a typical proliferation marker) and an increase in TUNEL positivity (an apoptosis marker) following SL0101 treatment (Fig. 7J). Flow cytometry further validated the reduction of Tregs and enhancement of CD8 + T cells.

Fig. 7

Role of SL0101 on the malignant phenotype and microenvironment of thyroid cancer. A The mRNA expression levels of Snail1, vimentin, N-cadherin, and ZEB1 in each group were determined by qPCR. B The protein expression levels of Snail1, vimentin, N-cadherin, and ZEB1 in each group were detected through Western blot. C CCK-8 assay was carried out to ascertain the cell viability following SL0101 treatment in FTC-238. D Schematic illustration of the thyroid cancer cell-PBMCs co-culture system. E - F The proportions of Tregs (E) and CD8+ T cells (F) were determined by flow cytometry. G Schematic representation of the in vivo experimental design. H Growth curves were plotted every three days based on the tumor size. I Tumor weight was analyzed when the tumors were harvested 28 days after subcutaneous implantation. J Immunofluorescence staining of Ki67 and TUNEL in subcutaneous tumors of each group. Representative images were presented. Ki67+ and TUNEL+ tumor cells were analyzed using Image-J software. K - L The proportions of tumor-infiltrating Tregs (K) and CD8+ T cells (L) in subcutaneous tumors were determined by flow cytometry

Effect of SL0101 on Aerobic Glycolysis and FOXP3 + T cell-mediated Remodeling of the Immune Environment in THCA

Given the suspected involvement of APOBEC2 in cellular energy metabolism, we investigated SL0101’s impact on aerobic glycolysis. SL0101 downregulated both mRNA and protein levels of key glycolytic biomarkers (MCT1, MCT4, GLUT1, GLUT4, LDHA, HK2, PKM2, and PFKP) (Fig. 8 A-B) and significantly reduced lactate production(Fig. 8C). It also decreased the extracellular acidification rate (ECAR) while increasing the oxygen consumption rate (OCR), suggesting a shift toward oxidative phosphorylation (Figs. 8D-E). To further probe the relationship between glycolysis and Treg infiltration, we blocked glycolysis using the LDHA inhibitor GSK2837808A. Co-culturing PBMCs with thyroid cancer cells revealed that SL0101 decreases Treg proportions and increases CD8 + T cell proportions, which is consistent with previous findings. However, co-treatment with GSK2837808A abolished these immune changes (Figs. 8F - H). In vivo experiments using the mouse thyroid cancer cell line TBP3743 corroborated that SL0101’s enhancement of the immune microenvironment is glycolysis dependent (Figs. 8I - O).

Fig. 8

Effect of SL0101 on aerobic glycolysis and FOXP3 + T cell-mediated remodeling of the immune environment in thyroid cancers. A The expression levels of glycolysis-related biomarkers in each group were detected via qPCR. B The expression levels of glycolysis-related biomarkers in each group were determined through Western blot. C The relative lactate production level in each group was measured. D ECAR between each group was compared, and glycolytic rate and glycolytic capacity were analyzed. E OCR between each group was compared, and basal respiration and maximum respiration levels were evaluated. F Schematic illustration of the thyroid cancer cell-PBMCs co-culture system treated with SL0101 or GSK2837808A. G - H The proportions of Tregs (G) and CD8 + T cells (H) were determined by flow cytometry. I Schematic representation of the in vivo experimental design. J Growth curves were plotted every three days based on the tumor size. K Tumor weight was analyzed when the tumors were harvested 28 days after subcutaneous implantation. L - M The proportions of Ki67 + (L) and TUNEL + (M) tumor cells in each group were assessed. N - O The proportions of tumor-infiltrating Tregs (N) and CD8 + T cells (O) in subcutaneous tumors were determined by flow cytometry

Synergistic Therapeutic Potential of SL0101 and FOXP3 Inhibitor in THCA

Since SL0101 reduces Treg infiltration, we postulated that it could sensitize tumors to Treg-targeted therapies. We tested the synergy between SL0101 and the FOXP3 inhibitor Alkyne-P60 using Synergyfinder-2.0. The combination achieved a synergy score of 23.399, indicating strong synergistic cytotoxicity (Fig. 9A). In vivo thyroid cancer-bearing C57BL/6J mice were randomized to receive vehicle, SL0101, Alkyne-P60, or their combination (Fig. 9B). Combination treatment significantly reduced significantly reduced tumor burden and proliferation, increased tumor apoptosis, and significantly altered the proportion of tumor-infiltrating immune cells compared to single treatments (Fig. 9C-G). These findings highlight the promising therapeutic potential of combining SL0101 with FOXP3 inhibition for THCA. Consistently, the safety profile of SL0101 in conjunction with Alkyne-P60 was deemed acceptable, as confirmed by the assessment of hepatorenal toxicity, peripheral blood serum markers, and histopathological characteristics of various organ tissues (Fig.S6).

Fig. 9

Exploring the Synergy of SL0101 and FOXP3 Inhibitor in Thyroid Cancer Therapeutics. A The synergy score of SL0101 and Alkyne-P60 was analyzed using Synergyfinder. B Schematic illustration of the in vivo experimental design. C Growth curves were plotted every three days based on the tumor size. D Tumor weight was analyzed when the tumors were harvested 28 days after subcutaneous implantation. E Immunofluorescence staining of Ki67 and TUNEL in subcutaneous tumors of each group. Representative images were presented. Ki67 + and TUNEL + tumor cells were analyzed by Image-J software. F-G The proportions of tumor-infiltrating Tregs (F) and CD8 + T cells (G) in subcutaneous tumors were determined by flow cytometry