Exploring the relationship between RMRs and treatment response in Pan-cancer: insights from single-cell RNA datasets

The present study outlines our investigation into the relationship between RMRs and cancer treatment response. We conducted an extensive analysis, as illustrated in Fig. 1, which delves into the main findings of our research. To gather data, we curated seven single-cell RNA datasets that encompassed neoadjuvant therapy information. Details of these datasets can be found in Table S1.

Fig. 1

The whole flow chart and the analysis content schematic

Utilizing established markers (Fig. 2A and Fig. S1), we classified tumor parenchymal and microenvironmental cells into distinct categories, including epithelial cells, T cells, B cells, myeloid cells, mast cells, fibroblasts, endothelial cells, pericytes, a few endocrine cells, and Schwann cells. Notably, except for breast cancer, wherein the analysis focused solely on tumor cells, we observed varying proportions of these cell types between responders and non-responders (Fig. 2B). To further investigate the molecular mechanisms underlying treatment response, we conducted differential expression analysis between the two groups and performed pathway enrichment analysis on the differentially expressed genes. Interestingly, we discovered that these genes consistently participated in RNA methylation modification (Fig. 2C-D). Through an extensive literature search, we identified a total of 46 RMRs, comprising 23 m6A regulators, 6 m1A regulators, 13 m5C regulators, and 4 m7G regulators. Among these, the m6A regulators, particularly FTO and FMR1, have been extensively studied (Fig. 2E, Table S2).

Fig. 2

Differences between treatment responders and non-responders are associated with RMRs. (A) UMAP plot displaying diverse cell types that have been identified. (B) Bar plot showing the percentage of different cell subpopulations in treatment responders and non-responders. (C) Volcano illustrating different expression genes between responders and non-responders. (D) Demonstrating pathways for differential gene enrichment with P < 0.05. (E) 46 RNA methylation regulators classified as writers, readers, and erasers from m6A, m5C, m1A, and m7G. (F) Heatmap displaying the correlation between RMRs expression and drug IC50 values of the GDSC database

Moreover, we evaluated the expression patterns of these 46 RMRs across approximately 1000 cell lines and found a significant association between their expression levels and drug sensitivity. This crucial finding further highlights the indispensable role played by RMRs in drug therapy (Fig. 2F).

Comprehensive analysis reveals landscape of genetic and transcriptional aberrations in RMRs across diverse cancer types

Subsequently, an in-depth analysis of 38 frequently occurring RMRs that showed copy number variations (CNVs) in at least one cancer type was conducted. Notably, we observed that the m6A writer WTAP was located within deletion peaks in 12 cancer types, whereas the m5C reader ALYREF was amplified in 7 cancer types (Fig. 3A). Among a comprehensive dataset of 9991 TCGA samples, we found that m6A writers and readers exhibited the highest frequencies of CNVs (Fig. S2).

Fig. 3

Genomic aberrations, transcriptional alterations, proteomic discrepancy, and survival relevance of RMRs across cancer types. (A) Bubble plot showing the Gscore and CNV status of RMRs across cancers. The sizes and colors of the bubbles represent Gscore and CNV status, respectively. (B) Heatmap depicting the somatic mutation frequencies of each RMR in each cancer type. (C) Waterfall plot showing somatic mutations of RMRS in 2876 tumor samples. The mutation frequencies were shown on the right barplot, and the cancer types of samples were shown at the bottom. (D) Dot plot indicating the differentially expressed RMRs between tumor and adjacent normal tissues in each cancer type, with up or down-regulating and the P-value being annotated. (E) The bar plots of RT-PCR showing mRNA expression levels of RMRs in tumor and normal tissues (n = 6). * P < 0.05, **P < 0.01 or *** P < 0.001. (F) The protein expression of RMRs according to HPA. (G) Hazard ratios of overall survival between high and low expression groups regarding each RMRs in each cancer type, with the size and color of the bubble denoting the P-value and Hazard Ratio (HR) of overall survival, respectively

Furthermore, the non-silent somatic mutations in RMRs were investigated. These mutations were particularly prevalent in specific cancer types, including skin cutaneous melanoma (SKCM), bladder urothelial carcinoma (BLCA), and uterine corpus endometrial carcinoma (UCEC) (Fig. 3B). Out of the total 9834 pan-cancer samples examined, 2876 (29.24%) samples harbored at least one mutation in an RMR. Interestingly, among these regulators, the m5C erasers displayed the highest mutation frequency across the 2876 samples (Fig. 3C).

