The batch-corrected, normalized pan-cancer and normal tissue datasets were obtained from the University of California, Santa Cruz (UCSC) datasets, including The Cancer Genome Atlas (TCGA; containing 33 cancer types), Therapeutically Applicable Research to Generate Effective Treatments (TARGET; containing 7 pediatric cancers), and the Genotype-Tissue Expression (GTEx), composed of 54 normal tissues). Cancer types with less than 3 samples were excluded. The simple-nucleotide variation (SNV) data processed with Mu Tect2 software (22) were downloaded from the GDC data portal (https://portal.gdc.cancer.gov/).

The single-cell sequencing datasets of bladder urothelial carcinoma (BLCA) (GSE145137), cholangiocarcinoma (CHOL) (GSE125449), glioblastoma multiforme (GBM) (GSE138794), head and neck squamous cell carcinoma (HNSC) (GSE103322), kidney renal clear cell carcinoma (KIRC) (GSE121636 and GSE171306), liver hepatocellular carcinoma (LIHC) (GSE125449), and ovarian serous cystadenocarcinoma (OV) (GSE118828) were collected from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), and lung adenocarcinoma (LUAD) single-cell dataset was downloaded from the National Center for Biotechnology Information (NCBI) BioProject #PRJNA591860. The FASTA sequences of JMJD8 protein were obtained from the NCBI protein database. The molecule structure data were downloaded from the free dataset ZINC (https://zinc15.docking.org/).

