The N6-Methyladenosine-Modified Pseudogene HSPA7 Correlates With the Tumor Microenvironment and Predicts the Response to Immune Checkpoint Therapy in Glioblastoma

Human GBM and normal brain tissues for m⁶A-seq were obtained from patients admitted to Qilu Hospital. All participants provided written informed consent, and the research was approved by the Ethics Committee on Scientific Research of Shandong University Qilu Hospital (approval number: KYLL-2018-324).

The RNA sequencing (RNA-seq) transcriptome, somatic mutation data, and corresponding clinicopathological parameters of the TCGA GBM cohort were obtained from the TCGA database (http://cancergenome.nih.gov/↗). Two Chinese Glioma Genome Atlas (CGGA) GBM RNA-seq datasets and the corresponding clinicopathological parameters were obtained from the CGGA database (http://www.cgga.org.cn/↗). For ICB data, genomic and clinical information from the IMvigor210 cohort, complete expression data, and detailed clinical annotations were obtained from http://research-pub.Gene.com/imvigor210corebiologies↗ based on the Creative Commons Attribution 3.0 license. For data from GBM patients treated with PD-1 inhibitors (pembrolizumab or nivolumab), clinical information was obtained from Supplementary paper data, and the sequencing data were obtained from SRA PRJNA482620 (2). For patients with melanoma treated with anti-CTLA4 therapy, the expression data were downloaded from the cBioPortal database (http://www.cbioportal.org/↗), and the detailed clinical characteristics of individual patients were obtained from the supplementary data of a previous paper (23). In addition, the somatic mutation data and information in Figure S11A were obtained from the cBioPortal database. The m⁶A-seq sequencing data have been deposited in SRA PRJNA661159 (the data are being processed, submission ID: SUB8069560, released when the paper is published). The processed data are available from the corresponding author upon reasonable request.

We used the single-sample gene set enrichment analysis (ssGSEA) algorithm to calculate the enrichment score of immune cell infiltration into the GBM TME for each sample. Immune cell-related genes were obtained from Bindea et al. (24) and Robert L (25). (Table S1) and included genes related to immune cell types, immune-related pathways and functions. Based on the ssGSEA results, samples from the TCGA GBM cohort were classified into the high immune cell infiltration (immune-H) group or low immune cell infiltration (immune-L) group by using the “hclust” R package.

To explore the differences in biological behavior among the samples with distinct HSPA7 expression levels, we used selected HALLMARK (26) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (27) from gene sets from the Molecular Signatures Database (MSigDB) and other commonly used gene signatures (28–30) to estimate pathway enrichment scores for each sample (Table S2) by gene set variation analysis (GSVA) using the “GSVA” R package.

The terms enriched with genes positively correlated with the expression of HSPA7 and genes interacting with HSPA7 detected in NCBI (https://www.ncbi.nlm.nih.gov/gene/3311↗) and the starBase database (http://starbase.sysu.edu.cn/index.php↗) were analyzed via the Metascape resource (http://metascape.org/gp/index.html#/↗).

All patient-derived GSC cell lines, including mesenchymal (MES) subtype GSC cell lines (GSC 20 and GSC 267), proneural (PN) subtype GSC cell lines (GSC 8–11), and neural progenitor cells (NPCs) were kindly donated by Dr. Frederick F. Lang and Dr. Krishna P.L. Bhat (The University of Texas, M.D. Anderson Cancer Center, Houston, TX, USA). The cells were cultured in DMEM/F12 supplemented with B27 (Invitrogen,California, USA), 20 ng/ml EGF (R&D Systems, USA), and 20 ng/ml bFGF (R&D Systems, California, USA). The human glioma cell lines U87MG, U251MG, A172, and LN229 and the human monocyte cell line THP-1 were obtained from the Chinese Academy of Sciences Cell Bank. U87MG, U251, A172, and LN229 cells were cultured in DMEM (Thermo Fisher Scientific, USA) supplemented with 10% FBS (Thermo Fisher Scientific; Waltham, MA, USA). THP-1 cells were cultured in RPMI-1640 (Thermo Fisher Scientific) supplemented with 10% FBS (Thermo Fisher Scientific). THP-1 cells were incubated with 100 ng/ml PMA (Sigma-Aldrich; St. Louis, MO, USA) for 24 h in vitro to induce their differentiation into macrophages. Cells were cultured in a standard humidified atmosphere of 5% CO₂ at 37°C.

Protein was extracted from GSC cells or glioma cells. The following primary antibodies were used: GAPDH (Cell Signaling Technology, Boston, USA, 5174), YAP1 (Cell Signaling Technology, 14074), CD44 (Cell Signaling Technology, 3570), LOX (Abcam, Cambridge, UK, ab174316), and YKL40 (Cell Signaling Technology, 47066).

ASOs were synthesized by RiboBio (Guangzhou, China). HSPA7/FTO overexpression, METTL3 shRNA and corresponding control lentiviruses were synthesized by GeneChem (Shanghai, China). Target sequences for HSPA7 were as follows: ASO#1: 5′-GGAAGCGGAGCTGAGCAGAT-3′; ASO#2: 5′-CTAACAAGATCACCAATGAC-3′. Target sequences for METTL3 were as follows: 5′-GCCAAGGAACAATCCATTGTT-3′.

The slides were washed with PBS for 15 min and blocked with 10% goat serum in PBS. The slides were incubated overnight in a humidified chamber at 4°C with the following primary antibodies: LOX (Abcam, ab174316); YAP1 (Cell Signaling Technology, 14074), CD68 (Abcam, ab213363), CD44 (Cell Signaling Technology, 3570), Ki67 (Cell Signaling Technology, 9449), PD-L1 (Abcam, ab213524), and SPP1 (Abcam, ab8448). After primary antibody incubation, the samples were washed with PBS and incubated with the matching fluorescent-conjugated secondary antibody (1:500 dilution, Thermo Fisher) at room temperature for 1 h. Images were captured using a LeicaSP8 confocal microscope (Leica Microsystems, Wetzlar, Germany)

RNA-FISH was performed according to the instructions of the manufacturer (GenePharma, Shanghai, China). Incubated with a cy3-labeled HSPA7 probe, nuclei were counterstained with DAPI. Images were captured using a LeicaSP8 confocal microscope. Probe sequence: ATCCTTTTGCACCTCCCCGACCC; AACCTTCCCGCACCTTCCCGCCCAGTC.

Glioblastoma organoid models were generated as previously described (31). GBO medium containing 50% DMEM:F12 (Thermo Fisher Scientific), 50% Neurobasal (Thermo Fisher Scientific), 1× GlutaMax (Thermo Fisher Scientific), 1× NEAAs (Thermo Fisher Scientific), 1× PenStrep (Thermo Fisher Scientific), 1× N2 supplement (Thermo Fisher Scientific), 1× B27 w/o vitamin A supplement (Thermo Fisher Scientific), 1× 2-mercaptoethanol (Thermo Fisher Scientific), and 2.5 μg/ml human insulin (Sigma) per well was placed on an orbital shaker rotating at 120 rpm within a 37°C, 5% CO₂, and 90% humidity sterile incubator.

Transwell assays were performed in 24-well multiwell insert systems according to the protocol of the manufacturer. THP-1 cells were incubated with 100 ng/ml PMA (Sigma-Aldrich) for 24 h in vitro to induce their differentiation into macrophages and then added to the top chamber in serum-free media. The bottom chamber was filled with 10% FBS 1640 and GSC growth media. After 24 h of incubation, the top chamber cells were removed using a cotton swab, and the membrane was fixed in 4% paraformaldehyde for 15 min and stained with crystal violet for 15 min. Five fields of adherent cells in each well were photographed randomly.

Kaplan–Meier survival analysis was performed using GraphPad Prism 7.04, and significant differences between two groups were compared by the log-rank (Mantel–Cox) test. The waterfall function in the “maftools” package was used to visualize the mutational landscape in patients in the high- and low-HSPA7 groups or the high- and low-tumor mutation burden (TMB) groups. Student’s t-test was used for two-group comparisons. For comparisons among more than two groups, the Wilcoxon test and one-way ANOVA were used for non-parametric and parametric data (32). P > 0.05 was considered nonsignificance (ns), P ≤ 0.05 was considered statisticallysignificant (*P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001). All data processing with R packages was performed using R Studio (version 3.6.3)

To understand the pattern of the lncRNA m⁶A modification profiles in GBM, three human GBM tumor tissues and three normal brain tissues were subjected to transcriptome m⁶A-seq, which revealed that a considerable proportion of lncRNAs were extensively m⁶A-modified. In addition, 2,113 m⁶A peaks were identified in the normal group, corresponding to the transcripts of 1,538 genes, including 544 long intergenic ncRNAs (lincRNAs), 720 antisense lncRNAs, 262 pseudogenes, and 12 others. In the GBM group, 2,412 m⁶A peaks were identified, corresponding to the transcripts of 1,508 genes, namely, 445 lincRNAs, 771 antisense lncRNAs, 283 pseudogenes, and 9 other genes (Figure 1A). Further analysis showed that most of the lncRNAs (74.5% of the m⁶A-methylated genes in the normal group but 65.3% of the methylated genes in the GBM group) in both groups contained only one peak, while a relatively small number of lncRNAs contained two peaks, and very few lncRNAs contained three or more peaks (Figure 1B). Upon further analysis of the distribution profiles of m⁶A peaks within lncRNAs, we found that m⁶A sites were distributed in almost the same proportion in intronic and exonic regions in the normal group but exhibited a slightly increased tendency to be distributed in exonic regions in GBM tissues (Figure 1C).

To reveal the significance of m⁶A-methylated lncRNAs in GBM, the differences and overlaps in the m⁶A-modified lncRNAs between the GBM and normal brain groups were analyzed by constructing a Venn diagram. As shown in Figure 1D, 781 m⁶A-modified lncRNAs were common to both groups. A total of 727 new genes were expressed, and the expression of 757 genes was lost in the GBM group compared with the normal group, indicating a significant difference in global lncRNA m⁶A modification patterns between the GBM and normal groups. To explore the effect of m⁶A on lncRNA expression, the differentially expressed lncRNAs between 153 primary GBM and five normal brain tissues in the TCGA GBM cohort were compared. Compared with normal samples, GBM tissues exhibited 4,375 differentially expressed lncRNAs (LogFC ≥ 1 and padj ≤ 0.05), with 2,614 upregulated and 1,761 downregulated (Table S3). In addition, the global abundance of m⁶A peaks between GBM and normal brain tissues was also compared. Furthermore, integrated analysis of these differentially m⁶A-modified lncRNAs and differentially regulated lncRNAs from the TCGA dataset was conducted. The distribution of genes with a significant change in both the m⁶A level (|FC| ≥ 1.2, P ≤ 0.05) and the overall transcript expression level (|FC| ≥ 2, padj ≤ 0.05) is shown in Figure 1E. These genes were divided into four main groups: 208 were hypermethylated and upregulated (“hyperup”), 394 were hypomethylated and downregulated (“hypodown”), 56 were hypermethylated but downregulated (“hyperdown”) and 64 were hypomethylated but upregulated (“hypoup”) in GBM tissues relative to normal brain tissues. We also discovered a positive correlation between differentially methylated m⁶A peaks and the expression levels of their corresponding genes [Figure 1E, Pearson correlation coefficient (r) = 0.495, P<0.0001], and the m⁶A modification sites in all of the above four groups of genes were distributed in exonic regions, reconfirming that m⁶A modification can regulate the expression of mature lncRNAs. These results revealed obviously distinct m⁶A modification patterns between GBM and normal brain tissues. Moreover, m⁶A modification could regulate numerous lncRNAs, possibly by regulating their stability, degradation or other functions, as reported previously.

To investigate the effects of these m⁶A-regulated lncRNAs on cell infiltration into the TME, we assessed the tumor purity, stromal score, and immune score of 153 primary GBM cases in the TCGA GBM cohort using the ESTIMATE algorithm (see Materials and Methods, Table S4). Then, we performed differential analysis of all RNA-seq data from these 153 GBM samples in the TCGA database based on the median cutoff immune/stromal scores. The volcano plot of the high/low stromal/immune scores revealed differential gene expression profiles between the samples. A total of 791 upregulated lncRNAs and 459 downregulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were identified based on the difference in immune scores (Supplementary Figure S1A and Table S5). Simultaneously, 734 upregulated lncRNAs and 269 downregulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were identified based on the differential analysis of stromal scores (Figure S1B and Table S6). As the Venn diagram (Supplementary Figures S1C, D) indicates, four identical upregulated genes and 11 identical downregulated genes were related to immune activation, stromal activation, and m⁶A modification. Then, we performed Kaplan–Meier analysis on patients stratified by the expression levels of these 15 differentially expressed genes, and only HSPA7 (Figure 2B) and AC011899.9 (Supplementary Figure S2B) were found to have prognostic significance (P ≤ 0.05) in the 153 TCGA GBM samples based on the median expression levels. The Gene Expression Profiling Interactive Analysis (GEPIA) database was used to show that HSPA7 was significantly overexpressed in GBM tissues compared with normal brain tissues in the Genotype-Tissue Expression (GTEx) database (Figure 2A), while AC011899.9 was not (Supplementary Figure S2A). Thus, we focused only on HSPA7 in the following work.

We then tested the RNA expression level of HSPA7 in GBM cell lines (U87MG, U251MG, A172, and LN229 cells), and discovered that its expression was significantly higher in GBM cell lines than in normal human astrocyte (NHA) cells (Figure 2C). The m⁶A methylation peak distribution and abundance in HSPA7 transcripts from GBM and normal brain tissues, as detected by m⁶A-seq, were visualized using IGV software (Figure 2D). We found that HSPA7 was highly enriched in the m⁶A-precipitated fraction, and the m⁶A modification enrichment level could be regulated by methyltransferase-like 3 (METTL3), which contains a catalytic activity domain to catalyze m⁶A formation (Figure 2E). Additionally, the expression of HSPA7 was significantly inhibited by knocking down METTL3 (Figure 2F) and overexpressing FTO, an m⁶A demethylase (Figure 2G). These results provided evidence that HSPA7 harbored high m⁶A modification levels and thus could be regulated in an m⁶A-dependent manner.

To explore the association between HSPA7 expression and clinical characteristics, we first compared HSPA7 expression levels in patients in the TCGA GBM cohort stratified separately by molecular subtype, IDH1 status, CpG island methylator phenotype (G-CIMP) status, MGMT promoter status, age, and sex. As shown in Figure 2H, HSPA7 expression in the PN subtype was significantly lower than the HSPA7 expression in the classical (CL) and MES subtypes and was highest in mesenchymal samples. HSPA7 expression was lower in samples with IDH mutations than in samples with wild-type IDH1. Regarding the G-CIMP status, HSPA7 expression was lower in patients with G-CIMP tumors than in those without G-CIMP tumors. No obvious correlations between HSPA7 expression and MGMT promoter methylation, age, or sex were observed. Furthermore, univariate Cox regression analysis of the overall survival of GBM patients in the TCGA cohort showed that high HSPA7 expression (HR: 1.511, P = 0.056) was an independent risk factor associated with the prognosis of GBM. Moreover, high HSPA7 expression (HR: 1.418, P = 0.034) remained a statistically significant factor in GBM patients after adjustment for age, sex, IDH status, MGMT promoter methylation status, and G-CIMP status in subsequent multivariate Cox regression analysis (Figure 2I). These results indicated that the pseudogene HSPA7 is a novel risk prognostic biomarker and indicates therapeutic outcomes of GBM patients.

To gain further insight into the exact role of HSPA7 in immunophenotype determination, we analyzed 31 immune-associated gene sets representing diverse immune cell types, functions, and pathways (see Materials and Methods). As shown in Figures 3A, B, the HSPA7-high expression group had significantly greater cell infiltration into the TME, higher immune and stromal scores, and lower tumor purity than the HSPA7-low expression group, confirming that HSPA7 could indeed regulate immune cell infiltration and immune-related gene expression. We then explored the specific differences in 31 immune cell phenotypes with high and low HSPA7 expression. Compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased infiltration of immunosuppressive cell populations, such as macrophages, neutrophils, and Tregs. However, some immune-activating cells, including activated dendritic cells (aDCs), immature DCs (iDCs), plasmacytoid DCs (pDCs), and tumor-infiltrating lymphocytes (TILs), were also enriched (Supplementary Figure S3A), indicating the complexity of the TME, in which GBM cells elicit multiple biological behavioral changes through direct or indirect interactions with other TME components. Two recent publications (7, 8) on single-cell mapping of human brain cancer reported that tumor-associated macrophages (TAMs) are the most abundant cellular components in brain TMEs, which can be subdivided ontogenetically into tissue-resident microglia (MGs) and macrophages of embryonic origin, and bone marrow-derived macrophages (BMDMs), showing tumor-promoting and immunosuppressive functions. To identify specific macrophage cell populations linked to HSPA7 expression in GBM, we examined the TCGA GBM dataset for TAM MG and TAM BMDM using validated gene set signatures (25), which demonstrated that tumors with high HSPA7 expression exhibited significantly increased infiltration of TAM BMDMs, while TAM MG was significantly decreased (Figure 3A and Supplementary Figure S3A), suggesting that HSPA7 enhanced recruitment of tumor-promoting BMDMs into the GBM TME rather than MG. These studies also showed that both activation and exhaustion of lymphocytes were prevalent, with increased relative frequencies of Tregs in the brain tumor TME. Moreover, we found elevated neutrophil infiltration in tumor tissues, revealing the complexity and multifaceted functions of the brain TME. We found that HSPA7 had a significant positive correlation with the enrichment scores of the three immunosuppressive immune cells and immune checkpoints (Figure 3C). We next investigated chemokines and immune modulators associated with immune suppression states. We then explored the specific differences in markers of myeloid lineages with suppressive functions (CD33, NOS2, CD163, CD68, SPP1, CD14, CD206), immune inhibitory checkpoints (LAG3, CTLA4, HAVCR2, PDCD1, PDCD1LG2, CD27, TIGIT, CD274, and TNFRSF9), and major neutrophil-recruiting chemokines and their receptors (ITGA3, CD177, MET, and CXCL8) between patients with high and low HSPA7 expression. Tumors with high HSPA7 expression exhibited significantly increased expression levels of these markers compared to tumors with low HSPA7 expression (Figures 3D–F), and positive correlations were found between HSPA7 expression and these molecules (Supplementary Figures S4A–C). Furthermore, to explore the direct involvement of the HSPA7 in the biological pathway causing immune suppression, we then analyzed myeloid cell-derived macrophage-restricted chemokines, which were the main factors that cause immunosuppression in GBM. We found that compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased myeloid cell-derived macrophage-restricted chemokines (8, 25) (Supplementary Figures S4D, E), including CCL17, CXCL2, CXCL3, and CXCL16 (involved in wound healing); immunosuppressive cytokines IL-10, CTSB, and CTSW (participating in multiple tumor-promoting processes, including invasion and metastasis); CCL2 and CCR2 (involved in macrophage chemotaxis); and other chemokine contributions of myeloid cell populations to the inflammatory TME milieu, indicating that HSPA7 could promote immunosuppressive phenotypes and suppress intratumoral antitumor immune responses.

To further understand the exact role of HSPA7 in determining the TME profile and immunophenotype, we performed unsupervised consensus clustering based on the TME cell populations and immune-related functional gene sets identified by Bindea et al. (24). This analysis identified a TME pattern with two clusters corresponding to a TME immune-low (immune-L) and a TME immune-high (immune-H) phenotype (Supplementary Figure S3B). Patients in the immune-H group (n = 51) exhibited poorer prognosis (Figure 3G) and higher HSPA7 expression levels (Figure 3H) than patients in the immune-L group (n = 102). The PN, CL, and MES subtypes of GBM have been most consistently described in the literature; the PN subtype is related to a more favorable outcome, and the MES subtype is related to poorer survival (23). The association between the MES gene expression signature, characterized by NF1 mutation, with reduced tumor purity, elevated invasion, enhanced migration capacity, and infiltration of immunosuppressive cells (macrophages, microglia, mesenchymal stem cells, or other cells) has been identified as a common theme across cancers (33, 34). Our results showed that most MES GBM samples had high expression of HSPA7 as well as the immune-H phenotype (Figures 3I, J). Moreover, we observed the same results as all GBM samples by analyzing MES glioblastoma alone, in which tumors with high HSPA7 expression exhibited significantly increased infiltration of immunosuppressive cell populations and regulators, such as macrophages, TAM BMDMs, neutrophils, Tregs, and immune checkpoints (Supplementary Figures S5A, B). As MES glioblastomas are strongly associated with higher immunosuppressive cell infiltration, we hypothesized that HSPA7 activated the immune microenvironment by promoting the phenotypic transformation of the GBM MES subtype. As shown in Supplementary Figure S5C, compared to tumors with low HSPA7 expression, tumors with high HSPA7 expression exhibited significantly increased expression levels of genes in the MES phenotype signature genes (33). Furthermore, HSPA7 positively correlated with CD44, a marker of the MES subtype (Figure 3K). GSEA also showed that GBM samples with high HSPA7 expression were enriched in the MES subtype compared to GBM samples with low HSPA7 expression (Figure 3L). In addition, we found that HSPA7 expression was higher in MES GSCs than in PN GSCs (Figure 3M). Moreover, Western blot assays revealed that knockdown of HSPA7 reduced the expression of CD44 (Figures 3N, O). Overall, these results indicated that in GBM, HSPA7 may be a robust indicator of the immunophenotype and be significantly correlated with a poorer immune response.

To explore the differences in biological behaviors among these samples with distinct HSPA7 expression levels, we used the GSVA algorithm to estimate pathway enrichment scores for each sample (see Materials and Methods). Compared to the low HSPA7 expression group, the high HSPA7 expression group exhibited marked enrichment of stromal activation pathways [angiogenesis, epithelial–mesenchymal transition (EMT), VEGF, and TGF-β signaling pathways], oncogenic signaling pathways (hypoxia, apoptosis, PI3K Akt MTOR signaling, glycolysis, KRAS signaling, and other pathways), and immune responses (IL6 Jak stat3 signaling, which exerts immunosuppressive effects on T cell function and mediates ICB resistance in cancers (35), inflammatory response; interferon response; complement; and other pathways). However, the high HSPA7 expression group exhibited lower enrichment in pathways related to DNA replication- and DNA damage response-related functions (Figure 4A). Previous studies demonstrated that activation of stromal cells in the TME can induce T cell suppression and oncogenic signaling pathway activation via complex cellular and biological reconfiguration mechanisms (36, 37), attenuating the tumor response to PD-L1 blockade (28). As shown in Figure 4B, stromal activation pathways were positively correlated with carcinogenic signaling pathways and immune suppression pathways but negatively correlated with DNA damage response and repair pathways in TCGA GBM samples. Further analyses of the activity of stroma-related pathways indicated that high HSPA7 expression was significantly associated with higher stromal activation signatures (Figure 4C), as constructed by Mariathasan et al. (28), positively correlated with fibroblast enrichment (Figure 4D), as calculated by the MCP-counter method (38), and positively correlated with typical stromal cell activation-related pathways (Figure 4E). Furthermore, TGF-β is a pleiotropic cytokine associated with poor prognosis in GBM, playing a protumorigenic role by promoting immunosuppression, angiogenesis, metastasis, mesenchymal transition, and fibroblast activation (28, 39–41), and has been proven to promote the progression of GBM via an autocrine signaling loop (42). In our analysis, a TGF-β ligand (TGFB1) and a TGF-β receptor (TGFBR2), two key regulators in stromal activation and EMT pathways (28), and other genes encoding ECM and matricellular proteins (COLA4A1, COLA4A2, ECM1, and FN1) (8), which form a barrier to lymphocyte infiltration, showed increased expression in the high-HSPA7 group compared to the low-HSPA7 group (Figure 4F). These results suggested that immune cells and stromal cells in the TME can cooperate to synergistically regulate the immunosuppressive microenvironment, thereby promoting tumor immune escape.

Then, we explored the pathways enriched with genes positively correlated with HSPA7 (Pearson’s r ≥ 0.3, P ≤ 0.05, Table S7) via the Metascape database (see Materials and Methods). These genes were significantly enriched in pathways related to the immune response, including cytokine-mediated signaling pathways, leukocyte migration, differentiation, and other immune-related pathways. Additionally, they were enriched in pathways involving stromal activation such as extracellular structure organization and response to wounding. These enriched pathways interacted with each other to form a protein–protein interaction network (Figure 5A). Furthermore, enrichment analysis in the PaGenBase (43) showed that these genes were almost completely specifically expressed in spleen, blood, bone marrow, and some other tissues with peripheral immune cell aggregation (Figure 5B), suggesting that HSPA7 can indeed promote the infiltration of immune cells into tumor tissues. In addition, enrichment analysis in the TF-target interaction database Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) (44) showed that the genes were substantially transcriptionally regulated by NF-κB1 and RELA, two important transcription factors involved in NF-κB signaling pathways (Figure 5C). Further GSEA using the Gene Ontology (GO) and KEGG databases also showed that HSPA7 expression was positively correlated with the immune response and extracellular structure organization (Figure 5D).

To further validate the function of HSPA7 in GBM, we explored its expression pattern in two CGGA RNA-seq cohorts. The expression of HSPA7 was highest in GBM (WHO IV) among the glioma specimens with three different WHO grades (,). Kaplan–Meier survival analysis showed that GBM patients in both cohorts with high HSPA7 expression had poor survival outcomes (,). Further GSVA enrichment analysis also showed that HSPA7 expression was correlated with the immunophenotype, stromal activation, and oncogenic pathway activation in GBM samples (,). These results confirmed that HSPA7 can be used as an independent prognostic risk biomarker in GBM patients. Moreover, HSPA7 can regulate the TME immune response and stromal activation, promoting malignant progression of GBM tumors. 1 1 1 1 1 1

To explore the potential molecular mechanism regulating HSPA7 RNA metabolism, we first screened HSPA7-interacting proteins via the NCBI database (Table S8). We verified that WTAP (45), a regulatory subunit of the m⁶A methyltransferase (46, 47), is required for the localization of the m6A methyltransferase into nuclear speckles as a putative HSPA7-binding regulator. To better map the association between WTAP and HSPA7, we explored the expression pattern in different grades of glioma in the TCGA, CGGA, Gravendeel and Rembrandt datasets. As shown in Supplementary Figure S8A, the expression of WTAP was highest in the WHO IV (GBM) group among glioma samples of three WHO grades and normal brain samples. Furthermore, overexpression of WTAP correlated with poor overall survival of GBM patients in all four datasets (Supplementary Figure S8B). Moreover, the expression of WTAP was positively correlated with the expression of HSPA7 in the TCGA GBM dataset (as calculated via the GEPIA database) and two CGGA GBM datasets (Supplementary Figure S8C). Further GSVA enrichment analysis showed that compared to the WTAP-low group, the WTAP-high group showed increased immune cell infiltration into the TME, stromal activation and oncogenic pathway activation, accompanied by higher immune and stromal scores and lower tumor purity, in the TCGA GBM cohort (Supplementary Figure S8D), CGGA GBM cohort 1 (Supplementary Figure S9A) and CGGA GBM cohort 2 (Supplementary Figure S9B). These patterns were the same as the patterns observed for HSPA7 expression. The above findings indicate that m⁶A methylation of HSPA7 is regulated by WTAP, although this conclusion and the functional mechanisms need further demonstration. Subsequent pathway and process enrichment analyses of HSPA7-interacting proteins showed that they can mediate many biological behaviors, such as mitotic cell cycle phase transition, neddylation, proteasome-mediated ubiquitin-dependent protein catabolism, HIF1α pathway activity, and other functions (Supplementary Figure S10A, Table S8). Then, we further explored the enriched terms across the RNA-binding proteins (RBPs) detected by crosslinking immunoprecipitation coupled with sequencing (CLIP-seq) from data deposited in the starBase database. These RBPs were enriched mainly in terms such as mRNA processing, nucleic acid transport, and RNA catabolic process (Supplementary Figure S10B and Table S9). However, the exact mechanism by which HSPA7 regulates the GBM TME requires further clarification.

We desired to further investigate the capacity of HSPA7 to predict the response to immune checkpoint therapy in GBM. However, few published immunotherapy datasets include GBM patients who received ICB therapy. Thus, we used a cohort of melanoma patients who received anti-CTLA4 therapy (23) and urothelial cancer cohorts of patients who received anti-PD-1 therapy (IMvigor210) to perform a complementary evaluation of the ability of HSPA7 to predict the immunotherapy response. Patients in the anti-CTLA4 cohort with low HSPA7 expression exhibited significant clinical benefits [Figure 6A, anti-CTLA cohort, HR 1.927 (0.8997–4.129)]. Furthermore, compared to the HSPA7-high group, the HSPA7-low group exhibited a significant therapeutic benefit and clinical response to antiCTLA-4 immunotherapy (Figures 6B, C). Similar results were found in the anti-PD-L1 cohort. Kaplan–Meier analysis showed that the patients with the 50 lowest HSPA7 expression levels exhibited markedly prolonged survival compared to the patients with the 50 highest HSPA7 expression levels [Figure 6D, HR 1.670 (1.022–2.792)]. Additionally, the HSPA7-low group exhibited significantly better therapeutic and clinical responses to anti-PD-L1 immunotherapy than the HSPA7-high group (Figures 6E, F).

Further GSVA enrichment analysis showed that compared to the HSPA7-low group, the HSPA7-high group exhibited higher TME immune cell infiltration with higher immune and stromal scores, lower tumor purity (Supplementary Figure S11A), higher stromal activation, higher oncogenic pathway activation, and lower MMR pathway activation (Supplementary Figure S11B), completely consistent with the functional enrichment patterns in GBM (Figures 3A, 4A and Supplementary Figures S5C, S6C). PD-L1 expression in immune cells (ICs) and tumor cells (TCs) was also assessed in the IMvigor210 cohort, and we examined the difference in HSPA7 expression among groups with different PD-L1 expression levels. As indicated in Figures 6G, H, patients with higher PD-L1 expression levels in either immune cells or tumor cells exhibited higher HSPA7 expression levels (ANOVA summary IC: P = 0.0001, TC: P = 0.0009), indicating that HSPA7 can upregulate PD-L1 expression, suppressing immune activation. Moreover, previous studies indicated that F-box and WD repeat domain containing 7 (FBXW7) is a vital tumor suppressor in various cancers, controlling proteasome-mediated degradation of oncoproteins such as cyclin E, c-Myc, Mcl-1, mTOR, Jun, and Notch (48) and that its loss-of-function mutation promotes resistance to anti-PD-1 therapy through downregulation of vital sensing pathways (49). Compared to the FBXW7 wild-type group, HSPA7 expression was obviously decreased in the FBXW7 mutant group (Figure 6I). Fibroblast growth factor receptor 1 (FGFR1) is frequently mutated in various tumors, and inhibitors of FGFR1 have shown promising therapeutic value in several preclinical models (50). Palakurthi S. et al. (51) found that the combination of FGFR inhibition and PD-1 blockade can promote tumor-intrinsic induction of antitumor immunity. We also found that HSPA7 expression was significantly higher in the FGFR mutant group than in the FGFR wild-type group (Figure 6J). Collectively, our work strongly indicates that HSPA7 expression contributes to predicting the response to immune checkpoint therapy.

We then analyzed the differences in the distribution of somatic mutations between the low and high HSPA7 expression groups in the TCGA GBM cohort using the “maftools” package and found that the PTEN mutation rate was significantly increased in the HSPA7-high group compared to the HSPA7-low group (Figure 6K; low: 20%, high: 36%). Additionally, the NF1 mutation rate, which often occurs in the MES subtype and drives recruitment and activation of TAMs (33), was also obviously elevated in the HSPA7-high group compared to the HSPA7-low group (Figure 6K; low: 8%, high: 11%). Chen et al. (52) found that PTEN deficiency in GBM increases macrophage infiltration by activating YAP1 signaling, which directly upregulates lysyl oxidase (LOX) expression, an MES subtype marker and a macrophage chemoattractant (33), in turn inducing SPP1 secretion to support GBM survival. We found that YAP1 signaling was significantly upregulated in the HSPA7-high group compared to the HSPA7-low group (Figure 7A). We also confirmed that LOX and YAP1 were higher in MES GSCs than in PN GSCs (Figure 7B). Moreover, LOX and YAP1 expression levels were enhanced upon overexpression of HSPA7 in PN subtype GSCs 8–11 (Figures 7C, D) but decreased significantly upon knockdown of HSPA7 in MES subtype GSCs 267 (Figure 7D), suggesting that HSPA7 may facilitate tumor promoting macrophage infiltration by enhancing LOX expression, which could promote macrophage migration into the GBM TME and enhance angiogenesis by upregulating TAM-derived SPP1 (52). Finally, to further confirm the role of HSPA7 in facilitating macrophage migration by enhancing the YAP1-LOX axis, using the transwell assay, we found that conditioned medium from HSPA7 knockdown GSCs 267 inhibited THP1-differentiated macrophage migration significantly compared to the NC group (Figure 7E). Then, THP1-differentiated macrophages were cultured in conditioned medium collected from GSCs 267 transfected with HSPA7 ASO or NC. As shown in Figure 7F, conditioned medium from HSPA7-knockdown GSCs inhibited SPP1 expression in macrophages compared to the NC group. The immunofluorescent staining of our local clinical tumor tissue section analyses also showed that the expression of YAP1, LOX, CD68 (human macrophage marker), SPP1, and PD-L1 was enhanced in HSPA7-high GBM tissues compared to HSPA7-low tissues (Figure 7G).

Moreover, Robert M. Samstein et al. (53) found that TMB predicted survival after immunotherapy across multiple cancer types (including GBM). Thus, we downloaded the mutation data of 82 GBM samples and then analyzed the differences in the distribution of somatic mutations between the samples with the 30 lowest TMBs and the 30 highest TMBs. The PTEN and NF1 mutation rates were significantly increased in low-TMB samples compared to the high-TMB samples, and the same results were found for the HSPA7 expression level (Supplementary Figure S12A). We further explored the association between the HSPA7 expression level and TMB. As shown in Supplementary Figure S12B, HSPA7 expression was statistically negatively correlated with TMB. In addition, Touat M. et al. (54) recently found that mismatch repair (MMR)-deficient gliomas were characterized by a lack of prominent T cell infiltration, extensive intratumoral heterogeneity, poor patient survival, and a low rate of response to PD-1 blockade therapy. Our data also showed that the high HSPA7 expression group showed reduced activation of MMR-associated pathways (Figure 4A), which was negatively correlated with immunosuppression, stromal activation, and oncogene pathway activation (Figure 4B). Other studies have also shown the same results (28, 55). Furthermore, HSPA7 expression was negatively correlated with the expression of MMR genes such as MLH1, MSH2, MSH6, and PMS2 (Supplementary Figure S12C). Previous studies indicated that GBM patients exhibit resistance to immune checkpoint inhibitors (ICIs) due to the low mutation rate, the PTEN-deficient immunosuppressive microenvironment, infiltration by myeloid-derived suppressor cells (including TAMs, prominent players in brain cancers), and activation of tumor stromal cells (56, 57), indicating that knockdown of HSPA7 may enhance the efficacy of ICIs such as PD1 inhibitors in GBM. To illustrate this function, we generated patient-derived glioblastoma organoids (GBOs), which maintain cell-type heterogeneity (including macrophages, T cells, and vascular cells) and molecular signatures of their respective parental tumors, according to experimental protocols already reported by Jacob F. et al. (31). Immunofluorescence assays of our GBO sections also confirmed the presence of macrophages and vascular cells, as detected by CD68 and CD31 markers, respectively (Supplementary Figure S13A), supporting the use of three-dimensional GBO models for holistic study of the TME of GBM and evaluating the efficacy of PD1 inhibitors. We next applied our GBO model to test treatment responses in vitro. To mimic the postsurgical standard of care treatment, we subjected GBO samples from patient whose GBM tissues expressed high CD68 and SPP1 (Supplementary Figure S13B) to a single exposure conditioned medium from HSPA7-knockdown GSCs or with PD1 inhibitor (5 µM) treatment for 5 days. The therapeutic response was evaluated by quantifying the percentage of cells expressing Ki67 (proliferation index) and CD44 (invasion index). As shown in Figure 7H, knockdown of HSPA7 significantly enhanced the efficacy of the PD1 inhibitor in GBM. In summary, HSPA7 facilitated tumor promoting macrophage migration into the GBM TME by activating the YAP1–LOX axis, and our results indicated that HSPA7 might be a potential immunotherapy target for GBM patients.

We investigated whether HSPA7 expression is correlated with prognosis in a pancancer patient cohort. The expression of HSPA7 in tumor tissues and GTEx normal brain tissues was determined in the GEPIA database. As shown in Supplementary Figure S14, the expression of HSPA7 was significantly lower in tumor tissues in the ACC, COAD, DLBC, LAML, LUAD, LUSC, READ, THCA, and THYM cohorts but significantly higher in tumor tissues in the GBM, KIRC, KIRP, and PAAD cohorts than in the corresponding normal tissues. Genetic alterations were also observed in many other tumors, although these alterations were not significant. In addition, Cox regression analysis of survival rates and Kaplan–Meier analysis were conducted using the SangerBox tool. The relationships between HSPA7 expression levels and the prognoses of different tumors are shown in Supplementary Figures S15-S17. Statistically, HSPA7 is a risk factor in most cancers (such as ACC, LDBC, and CESC). However, we found that it is a relatively favorable factor in several cancer types such as SKCM, KIRP, and SARC. We next evaluated the immune score and stromal score across 33 cancers in the TCGA database. HSPA7 expression was significantly positively correlated with the immune score (Supplementary Figure S18) and stromal score (Supplementary Figure S19) in all 33 cancer types. Moreover, HSPA7 was significantly positively correlated with immune checkpoint expression (Supplementary Figure S20A) and the immune response (Supplementary Figure S20B) in most cancer types such as GBM, LGG, OV, and LUAD. However, no relationship was found in other cancers such as SKCM, UVM, TGCT, and PAAD, indicating that HSPA7 may act via a mechanism other than TME regulation in these cancers. Thus, the molecular mechanisms of this gene need further study in specific tumors. Collectively, these data emphasize that HSPA7 may be a key factor in facilitating the acquisition of various immunophenotypes in various cancers and may be considered in the development of more effective immunotherapies for GBM and other immunotherapy-resistant tumors. Because different immune cells may play different roles in different tumors, the functions of HSPA7 in specific tumors must be explored.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/↗, SRA PRJNA661159.