To explore the transcriptional levels of RMRs, we compared tumor tissues with adjacent normal tissues obtained from the TCGA and GTEX databases. Despite being frequently affected by CNVs leading to deletions, most RMRs exhibited significant upregulation in tumor tissues. This finding suggests that there is a transcriptional upregulation mechanism governing RMRs in cancers (Fig. 3D). We further validated these findings using RT-PCR, which confirmed the upregulation of most RMRs in tumor tissues across various cancer types (Fig. 3E, Table S3). Additionally, immunohistochemical data from the HPA supported our observations, indicating higher protein levels of RMRs in tumor tissues compared to normal tissues (Fig. 3F).

Moreover, we investigated the prognostic relevance of RMRs by conducting survival analysis across diverse cancer types. Our results revealed that RMRs had a broad impact on cancer survival, with kidney cancers exhibiting a particularly significant association (Fig. 3G). Overall, our comprehensive analysis highlights the crucial role of RMRs in regulating tumor progression.

Machine learning constructed distinct regulatory RMR clusters in Pan-cancer and their association with genomic alterations

Subsequently, a total of over ten thousand samples were classified into three distinct RMR clusters through unsupervised clustering after normalization: Cluster1 (3388 samples), Cluster2 (2346 samples), and Cluster3 (4593 samples) (Fig. 4A, Table S4). The three RMRs clusters were clearly distinguished based on the principal component analysis (PCA) results and the distribution of the clusters varied across specific cancer types, with Cluster3 being the predominant group in most cases (Fig. 4B-C). We utilized the normalized expression matrix of 46 RMRs from TCGA samples to train a random forest model. The predictive accuracy of the model was determined to be 0.91, and the area under the receiver operating characteristic curve in the tenfold cross-validation exceeded 0.9 for the classification of the three pairs of clusters. These results signified that the model demonstrated both accuracy and robustness (Fig. 4D-E). To assess the value of each regulator’s contribution to the model, we calculated its importance score. Notably, the m5C writer NOP2 ranked first among these regulators and was reported to methylate the C(5) position of cytosine 4447 in 28 S rRNA [42] (Fig. 4F). The expression levels of the regulatory RMRs exhibited differential patterns across the different clusters. Specifically, nine regulators displayed a decreasing trend in expression from Cluster1 to Cluster3, indicating a gradient of high to low expression among these clusters (Fig. 4G). Furthermore, their expression profile in the three clusters of individual cancer types resembled the profile of the pan-cancer clusters (Fig. S3). Subsequently, the analysis of three clusters revealed distinct characteristics. Cluster1 displayed a higher frequency of oncogene mutations, particularly TP53, and showed genomic instability with copy number variations. This cluster was also associated with increased genetic alterations and enriched pathways related to cell proliferation and carcinogenic activation, suggesting a more aggressive tumor phenotype (Fig. 4H-J).

Fig. 4

RMRs clusters construction and different genomic characteristics. (A) Consensus matrix plot of 46 RNA methylation regulators in pan-cancer. (B) Scatter plot of PCA showing distinct three clusters. (C) Bar plot showing the percentage of samples of each cluster in each cancer type. (D) Confusion matrix of predicted RMRs clusters and true RMRs clusters (E) ROC curves of ten-fold cross-validation tests performed on the model built with RMR expression data of the TCGA samples. (F) Lollipop plot showing the importance scores of 46 RMRs contributing to the trained random forest model. (G) Heatmap of expression levels of 46 RMRs in three clusters. (H) Waterfall plot depicting the top 15 genes with the highest mutation rates. (I) Heatmap of CNV across all cancer types in three clusters, blue means copy number loss while red means copy number gain. (J) Bubble plot visualizing the tumor progression-related features and Hallmark pathways from MSigDB. The size of the bubble denoted the P-value. Prune and blue denoted higher and lower pathway enrichment scores in each comparison

Clinical characteristics and prognosis associated with RMRs clusters across Cancer types