The JMJD8-associated diseases or phenotypes were first presented with a bubble graph from the Open Target Platform (https://platform.opentargets.org/). The JMJD8 mRNA expression differences between cancer tissues and adjacent normal tissues were analyzed in the Gene_DE module of TIMER2.0 (23) (http://timer.cistrome.org/). For the cancer types that lack enough normal adjacent samples, the GTEx and TCGA data were adopted and analyzed by the Box Plot module in the Expression DIY function of GEPIA2.0 (5) (http://gepia2.cancer-pku.cn/); the absolute log2FC cutoff and p-value cutoff were set as 0.585 and 0.01, respectively. Moreover, mRNA expression differences in cancer stages were also presented in the Stage Plot module. The Clinical Proteomic Tumor Analysis Consortium (UPTAC) data, obtained from UALCAN (1) (http://ualcan.path.uab.edu/), were used to compare the JMJD8 protein expression differences between cancer and normal tissues.

The package “survival” was run to analyze the pan-cancer prognostic risks of JMJD8, and the “survfit” function was used to analyze the survival differences and visualized by the Kaplan–Meier curves.

Three genomic alteration types (mutation, amplification, and deep deletion) frequency in pan-cancer was analyzed via the Cancer Types Summary module of the online web tool cBioPortal (24) (https://www.cbioportal.org/). The processed pan-cancer SNV data were collected and combined with the protein domains retrieved from the R package “maftools” to show the mutation landscape of JMJD8. The Kaplan–Meier curves were obtained from the Copy_Number module of Tumor Immune Dysfunction and Exclusion (TIDE) (25), a web tool to study the correlations between cancer genomic or transcriptional alterations and immunotherapy responses, to explore the prognostic significance of JMJD8 copy number variations (CNVs). Moreover, the package “maftools” was used to calculate the tumor mutation burden (TMB), the microsatellite instability (MSI), homologous recombination deficiency (HRD), and neoantigen data obtained from previous studies (26, 27); the correlations between TMB, MSI, neoantigen, and the expression of JMJD8 were then presented.

The association between 5 mismatching repair genes (28), 4 DNA methyltransferase (29), and JMJD8 expression in pan-cancer was visualized. To further analyze the HRR signature in pan-cancer, 30 HRR-related genes were retrieved from an ARIEL3 clinical trial (30) and input into GEPIA2.0 to calculate their correlations with JMJD8. The “differentially methylated probes-based stemness index” (DMPsi) of each cancer type was also obtained from a previous study (31), and their intercorrelation with JMJD8 mRNA expression was investigated. The association between JMJD8 promoter methylation and cytotoxic T lymphocytes (CTLs) and patients’ survival rates were subsequently presented in the Methylation module of TIDE. Finally, a heatmap was plotted to display the correlation between JMJD8 and 44 N1‐methyladenosine (m1A), 5‐methylcytosine (m5C), and N6‐methyladenosine (m6A) modification genes (32–34).

To seek the clinically relevant alternative splicing (AS), the ClinicalAS module of the OncoSplicing server (4) (http://www.oncosplicing.com/) was searched for the AS events of JMJD8 included in both the SplAdder and the SpliceSeq projects. The PanPlot was displayed to show the percent spliced-in (PSI) of TCGA cancers and GTEx tissues. The PanDiff plots were presented to compare the PSI differences of the AS events (detected in more than 3 cancers) between cancers and adjacent or GTEx normal tissues. Finally, the Kaplan–Meier curves were plotted to explore the prognostic significance of the AS events in pan-cancer.

The current JMJD8 Protein–Protein Interaction network with known experiment validations was searched via the online web tool String (35) (https://www.string-db.org/). The pan-cancer pathway-level somatic alterations of several vital pathways were explored on UALCAN, and the expression correlations between pathway-related signature (36) and JMJD8 were investigated on GEPIA2.0. Also, the top 100 co-expression genes of JMJD8 in pan-cancer were obtained from the Similar Gene Detection function of GEPIA2.0. The top 100 genes were adopted for functional analysis of Gene Ontology (GO) with a false discovery rate (FDR) < 0.05 using the R package “clusterProfiler” the GO annotation was retrieved by R package “org.Hs.eg.db”. For the top 5 co-expression genes, JMJD8’s correlations with them were visualized by the heatmap plotted in TIMER2.0 and the scatter plots in GEPIA2.0. The quantification results of JMJD8 functional enrichment were further investigated by Gene Set Enrichment Analysis (GSEA) using the GSEA software (http://software.broadinstitute.org/gsea/index.jsp) (37). According to the median expression of JMJD8, all cancer samples were divided into low- and high-JMJD8 groups, and the gene sets c2.cp.kegg.v7.4.symbols.gmt and h.all.v7.4.symbols.gmt were downloaded from the Molecular Signatures Database (38) of http://www.gsea-msigdb.org/gsea/downloads.jsp for Kyoto Encyclopedia of Genes and Genomes (KEGG) and hallmark pathways GSEAs (p < 0.05).

The roles of JMJD8 in pan-cancer microenvironment infiltration were first investigated by calculating the Estimation of STromal and Immune cells in MAlignant Tumour tissues using Expression data (ESTIMATE), stromal, and immune scores using the R package “ESTIMATE” (version 1.0.13) (39).

The immune checkpoint markers from a previous study (27) were obtained to analyze their correlations with JMJD8, JMJD8’s expression relevance with immune subtypes was analyzed, and its levels were compared between these subtypes in pan-cancer in the Subtype module of TISDB, a web portal for tumor and immune system interaction (3) (http://cis.hku.hk/TISIDB/). Additionally, the heatmaps were plotted to show the association between JMJD8 expression and chemokines, chemokine receptors, and immunostimulators in the Chemokine and Immunomodulator modules, respectively. To investigate cytokine treatment’s effects on JMJD8 expression, the Tumor Immune Syngeneic MOuse (40) (TISMO, http://tismo.cistrome.org/) web tool was used to compare gene expression levels across cell lines between pre- and post-cytokine-treated samples.

TIMER2.0 was used to run 5 different algorithms on the intercorrelation between M2 macrophage and JMJD8. Moreover, the SpatialDB online tool (41), a database for spatially resolved transcriptomes (https://www.spatialomics.org/SpatialDB/), was applied to analyze the spatial expression levels and overlapping of JMJD8 and M2 macrophage markers CD68 and CD163 in breast cancer and prostate cancer. For single-cell resolution, the JMJD8 expression among various cell subtypes in pan-cancer was compared with the aid of single-cell datasets collected from GEO and the Tumor Immune Single-cell Hub (TISCH) (42) (http://tisch.comp-genomics.org/). Dataset integration for KIRC was performed by the RunHarmony function from the R package "harmony". Subsequently, the immunocyte infiltrating correlations of JMJD8 were calculated by the CIBERSORT algorithm (43).

Multiple fluorescence staining was performed on pan-cancer paraffin sections to validate the M2 macrophage biomarker potential of JMJD8. The sections in this study contained 9 cancer types. These sections were deparaffinized and blocked with 3% H₂O₂ and 2% bovine serum albumin (BSA) after antigen retrieval. Subsequently, they were incubated sequentially with the three primary antibodies, JMJD8 (mouse, 1:100, Santa Cruz, Dallas, TX, USA), CD68 (rabbit, 1:3,000, Servicebio, Wuhan, China), CD163 (rabbit, 1:3,000, ProteinTech, Wuhan, China). After primary antibody labeling, the sections were incubated by horseradish peroxidase (HRP)-conjugated secondary antibody (GB23301, GB23303, Servicebio, China), followed by tyramide signal amplification (TSA) (fluorescein isothiocyanate (FITC)-TSA, CY3-TSA, and CY5-TSA) (Servicebio, China). Next, 4′,6-diamidino-2-phenylindole dihydrochloride (DAPI) counterstaining of nuclei, the antifade mounting medium was applied, and the Pannoramic Scanner was used (3DHISTECH, Budapest, Hungary) to obtain multispectral images of the stained sections.

For fluorescence spectra, the excitation wavelength and emission wavelength for different fluorescence dyes are listed respectively as follows: DAPI (blue, 330–380 and 420 nm), CY3 (red, 510–560 and 590 nm), CY5 (pink, 608–648 and 672–712 nm), and FITC (green, 465–495 and 515–555 nm). Caseviewer (C.V 2.4) and Pannoramic viewer (P.V 1.15.3) image analysis software were used to quantify the cells with positive staining at single-cell levels in the multispectral images.

TIMER2.0 was used to present the JMJD8 expression intercorrelations with regulatory T cells (Tregs), cancer-associated fibroblasts (CAFs), myeloid-derived suppressor cells (MDSCs), and CD8+ T cells via different immune algorithms. Also, JMJD8’s roles in T-cell dysfunction and CTL-related prognosis in pan-cancer subtypes were investigated on TIDE.

To investigate the effects of JMJD8 on the routine therapeutic response for glioblastomas and breast cancers, the JMJD8 expression differences between responders and non-responders and the receiver operating characteristic curves (ROC) of therapy-related survivals were searched on ROC plotter (44) (www.rocplot.org), an online tool to link gene expression and response to therapy using transcriptome-level data of four cancer types. The “query” tool of cMap (45) (https://clue.io/) was applied to screen out the anti-JMJD8 chemical compounds, the heatmap was plotted to show the top 30 compounds against the JMJD8-related differentially expressed gene signature, and their mechanisms of action (MoA) were also displayed. Moreover, JMJD8’s correlations with drug sensitivities were demonstrated on RNAactDrug (46) (http://bio-bigdata.hrbmu.edu.cn/RNAactDrug/), a comprehensive resource for querying associations between drug sensitivity and RNA molecules, and drugs with FDR < 0.05 were selected. The JMJD8 expression and the concentrations that cause 50% growth inhibition (GI50) of the top four cMap compounds among cell lines were analyzed using the COMPARE tool (https://nci60.cancer.gov/publiccompare/) in the Developmental Therapeutics Program (DTP) of the US National Cancer Institute (NCI).

For protein–compound interactions, homology modeling of JMJD8 protein was performed with the AlphaFold2 software (47). The rank_1 unrelaxed protein structure was estimated on SAVES v6.0 (https://saves.mbi.ucla.edu/) and applied to molecular docking. The docking was performed with the Discovery Studio software (version 4.5). After Auto Preparation of JMJD8 and ligand preparation of compounds, the binding sites, and all compound conformations were identified, the LibDock was selected for docking. The site and molecule conformation with the highest LibDockScore were determined for final interaction. The binding pocket 3D view and the intermolecular forces distance 2D view were displayed.

All the bioinformatics analyses were conducted on the R software. A log-rank test assessed the survival significance. The correlations were quantified by Pearson's or Spearman's correlation coefficients. The correlations between JMJD8 and methyltransferases were considered positive/negative if the statistical significance of the correlation was observed in any one of the 4 methyltransferases, and no opposite statistical result in other methyltransferases was presented. p-Values less than 0.05 (*p < 0.05) were considered significant.

The JMJD8-related disease was explored on OpenTarget, and the bubble graph shows that JMJD8 was associated with non-small lung carcinomas (Figure 2A). We then analyzed the JMJD8 mRNA level differences between pan-cancer and corresponding normal tissues on TIMER2.0; JMJD8 mRNA was significantly upregulated in 11 cancer types (BRCA, CHOL, COAD, ESCA, GBM, HNSC, KIRP, LIHC, LUAD, PRAD, and STAD) and downregulated in CESC, KICH, KIRC, and THCA (Figure 2B). For cancers lacking normal tissues, we compared their JMJD8 expression differences on GEPIA2.0 and UALCAN. The GEPIA2.0 results show that JMJD8 mRNA in cancers was highly expressed in DLBC, THYM, LGG, and PAAD but poorly expressed in TGCT and UCS (Figure 2C). Additionally, JMJD8 was negatively correlated with the high stage of ESCA, PAAD, and THCA (Figure 2D). The UALCAN results exhibited that JMJD8 protein was upregulated in cancers of BRCA, UCEC, and GBM and downregulated in LIHC and HNSC (Figure 2E). For the prognostic significance of JMJD8, forest plots of the risks (Supplementary Material S1) and the Kaplan–Meier curves were plotted on TCGA data. We noticed that high expression of JMJD8 was associated with lower overall survival (OS) percentages in GBM and LGG; lower disease-specific survival (DSS) percentages of GBM, LGG, and STAD; and shortening of the progression-free interval (PFI) of ACC, LGG, and STAD. By contrast, it correlated with higher OS in ESCA, PCPG, THYM, PRAD, and SARC and higher DSS in PCPG and elongation of PFI in SARC and THYM (Figure 2F). These results indicated that JMJD8 might be a cancer driver gene in gliomas and promote ACC progress but might be a protective gene in THYM.

Cancers characterized genomic alterations. To seek whether JMJD8 gene is altered at the genome level, we displayed the JMJD8 pan-cancer CNV and SNV analysis results and found high JMJD8 amplification in BRCA and high deep deletion rates in diffuse large B-cell lymphoma and UCS (>3%), while no high SNV rates were observed (Figures 3A, B). When the CNV levels were applied for patient grouping on TIDE, high JMJD8 CNV group patients showed higher survival rates in AML, KIRC, COADREAD, and LIHC but lower survival rates in UCEC, BRCA (HER2), HNSC (HPV+), and PADD (Figure 3C). Moreover, we compared the TMB, MSI, neoantigens, and ploidy correlations with JMJD8 in pan-cancer since these genomic alterations frequently appeared in cancers and affected patient prognosis and therapeutic responses (48–50). As displayed in Figure 3D, JMJD8 was positively correlated with TMB in 2 cancers (LGG and UCEC) and with MSI in 7 cancers (COAD, KICH, KIRC, LIHC, LUSC, TGCT, and UCEC). By contrast, it was negatively correlated with TMB in 6 cancers (BRCA, CESC, LUAD, PCPG, PRAD, and THCA) and with MSI only in BRCA and SARC. As for HRD, JMJD8 showed a negative correlation with it in BRCA (correlation coefficient is nearly −0.4), followed by SARC, LUAD, and BLCA, while their positive correlations were observed in THCA and HNSC. For aneuploidy, JMJD8 was positively correlated with it in 3 cancers (ESCA, HNSC, and UVM) and negatively associated with 7 cancers (TGCT, UCEC, KICH, SARC, KIRC, KIPAN, and BRCA) (Figure 3E). Cancer TMB and MSI often caused neoantigen presentation. As presented in Figure 3F, 5 cancers (BRCA, LUAD, PRAD, CESC, and STAD) showed negative associations between neoantigens and JMJD8 expression, and only KIRP presented a positive correlation. The significance threshold for all analyses was set at p < 0.05. All the results above strongly suggested that JMJD8 is a potential biomarker of genome stability in BRCA and LUAD.

Cancer genomic stability relied mainly on repairing the DNAs via different mechanisms, including DNA MMR (28) and HRR, which also contributed to stemness maintenance in cancers (12, 15). Hence, we analyzed the correlations between JMJD8 and MMR-related genes (EPCAM, MLH1, MSH2, MSH6, and PMS2), HRR signature, and cancer stemness.

We discovered that JMJD8 was positively correlated with multiple MMR genes in most cancers, including ACC, CESC, GBM, HNSC, KIRC, KIRP, PAAD, LIHC, PCPG, STAD, and especially THCA (Figure 4A). For cancer stemness, we noticed that JMJD8 obtained strong correlations with it in OV, followed by LGG, UVM, HNSC, and ESCA. In TGCT, KIRP, and KIPAN, they exhibited negative correlations (Figure 4B). Unsurprisingly, 9 cancers, which showed positive correlations between stemness and JMJD8, also presented consistent trends when it came to HRR signature (Figure 4C), demonstrating that JMJD8 interplayed with DNA repair-mediated cancer stemness.

Since the DMPsi reflexed the DNA methylation status of cancer, we subsequently sought the JMJD8’s influences on cancer epigenetic modulations. As depicted in Figure 4D, JMJD8 had negative correlations with methyltransferases significantly in BRCA and other 4 cancers (MESO, PRAD, SKCM, THYM). Inversely, their positive association was exhibited in 19 cancers, including LGG, LIHC, LUSC, OV, PAAD, PCPG, STAD, TGCT, THCA, UCEC, ACC, BLCA, CESC, ESCA, GBM, HNSC, KICH, KIRC, and KIRP. Their strong positive correlations in cancers suggest that high JMJD8 promoted the promoter methylation of its target genes and suppressed their expression.

We also searched TIDE for the interplays between JMJD8 promoter methylation and cancer subtypes, CTL, and risks. Figure 4E presents the list of the top 8 CTL-correlated cancer subtypes; JMJD8 promoter methylation was positively correlated with CTL infiltration in BRCA, PAAD, and MESO. Figure 4F exhibits the correlation scatter plots and the Kaplan–Meier curves, demonstrating that JMJD8 promoter methylation was associated with CTL infiltration and predicted more prolonged survival in three BRCA subtypes and MESO.

In addition to DNA methylation, we also investigated the correlation between JMJD8 and RNA modulator gene expression. Surprisingly, we discovered that high JMJD8 was associated with a majority of RNA modulator genes in many cancers, including m1A, m5C, and m6A (Figure 4G), indicating that JMJD8 was involved in RNA modifications.

AS is a common post-transcriptional modification type, producing various transcripts and subsequent proteins or non-coding RNAs. Its dysregulation frequently occurs in cancers and affects tumorigenesis (51). We analyzed the ASs on OncoSplicing, 5 clinical-relevant AS events were identified, we mainly displayed the Intron_Retention_51257 event here, and the other 4 events are presented in Supplementary Material S2. Figures 5A, B exhibit the splicing mode and the PSI of Intron_Retention_51257 in pan-cancer; cancers such as LUSC and READ showed higher PSI than the normal samples. Figure 5C summarizes the statistical results of the PSI differences between tumors and normal/adjacent tissues, and those with prognostic values were presented in Figure 5D by the Kaplan–Meier curves. High PSI predicted lower OS and DSS in KIRC and lower OS in MESO. Similarly, high PSI also predicted lower disease-free interval (DFI) and PFI in both LIHC and PRAD. However, in PAAD, high PFI was associated with longer DFI and PFI. These results implied the biological importance of regulated JMJD8 As events in cancer progression.

To investigate the functional roles of JMJD8 in cancers and the interactive or co-expressed genes, functional enrichment analyses were sequentially conducted. The interactive proteins with experimental validations were obtained from String web tool, and 10 proteins were displayed (Figure 6A). We then compared the JMJD8 expression between altered and non-altered pathways on UALCAN and noticed that JMJD8 was elevated in altered SWI/SNF complex, p53/Rb-related pathway, and chromatin modifiers status in HNSC and GBM while poorly expressed in BRCA (Figure 6B). Since somatic mutations or CNVs cannot directly explain the expression level, we explored the expression correlation between JMJD8 and these pathway-related signatures (36). We found that these signatures positively correlated to JMJD8 (Supplementary Material S3). The top 100 JMJD8 coexpressed genes in pan-cancer were analyzed on GEPIA2.0, and the top 5 genes (C16ORF58, IFT140, ITFG3, PIGQ, and WDR24) showed high correlations with JMJD8 in the majority of cancer types (Figure 6C). The functional enrichment of GO terms exhibited multiply cellular skeleton and ciliary transportation system-related activities (Figure 6D). Additionally, the GSEA results of GO and KEGG suggested the close association between metabolism, immune activities, and JMJD8 (Figure 6E).

To investigate the immunological roles of JMJD8 in the cancer environment, we calculated the ESTIMATE of JMJD8 in pan-cancer. As depicted in Figure 7A, JMJD8 was reversely correlated to ESTIMATEScore and ImmuneScore in many cancers, including TCGA tumors THCA, KIPAN, MESO, ACC, GBM, BRCA, CESC, THYM, and TARGET-WT of kidney tumors. However, JMJD8 was also positively associated with them in several cancers like UVM. It presented positive relevance with JMJD8 in UVM and TGCT et al for stromal cell infiltration and negative relevance in KIPAN, MESO, THCA, KIRP, and CHOL et al. The cancers showing negative JMJD8–ImmuneScore correlations also harbored negative correlations with most immune checkpoint genes, including THYM, THCA, TGCT, BRCA, LUAD, and MESO (Figure 7B). However, we noticed that several markers were positively correlated with JMJD8 in many cancers. CD276 showed the highest positive correlations with JMJD8 in 15 cancers, followed by LGALS9, VSIR, and TNFRSF4. We also noticed that in LGG, JMJD8 was positively associated with as many as 24 immune checkpoint genes, indicating the involvement of JMJD8-related immune checkpoint effects.

Subsequently, we explored whether JMJD8 was differentially expressed in diverse cancer immune subtypes via TISIDB. The histogram exhibits that JMJD8 was significantly associated with immune subtypes in 10 cancers (Figure 7C), the top 6 of which are presented in Figure 7D; JMJD8 expression increased in the C4 subtype in BLCA, UCEC, and LUAD, implying its reverse association with lymphocyte functions. Moreover, we analyzed the associations between JMJD8 and chemokines, receptors, and immunostimulators. As visualized in heatmaps (Figure 7E), JMJD8 was negatively associated with several chemokines (CXCL9, 10, 11, 12, and 13), many receptors, and immunostimulators in pan-cancer. We also noticed that high JMJD8 promoter methylation was positively correlated with most immunostimulators, demonstrating that JMJD8 expression affected chemokine-mediated immunostimulations against cancers.

Finally, we compared the JMJD8 expression differences between pre- and post-cytokine treatment in cancer cell lines on web tool TISMO (Figure 7F). We discovered that the JMJD8 expression decreased after the IFN-γ treatment in four cell lines, and it also decreased in one IFN-β and one TNF-α posttreatment cell line.

The results from multiple perspectives demonstrated that JMJD8 is a critical factor in immunosuppressive environment construction in many cancers, probably via suppressing immunostimulator function and immune checkpoint effects.

To further investigate the JMJD8 in cancer immune, we analyzed its expression in levels of immunocytes. We first ran the CIBERSORT algorithm to obtain 22 immunocyte correlations with JMJD8. We discovered that JMJD8 presented strong positive correlations with M2 macrophages in TGCT, BRCA, and LGG and negative correlations with M1 macrophages and activated CD4+ memory cells in many cancers. Also, Tregs were positively associated with JMJD8 in 13 cancers (Figure 8A). Focusing on M2 macrophages, we conducted multiple algorithms on TIMER2.0 to analyze the correlation between their infiltration level and JMJD8 expression in pan-cancer, and the association was observed in BLCA, HNSC, STAD, TGCT, UCEC, and UVM consistently presented. (Figure 8B). Spatial transcriptional data on SpatialDB were obtained to depict the spatial overlapping of JMJD8 and M2 macrophage biomarkers CD68 and CD163 on BRCA and PRCA cancer tissues (Figure 8C), and as expected, JMJD8, CD68, and CD163 markers presented similar spatial distributions, which implied potential co-expression of JMJD8, CD68, and CD163. What is more, we retrieved JMJD8 expression data in single-cell cellular subtypes from TISCH. As exhibited in Figure 8D, JMJD8 was expressed by M2 macrophages or malignant cells in AEL, BRCA, glioma, HNSC, LIHC, NSCLC, and OV. We also compared JMJD8 expression in the collected cancer single-cell datasets and discovered that M2 or undefined macrophages expressed it in BLCA, CHOL, GBM, HNSC, and LIHC (Figure 8E).

For experimental validation, we performed multiple fluorescence (H&E staining images are presented in Supplementary Material S4) staining in pan-cancer paraffin sections, and the photos of JMJD8, CD68, and CD163 staining showed their co-expression in 7 cancers (BLCA, PUC/PRUC, ureter urothelial cell carcinoma (UCC), UCEC, LGG, TGCT, and PRAD) (Figures 9A–I). Moreover, the JMJD8 fluorescence intensity seemed higher in cancer tissues in THCA and TGCT adjacent tissues (Figures 9G, H), consistent with the differential expression results from TIMER2 and GEPIA2 Figures 2B, C. Several cancers showed positive correlations between JMJD8 and CD163 intensity, such as LGG.

From the bulk, spatial, single-cell transcriptional data and the fluorescence staining results above, we confirm the close association between JMJD8 and M2 macrophages, and these suggested that JMJD8 is a potential cancer-specific marker.

Apart from M2 macrophages, we also sought other members contributing to immunosuppression. We used TIMER2.0 to display the correlations between JMJD8 and Tregs, CAFs, and MDSCs, JMJD8's positive correlations with at least two cell types were observed in several cancers, including CESC, COAD, HNSC, LUSC, PAAD, STAD, and THYM (Figure 10A). Their purity-adjusted correlations are presented in Figure 10B, showing CAFs with the highest correlations with JMJD8 in CESC, COAD, HNSC, LUSC, PAAD, and THYM.

Since CTLs were the main affected cells during immunosuppression, we also investigated JMJD8’s association with CD8+ T cells via 10 algorithms on TIMER2.0, and considerable negative relevance was observed in BRCA, HNSC, LUAD, and SKCM (Figure 10C). In addition, T-cell dysfunction was found related to JMJD8 in neuroblastoma and BRCA via TIDE web tool and a negative correlation was discovered between CTL and JMJD8 in BRCA (Figure 10D).

The extra exploration of the association between JMJD8 and other immunosuppressive cells in pan-cancer indicated that JMJD8 also functioned in CAFs, Tregs, and MDSCs and inhibited the anticancer immune by targeting CTLs.

To seek whether JMJD8 can predict therapeutic responses to cancers, we obtained the data from ROCplotter to show the association between the therapeutic outcomes and JMJD8 expression in four cancer types (BRCA, OV, GBM, and CRC). In GBM, JMJD8 was highly expressed in non-responders after chemotherapy, especially post-nitrosourea treatment, and the area under the curve (AUC) value of post-nitrosourea 16-month OS reached 0.7. However, in BRCA, responders of post-anti-HER2, chemotherapy, and endocrine therapy harbored higher JMJD8 expression, with the highest AUC value of 5-year PFS prediction of anti-HER2 therapy reaching 0.909 (Figure 11A).

Given the poor therapeutic effects of high-JMJD8 GBM patients receiving routine chemotherapy, we attempted to identify potential anti-JMJD8 drugs with higher effects, which may also improve GBM sensitivity to current chemotherapy. The cMap tool was utilized to filter the compounds, causing opposite transcriptional alterations to those that increased by high-JMJD8 expression, in 9 different tumor cell lines and the top 30 compounds with JMJD8-targeted potential were displayed (Figure 11B). Notably, the MoA of six compounds was the pattern of histone deacetylase (HDAC) inhibitors, suggesting the mechanism of JMJD8 functions in cancers (Figure 11C). We also searched the RNAactDrug for drugs correlated to JMJD8 mRNA expression. The top 26 drugs with FDR < 0.05 were shown, while only 5 drugs were identified opposite to JMJD8 expression (Figure 11D).

For the top cMap compounds, we compared their GI50 after treating pan-cancer cell lines using COMPARE tools. For no testing data of XMD-1150, XMD-892, genipin was found; we exhibited the GI50s of THM-I-94 and compared them with the JMJD8 expression in pan-cancer cell lines. The average −log10(GI50) of THM-I-94 was −6.65, and in central nervous system (CNS) tumor cell lines, high JMJD8 expression corresponded to higher GI50 (Figure 11E). To further investigate whether these compounds can bind to JMJD8 protein, we performed homology modeling of JMJD8 protein and its molecular docking with potential drugs. With the use of alphaFold2.0, 5 models were built using the FASTA sequence listed in Supplementary Material S5, and the rank_1 model was retained and estimated with an Overall Quality Factor of 75.7322. Subsequently, molecular docking was conducted via Discovery Studio (version 4.5); XMD-1150 and XMD-892 failed to dock with JMJD8, while genipin and THM-I-94 succeeded, with the highest LibDockScore of 131.25 and 104.947, respectively. The 3D structures of the docking pockets and the 2D graph showing the interactive forces and distances are presented in Figure 11F. THM-I-94 interacted with JMJD8 residues by many types of interactions.

Taken together, THM-I-94 and genipin were identified as potential JMJD8-targeted drugs and may be effective for temozolomide- and nitrosourea-resistant GBMs as alternative therapies.