In addition, the TCGA epithelial tumors were categorized into basal-like, luminal A, and luminal B subtypes using the PAM50 clustering algorithm. In our analysis, we compared the distribution of these subtypes within the three RMR clusters. It was observed that Cluster1 was predominantly composed of the basal-like subtype, while Cluster3 mainly consisted of the luminal A subtype. Furthermore, a comparison of the molecular subtypes revealed that Cluster1 was primarily associated with the basal subtype of breast cancer and the chromosomal instability (CIN) subtype of gastrointestinal tumors (Fig. 5A-B, Table S5-S6).

Fig. 5

The clinical characteristics of RMRs clusters. (A) The proportion of PAM50 subtypes in the RMRs clusters. (B) The proportion of typical molecular subtypes in RMR clusters. (C) The proportion of immune subtypes in the three RMRs clusters. (D) Kaplan-Meier curves showing overall survival in different clusters across 10 cancer types respectively. P-value was calculated by the two-sided log-rank test and p < 0.05 was considered statistically significant. (E) The fraction of patients with different clinical stages in three RMRs clusters

The TME subtype, including immune-enriched, nonfibrotic (IE), immune-enriched, fibrotic (IE/F), fibrotic (F), and immune-depleted (D), was defined by 29 functional gene expression signatures [43]. Our analysis further demonstrated that Cluster1 showed a strong association with TME subtype D, while Cluster3 was predominantly associated with TME subtypes F and IE/F. This indicates that Cluster1 exhibits an immunosuppressed phenotype resembling the basal-like subtype, which is known to be associated with poor prognosis. Importantly, patients in Cluster1 exhibited significantly worse overall survival in 10 cancer types, poorer disease-specific survival in 12 cancer types, and shorter progression-free survival in 7 cancer types, highlighting the unfavorable prognosis associated with this cluster. Additionally, Cluster1 displayed a higher proportion of stage III-IV cases in 6 cancer types (Fig. 5C-E, Fig. S4A-B, Table S7-8).

Link of RMRs clusters and cellular components of the tumor microenvironment in pan-cancer

On the other hand, the xCell method to analyze TME-infiltrating immune cells in three clusters was utilized. Cluster1 exhibited higher levels of Th2 cells, naive T cells, and pro-B cells, indicating a weaker tumor suppression ability. In contrast, Cluster3 displayed higher levels of CD4 + T cells, CD8 + T cells, NK cells, B cells, cancer-associated fibroblasts, and endothelial cells, suggesting an abundance of tumor-killing effector cells (Fig. 6A, Table S9). ImmuneScore, StromalScore, and MicroenvironmentScore were calculated to assess the immune and stromal components. Cluster3 demonstrated higher immune and stromal cell activity, while Cluster1 had the lowest MicroenvironmentScore (Fig. 6B). Moreover, we used the “ssGSEA” method to evaluate the infiltration of 28 immune cell types in RMRs clusters. Cluster3 showed increased immune cell infiltration, whereas Cluster1 exhibited lower abundance (Fig. S5A, Table S9).

Fig. 6

Tumor microenvironment components of RMRs clusters and groups. (A) Heatmap showing the abundance scores of immune cell types computed by the xCell algorithm. The color denoted the specific RMR cluster which had a higher abundance score than the other two clusters (B) Box plots showing the differences of the immune cell infiltration scores computed by xCell analysis among different clusters. * P < 0.05, **P < 0.01, *** P < 0.001, ****P < 0.0001 (Kruskal-Wallis test). (C) UMAP plot showing annotated cell types that have been defined. (D) Bar plot showing the proportion of different cell subpopulations in high-expression and low-expression groups and UMAP plot described the distribution of mean expression values across cell types. (E) Characterization of high-expression group and low-expression group differential genes in different cell populations. (F) The number of interactions and weight of interaction in high-expression group and low-expression group (G) Ligand-receptor interaction pairs genes in epithelial cells and other cell types. MIF and MK signaling pathway networks of high and low groups. The thickness of the interworking line segments represents the number of interworking pairs

Subsequently, we examined the expression levels of various immune-related genes among the three RMRs clusters. Most genes showed higher expression in Cluster2 and Cluster3 (Fig. S5B-F). By dividing the five single-cell datasets into high and low groups based on the average expression of RMRs, we observed distinct cell annotations through dimensionality reduction clustering, including epithelial cells, immune cells, and stromal cells in five cancer types (Fig. 6C, Fig. S6A). The high-expression group had a higher proportion of malignant tumor cells, while the low-expression group displayed greater T cell infiltration (Fig. 6D). Furthermore, the UMAP plot revealed elevated mean expression in epithelial cells, reinforcing the association between high expression of RMRs and tumor parenchymal cells (Fig. 6D).

Differential gene analysis and pathway enrichment demonstrated that these genes were involved not only in tumor proliferation pathways but also in immune cell differentiation and immune cytokine activity (Fig. 6E, Fig. S6B, Table S10). Cell communication analysis indicated the high-expression group exhibited more frequent and stronger interactions (Fig. 6F). Notably, the reciprocal molecules of macrophage migration inhibitory factor (MIF) and Midkine (MDK) were more pronounced in the high-expression group than in the low-expression group (Fig. 6G).

Based on our analysis, we concluded that high expression of RMRs exhibited characteristics resembling the basal-like subtype, including high proliferation, immune depletion, high malignancy, and poor prognosis. Conversely, low expression of RMRs correlated with increased infiltration of immune cells and large stromal cells, a lower degree of malignancy, and favorable survival outcomes.

Decoding the impact of RMR expression on tumor response and unveiling novel therapeutic avenues

Unveiling the impact of RMRs expression on tumor response, we collected multiple datasets of patients who had received immunotherapy or chemotherapy. Upon dividing them into high- and low-expression groups, a significant finding emerged: the high-expression group exhibited a higher proportion of non-responders (Fig. 7A-B). In addition, analysis of IC50 values for various drugs revealed that the high-expression group displayed insensitivity to conventional chemotherapeutic agents (Fig. 7C). We also used the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm [44] to predict immunotherapy response and found the high-expression group showed more immune exclusion scores than the low-expression group (Fig. 7D). Furthermore, in order to evaluate RMRs’ capacity to forecast resistance to chemotherapy and immunotherapy, we compared them with additional characteristics gathered from the published articles. Single-cell dataset of colorectal cancer includes cancer epithelial cells and various microenvironmental cells (Fig. S7A), from which we selected tumor cells to assess the predictive power of RMRs and other genes to predict drug response (Fig. S7B). RMRs were observed to score considerably higher in the non-respond group and to exhibit superior accuracy, AUC, precision, and F1 values when evaluating treatment response. (Fig. S7C-D). Since the scores of RMRs were higher in tumor cell clusters 2 and 4, we designated those clusters as RMRs + epithelial cells and the remaining of epithelial cells as RMRs-epithelial cells, which was in line with the expectation that RMRs + epithelial cells accounted for a higher percentage of non-response patients (Fig. S7E-G). To further explore the relationship between RMRs and immune cells, we did cellular communication between RMRs+/- epithelial cells and other microenvironmental cells. RMRs + epithelial cells were found to have more active MIF, MK, and APP signaling pathways with microenvironment cells than RMRs-epithelial cells. These pathways have been reported to be associated with tumor progression and immunosuppression (Fig. S7H). Similarly, we performed a validation analysis in ovarian cancer and found that RMRs in ovarian tumor epithelial cells were also significantly more potent in predicting treatment response than the genes in other studies (Fig. S7I-K). The RMRs scores of tumor epithelial clusters 1,2,3 were significantly higher than those of cluster 0, so we named clusters 1,2,3 as RMRs + epithelial clusters and the rest as RMRs-epithelial clusters (Fig. S7L-M). Consistent with colorectal cancer, the RMRs + epithelial cluster was concentrated in the non-response patient group (Fig. S7N). Cellular communication of RMRs+/- epithelial cells with other immune or stromal cells was performed, and it was evident that the MIF signaling pathway was more active in RMRs + epithelial cells (Fig. S7O-P). At the same time, we assessed the ability of RMRs and other genes to predict treatment effects at the level of pan-cancer cell lines and bulk samples. Notably, RMRs outperformed other genes in terms of AUC values for predicting the impact of response to immunotherapy and chemotherapy, indicating their superiority in determining therapeutic potential (Fig. S8A-N).

Fig. 7

The therapy response of RMRs to drugs. A. The fraction of patients with clinical response to immunotherapy. B. The fraction of patients with clinical response to chemotherapy. The difference was tested by the chi-square test. C. Volcano plot showing drugs that were differentially sensitive in the two groups. D. Boxplot revealing dysfunction and exclusion scores in different groups by TIDE methods, in which a higher score means ineffective immunotherapy. E. Correlation of RMRs expression with drugs. The table shows the top 10 positively and negatively correlated compounds from the connectivity map. The target score ranged from − 1 (negative connectivity) to + 1 (positive connectivity). F. GSEA of RMRs in various cancer types. The enrichment score > 0 demonstrated the positive correlation between RMRs expression and the activity of Hallmark pathways. G. Demonstration of JAK-STAT signaling pathways upstream and downstream genes. H. Chemical structural formula of JAK inhibitor TG-101,209

Then, we searched for potential compounds targeting RMRs expression in various cancer types by using connectivity map analysis (CMap) [38], which was employed to reveal functional links between small molecule compounds, genes, and disease states. Notably, the JAK inhibitors, TG-101,209 and TG-101,348 were effective compounds in negatively regulating RMRs expression, while the GSK inhibitor had the opposite effects (Fig. 7E, Table S11). Gene set enrichment analysis (GSEA) of RMRs also revealed positive enrichment in the JAK-STAT pathway (Fig. 7F). Among RMRs, the m5C writer NOP2 contributed the most to distinguishing RMRs clusters and was highly expressed in most cancer types and influenced their survival (Fig. 4F, Fig. S9A-B). We found high expression of NOP2 was related to tumor proliferation and JAK-STAT pathways and also responded to JAK inhibitors in multiple cancer types (Fig. S9C-D). The JAK/STAT pathway is an important cellular cascade that controls a wide range of processes, including cell differentiation, proliferation, and apoptosis [45] (Fig. 7G). The cytokine receptors activate Janus kinases, JAK inhibitors usually target these enzymes to intervene in tumor progression. We therefore hypothesized that TG-101,209, which is a selective JAK2 inhibitor could increase the therapeutic effect in the high expression group (Fig. 7H).

Enhancing chemo- and immuno-therapy efficacy in pancreatic cancer with TG-101,209

To gain deeper insights into the progression of pancreatic tumors, we conducted a comprehensive study involving 16 pancreatic cancer patients, we examined the impact of RMRs on tumor development (Fig. 8A, Table S12). Patients were divided into high and low RMR groups to investigate the DEG patterns related to cell cycle regulation, proliferation, and immune modulation (Fig. 8B). Our analysis identified 816 DEGs between the two groups, shedding light on the molecular alterations associated with RMR-mediated pancreatic tumor progression. Notably, the high-expression group showed a significant association with the JAK-STAT signaling pathway, suggesting the potential of using the JAK inhibitor TG-101,209 as a targeted therapy for these patients (Fig. 8C).

Fig. 8

TG-101,209 promoted pancreatic cancer patient response to chemotherapy and immunotherapy. (A) The heatmap depicting nine regulators expression of 16 PAAD patients (B) Volcano plot describing differential expression genes between high and low groups. Bar plots of KEGG and GO pathways of differential expression genes. (C) GSEA analysis of high-expression group and low-expression-group, the pathways were selected under P < 0.05. CMap results of predicted potential drugs of high expression of nine regulators. (D) Pancreatic tumor cells were subcutaneously transplanted into 4-6-week-old C57/BL mice. The mice were administered GEM, TG-101,209, or GEM combined TG-101,209. Box plot illustrating quantitative comparisons of tumor size across different treatment groups and control groups. (E) K-M survival curve of tumor-bearing mice during dosing. Death of mice as the endpoint of the survival record. The transplanted tumor cells were injected into the mouse pancreas in situ. (F) The tumors were stripped out from mice, and the frozen sections were stained with anti-Ki67. DAPI staining was included to visualize the nuclei. (G) Line graph telling proliferation ability in different treatment groups of KPC1 and KPC2 cell lines. (H) Heat map displaying the interaction of GEM and TG-101,209. (I) Mice were subcutaneously injected with pancreatic tumor cells. Tumor-bearing mice were injected with anti-PD1, TG-101,209, or a combination of anti-PD1 and TG-101,209. Tumor growth curve showing growth rate and volume size until mice sacrificed. (J) Survival curves of orthotopic pancreatic tumor-bearing mice in control, anti-PD1, TG-101,209, and combined therapy groups. The time at which the mice were sacrificed was the time at which the survival curve was finally recorded. K. Flow cytometry sorting CD8 + T cell types and the mathematical statistics were made in CD8 + T cells, effector CD8 + T cells, and exhausted CD8 + T cells of different immunotherapy groups

Subsequently, our in vivo and in vitro experiments demonstrated a significant reduction in tumor volume when TG-101,209 was administered in combination with GEM (Gemcitabine), a conventional chemotherapy drug. Moreover, the combination therapy resulted in improved overall survival compared to the use of either drug alone (Fig. 8D-E). Complementing these findings, mIHC analysis revealed a notable decrease in Ki67 levels, a well-established marker of cellular proliferation, within the combination chemotherapy group (Fig. 8F). Additionally, CCK-8 assays demonstrated that co-medication effectively suppressed tumor cell vitality, particularly in pancreatic cancer cell lines (Fig. 8G). Notably, our analysis revealed a strong synergistic effect with a ZIP Synergy score of 33.619 between GEM and TG-101,209, further supporting the potential clinical utility of this combination therapy (Fig. 8H).

We also investigated the sensitizing effect of TG-101,209 in combination with immunotherapy. Monotherapy using PD1 blockade alone displayed limited efficacy in tumor elimination; however, when combined with TG-101,209, it exhibited enhanced effectiveness (Fig. 8I). Co-administration of TG-101,209 and PD1 blockade significantly extended the survival of mice bearing experimentally induced tumors (Fig. 8J). Furthermore, flow cytometric analysis revealed that combination immunotherapy augmented the number of CD8+T cells and effector CD8+T cells secreting cytotoxic proteins while reducing the population of exhausted CD8+PD1+T cells (Fig. 8K).

In conclusion, these findings provide valuable insights into the treatment strategies for pancreatic cancer, suggesting the potential utility of targeting RMRs as a therapeutic approach to optimize patient outcomes.

TG-101,209 sensitized high-RMRs-expressing tumors to chemotherapy across cancer types

We further investigated the antitumor activities of TG-101,209 and other chemotherapeutic drugs in various cancer cell lines with high-expression RMRs, including colorectal, liver, breast, and lung cancers (Fig. 9A). Co-administration of TG-101,209 with L-OHP (Oxaliplatin) or CAPE (Caffeic acid phenethyl ester) showed the most effective tumor regression in colorectal cancer compared to mono- or dual-drug treatments (Fig. 9A). Similarly, in liver, breast, and lung cancers, combinations of TG-101,209 with 5-FU (5-Fluorouracil), L-OHP, ADM (Adriamycin), CTX (Cyclophosphamide), GEM, or DDP (Cisplatin) resulted in smaller tumor sizes compared to mono-chemotherapy groups (Fig. 9A).

Fig. 9

TG-101,209 increased tumor response to chemotherapeutic agents. (A) Subcutaneous growth of mice in the control group, single-drug, two-drug combination, and TG-101,209 in combination with two traditional drugs groups (n = 6). (B) The volume change of nude mice implanted tumors was recorded every 3 days. Each bar represents the mean ± SD for six animal measurements. (C) The K-M survival analysis in different chemotherapy drug treatment groups was performed until the sacrifice of mice. Tumor cell lines were transplanted to the liver in situ of mice. (D) CCK-8 assay was performed to determine the cell viability after treatment in MCF-7, A-549, SW480, and Hep-G2 cell lines. (E) Representative ki67 and 4-6‐Diamidino‐2‐phenylindole (DAPI) immunofluorescence staining of tumor sections from the tumor tissues in different treatment groups. (F) Three-dimensional stereo thermograms showing synergistic promotion between two-by-two drugs. * P < 0.05, **P < 0.01 or *** P < 0.001

Combination therapy significantly improved the survival of tumor-bearing mice and suppressed tumor growth during treatment (Fig. 9B-C). In vitro analysis using CCK-8 assay demonstrated that combined therapy effectively disrupted tumor cell viability (Fig. 9D). Additionally, DAPI and Ki67 immunofluorescence staining of tumors implanted in mice showed a significant reduction in cell proliferative activity with the combination therapy (Fig. 9E). Evaluation of combinatory effects using SynergyFinder revealed strong synergistic effects between traditional chemotherapeutic agents and TG-101,209 (Fig. 9F).

We examined the protein levels of JAK signaling pathway genes in the TG-101,209-treated and control groups and the activated JAK1/JAK2/STAT1/STAT3 proteins levels were lower in the treated cell lines, suggesting that TG-101,209 specifically inhibits the proteins of the JAK-STAT signaling pathway (Fig. S10A). By comparing the expression levels of RMRs genes in the blank control group and the TG-101,209 treatment group, we found that in the absence of other interfering agents, TG-101,209 significantly decreased the mRNA levels of RMRs genes. However, if the drug HJC0146, an inhibitor of the total JAK-STAT pathway was added in advance, the effect of TG-101,209 in decreasing the expression of RMRs was not obvious. Therefore, TG-101,209 functioned as a specific inhibitor of the JAK-STAT signaling pathway thereby reducing the transcriptome level of RMRs (Fig. S10B-J).

TG-101,209 sensitized high-RMRs-expressing tumors to Immunotherapy across cancer types

Furthermore, we observed that the combination of anti-PD1 and TG-101,209 had a greater impact on reducing tumor volume and prolonging the survival of mice with tumors compared to anti-PD1 or TG-101,209 alone (Fig. 10A-C). The combination therapy also resulted in greater inhibition of tumor cell viability (Fig. 10D). Flow cytometry analysis of harvested tumors indicated reduced infiltration of CD8+PD1+T cells but improved ratios of CD8+T cells and increased functional factors (IFNγ, TNFα, GZMB, and perforin) within CD8+T cells (Fig. 10E). We evaluated the magnitude of mutual reinforcement between anti-PD1 and TG-101,209, showing strong ZIP synergy scores (Fig. 10F).

Fig. 10

TG-101,209 enhanced tumor response to immunotherapy. (A) The tumor size and volume were documented in anti-PD1, TG-101,209, two-drug combination groups, and control groups. (B) The tumor growth curve was recorded every three days during the period of treatment. Each bar represents the mean ± SD for six animal measurements. (C) The survival time curve of tumor-bearing mice in different treatment groups. Death of mice as the endpoint of the survival record. Tumor cell lines were injected into the orthotopic liver of mice. (D) Cell viability curve of different immunotherapy groups. (E) Flow cytometry analysis was performed in the harvested tumors; Representative dot plots and bar plots of the percentage of CD8 + T cells, CD8 + TNFα + T cells, CD8 + IFNγ + T cells, CD8 + PD1 + T cells, CD8 + perforin + T cells and CD8 + Grzmb + T cells. (F) The heat map shows the synergistic effect of the anti-PD1 and TG-101,209. ZIP Synergy scores > 10 indicate synergism (red regions), scores < -10 indicate antagonism (green regions), and scores between − 10 and 10 mean the interaction between two drugs was likely to be additive. * P < 0.05, **P < 0.01 or *** P < 0.001

By evaluating the contribution of immune cells, such as effector CD8 + T cells, B cells, NK cells, regulatory T cells (Tregs), and myelosuppressive cells, to immune activation and clinical prognosis in tumor patients in the TCGA pan-cancer dataset, we were able to derive coefficients for each immune cell. This allowed us to further evaluate the relationship between TG-101,209 and immune microenvironment activation. (Fig. S11A). We further combined the weighted model coefficients of the number of infiltrations of each immune cell in the experimental mouse model to evaluate the TME-activated scores in the TG-101,209 group and the control group. In multiple cancers, the high immune score group had a higher proportion of the group treated with TG-101,209, and the immune activation function of the treatment group was also higher than the control group. The TME-activated score can well predict the TG-101,209 treatment group and control group with high AUC of the area under the ROC curve (Fig. S11B-F).

Overall, the combination therapy with TG-101,209 enhanced the killing function of CD8+T cells, reduced exhausted CD8+T cells, and decreased tumor burden in high-RMRs-expressing tumors. These findings highlight the potential of combined TG-101,209 therapy as an effective approach for improving treatment outcomes in various cancer types.