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

Background: Glioblastoma (GBM), one of the most aggressive tumors of the brain, has no effective or sufficient therapies. Identifying robust biomarkers for the response to immune checkpoint blockade (ICB) therapy, a promising treatment option for GBM patients, is urgently needed.
Methods: We comprehensively evaluated lncRNA m6A modification patterns in m6A-sequencing (m6A-seq) data for GBM tissues and systematically investigated the immune and stromal regulators of these m6A-regulated lncRNAs. We used the single-sample gene-set enrichment analysis (ssGSEA) algorithm to investigate the difference in enriched tumor microenvironment (TME) infiltrating cells and the functional annotation of HSPA7 in individual GBM samples. Further, we validated that HSPA7 promoted the recruitment of macrophages into GBM TME in vitro, as well as in our GBM tissue section. We also explored its impact on the efficacy of ICB therapy using the patient-derived glioblastoma organoid (GBO) model.
Results: Here, we depicted the first transcriptome-wide m6A methylation profile of lncRNAs in GBM, revealing highly distinct lncRNA m6A modification patterns compared to those in normal brain tissues. We identified the m6A-modified pseudogene HSPA7 as a novel prognostic risk factor in GBM patients, with crucial roles in immunophenotype determination, stromal activation, and carcinogenic pathway activation. We confirmed that HSPA7 promoted macrophage infiltration and SPP1 expression via upregulating the YAP1 and LOX expression of glioblastoma stem cells (GSCs) in vitro and in our clinical GBM tumor samples. We also confirmed that knockdown of HSPA7 might increase the efficiency of anti-PD1 therapy utilizing the GBO model, highlighting its potential as a novel target for immunotherapy. Conclusions: Our results indicated that HSPA7 could be a novel immunotherapy target for GBM patients.

Introduction

Glioblastoma (GBM), one of the most aggressive brain tumors, currently has no effective and sufficient therapies due to its intratumoral heterogeneity and molecular complexity. Immune checkpoint blockade (ICB) therapy is being actively pursued as a promising treatment option for GBM. However, very few patients respond to this therapy (1–5) partly because of the contribution of prominent immunosuppressive factors in the brain tumor microenvironment (TME) in GBM, including tumor-associated macrophages (TAMs), neutrophils, and regulatory T cells (Tregs) (6–8). The development of biomarkers and identification of the definite molecular mechanism underlying resistance to ICB therapy are urgently needed to identify effective therapeutic strategies for GBM.

N6-Methyladenosine (m6A), the most abundant reversible methylation modification of mRNA, critically affects processes in mRNA metabolism, including splicing, export, translation, and decay. Dysregulation of this modification is clearly linked to diverse pathological processes and disease progression (9, 10), including GBM tumorigenesis (11–15). Recent studies have described the role of m6A modification in regulating the immune response (16–19), prompting us to reveal the importance of the spectrum of m6A-regulated genes and m6A regulatory mechanisms in shaping the TME. Numerous studies have demonstrated that m6A is also present in numerous long non-coding RNAs (15, 20–22) (lncRNAs, transcripts longer than 200 nucleotides but lacking functional coding capacity), and lncRNA m6A modification has emerged as a fundamental player in cancer progression and immune regulation, suggesting a potential association between the tumor immune response and m6A lncRNA modification. However, the lncRNA m6A methylation profile has not been systematically clarified in GBM tumors. Additionally, none of these studies have specifically investigated the roles of m6A-modified lncRNAs in the overall TME landscape in GBM. Therefore, it is worthwhile to obtain comprehensive knowledge of the cellular TME infiltration characteristics mediated by m6A-modified lncRNAs, as this knowledge would contribute to our understanding of the role of m6A modification in immune regulation and guide the development of more effective immunotherapeutic strategies.

Here, we depicted the first transcriptome-wide lncRNA m6A methylation profile in GBM and normal brain tissues via m6A sequencing (m6A-seq) data analysis, revealing the highly distinct lncRNA m6A modification patterns between these two groups. Key immune-stromal-related lncRNAs were identified by differentially expressed gene (DEG) analysis in primary GBM cohorts from The Cancer Genome Atlas (TCGA). Integrating the m6A-regulated lncRNAs revealed in this study, we identified HSPA7 as a novel prognostic factor in GBM patients. Through detailed bioinformatic analyses of HSPA7, we identified its crucial role in immunophenotype determination, stromal activation and carcinogenic pathway activation and highlighted its robust capacity to predict the ICB response. We confirmed that HSPA7 facilitated macrophage infiltration via the YAP1–LOX axis in vitro. We also confirmed that HSPA7 enhanced the efficiency of anti-PD1 therapy utilizing GBM patient-derived glioblastoma organoids, an ex vivo model. These results demonstrated that HSPA7 could be a novel immunotherapy target for GBM patients.

Materials and Methods

Patient Specimens and Public Patient Cohorts

Human GBM and normal brain tissues for m6A-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 were obtained from the cBioPortal database. The m6A-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.

Estimation of TME Cell Characterization

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). ( ) 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.

Functional Annotation and Pathway Enrichment Analysis

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 ( ) 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#/).

Cell Lines and Reagents

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% CO2 at 37°C.

Western Blot

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).

Antisense Oligonucleotides, Lentivirus Transfection

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′.

Immunofluorescence

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)

Fluorescence In Situ Hybridization

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 Model

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% CO2, and 90% humidity sterile incubator.

Transwell Assay

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.

Statistical Analysis

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)

Results

Overview of Transcriptome-Wide m6A Methylation in lncRNAs in GBM and Normal Brain Tissues

To understand the pattern of the lncRNA m6A modification profiles in GBM, three human GBM tumor tissues and three normal brain tissues were subjected to transcriptome m6A-seq, which revealed that a considerable proportion of lncRNAs were extensively m6A-modified. In addition, 2,113 m6A 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 m6A peaks were identified, corresponding to the transcripts of 1,508 genes, namely, 445 lincRNAs, 771 antisense lncRNAs, 283 pseudogenes, and 9 other genes ( ). Further analysis showed that most of the lncRNAs (74.5% of the m6A-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 ( ). Upon further analysis of the distribution profiles of m6A peaks within lncRNAs, we found that m6A 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 1

Overview of m6A methylation within lncRNAs in GBM and normal brain tissues. (A) Numbers and gene types of lncRNAs identified in normal brain and GBM tissue samples as identified by m6A-seq. (B) The number of m6A-modified peaks per lncRNA transcript. (C) Distribution of m6A peaks across lncRNA transcripts. (D) Numbers of common and tissue-specific m6A lncRNAs in normal and GBM brain tissues. (E) Dot plot of Log2FC (lncRNA expression) versus Log2FC (differential m6A methylation) values showing a positive correlation between the overall m6A methylation level and the lncRNA expression level (Pearson’s r = 0.495; P < 0.0001) and the distribution of genes with significant changes in both the m6A (FC ≥ 1.2, P ≤ 0.05) and corresponding lncRNA expression levels in GBM samples compared with normal brain tissues (FC ≥ 2, padj ≤ 0.05).

To reveal the significance of m6A-methylated lncRNAs in GBM, the differences and overlaps in the m6A-modified lncRNAs between the GBM and normal brain groups were analyzed by constructing a Venn diagram. As shown in , 781 m6A-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 m6A modification patterns between the GBM and normal groups. To explore the effect of m6A 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 ( ). In addition, the global abundance of m6A peaks between GBM and normal brain tissues was also compared. Furthermore, integrated analysis of these differentially m6A-modified lncRNAs and differentially regulated lncRNAs from the TCGA dataset was conducted. The distribution of genes with a significant change in both the m6A level (|FC| ≥ 1.2, P ≤ 0.05) and the overall transcript expression level (|FC| ≥ 2, padj ≤ 0.05) is shown in . 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 m6A peaks and the expression levels of their corresponding genes [ , Pearson correlation coefficient (r) = 0.495, P<0.0001], and the m6A modification sites in all of the above four groups of genes were distributed in exonic regions, reconfirming that m6A modification can regulate the expression of mature lncRNAs. These results revealed obviously distinct m6A modification patterns between GBM and normal brain tissues. Moreover, m6A modification could regulate numerous lncRNAs, possibly by regulating their stability, degradation or other functions, as reported previously.

Identification of the Immune-Stromal-m6A-Related Pseudogene HSPA7 as a Novel Prognostic Risk Factor in GBM

To investigate the effects of these m6A-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, ). 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 ( and ). Simultaneously, 734 upregulated lncRNAs and 269 downregulated lncRNAs (|FC| ≥ 1.5, padj ≤ 0.05) were identified based on the differential analysis of stromal scores ( and ). As the Venn diagram ( ) indicates, four identical upregulated genes and 11 identical downregulated genes were related to immune activation, stromal activation, and m6A modification. Then, we performed Kaplan–Meier analysis on patients stratified by the expression levels of these 15 differentially expressed genes, and only HSPA7 ( ) and AC011899.9 ( ) 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 ( ), while AC011899.9 was not ( ). Thus, we focused only on HSPA7 in the following work.

Figure 2

Identification of the immune-stromal-m6A-related pseudogene HSPA7 as a novel prognostic factor in GBM. (A) The GEPIA database showed that HSAP7 was overexpressed significantly in GBM tissues compared with GETx normal brain tissues. (B) Kaplan–Meier survival curves show that HSPA7 is a prognostic risk factor in GBM. (C) qPCR assays showing that HSPA7 expression was significantly higher in GBM cells than in NHA cells. Data are presented as the mean ± SD, n = 3. Means were compared with one-way ANOVA, and the NC group is indicated as the control. (D) The gene IGV plots of HSPA7 in the m6A-seq. (E) MeRIP assay showing that HSPA7 was highly enriched by the m6A antibody, and the modification can be regulated by the m6A methyltransferase METTL3. Data represent mean ± S.D. from three independent experiments. qPCR assay showing that the expression of HSPA7 can be regulated by (F) m6A methyltransferase METTL3 and (G) demethylase FTO. Data represent mean ± S.D. from three independent experiments. (H) Distribution of HSPA7 expression in different cohorts stratified by molecular subtype (CL, n = 56; MES, n = 51; PN, n = 44; PN vs MES, P < 0.0001; CL vs MES, P = 0.0020; PN vs CL, P = 0.0013), IDH1 status (mutant, n = 8; WT, n = 141; P = 0.0002), G-CIMP status (G-CIMP, n = 9; non-G-CIMP, n = 142; P < 0.0001), MGMT promoter status (methylated, n = 55; unmethylated, n = 67; P = 0.7445), age (high, age > 55, n = 100; low, age ≤ 55, n = 53; P = 0.5114), and sex (female, n = 54; male, n = 99; P = 0.3996). (I) Univariate and multivariate Cox regression analyses of HSPA7 and other clinical features in the overall survival of GBM samples. The statistical significance is shown as: ns, P > 0.05; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

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 ( ). The m6A methylation peak distribution and abundance in HSPA7 transcripts from GBM and normal brain tissues, as detected by m6A-seq, were visualized using IGV software ( ). We found that HSPA7 was highly enriched in the m6A-precipitated fraction, and the m6A modification enrichment level could be regulated by methyltransferase-like 3 (METTL3), which contains a catalytic activity domain to catalyze m6A formation ( ). Additionally, the expression of HSPA7 was significantly inhibited by knocking down METTL3 ( ) and overexpressing FTO, an m6A demethylase ( ). These results provided evidence that HSPA7 harbored high m6A modification levels and thus could be regulated in an m6A-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 , 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 ( ). These results indicated that the pseudogene HSPA7 is a novel risk prognostic biomarker and indicates therapeutic outcomes of GBM patients.

HSPA7 Is Correlated With Immunophenotypes and TME Landscapes

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 , 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 ( ), 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 ( and ), 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 ( ). 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 ( ), and positive correlations were found between HSPA7 expression and these molecules ( ). 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) ( ), 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.

Figure 3

HSPA7 is correlated with immunophenotypes and TME landscapes. (A) The enrichment scores of immune cell types and immune-related function-related gene sets were calculated via the ssGSEA algorithm. A heatmap was used to visualize these immune characteristics between the HSAP7 high and low expression groups; yellow represents a high enrichment level, black represents a median enrichment level, and blue represents a low enrichment level. (B) High expression of HSPA7 presented significantly increased stromal, immune, and ESTIMATE scores and decreased tumorigenicity. (C) HSPA7 positively correlated with immunosuppression regulators (checkpoints, macrophages, neutrophils, and Tregs). High expression of HSPA7 presented significantly increased markers of (D) suppressor function of myeloid lineages, (E) immune inhibitory checkpoints and (F) major neutrophil-recruiting chemokines and their receptors. The asterisks indicate a significant statistical p-value calculated using the non-parametric Wilcoxon test (*P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001). (G) Compared to the immune-L group (102 patients), patients in the immune-H group (51 patients) experienced poorer prognoses. (H) HSPA7 expression was higher in the immune-H group than in the immune-L group. (I) The proportion of GBM molecular subtypes in the low and high HSPA7 groups. MES subtype, red; CL subtype, orange; PN subtype, blue. (J) Alluvial diagram showing the changes in immune phenotypes, GBM molecular subtypes, and HSPA7 expression. (K) HSPA7 positively correlated with CD44, a marker of the MES subtype. (L) GSEA of mesenchymal signatures showed that GBM samples with high HSPA7 expression were enriched in the MES subtype compared to GBM samples with low HSPA7 expression. NES, normalized enrichment score; FDR, false discovery rate. (M) qPCR assays showed that HSPA7 expression was higher in MES subtype GSCs (GSCs 20 and 267) than in PN subtype GSCs (GSCs 8–11) and neural stem cells (NPCs). (N) qPCR assays verified the knockdown efficiency of HSPA7. (O) Western blot assays showed that HSPA7 promoted the expression of CD44. The statistical significance is shown as: ns, P > 0.05; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

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 ( ). Patients in the immune-H group (n = 51) exhibited poorer prognosis ( ) and higher HSPA7 expression levels ( ) 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 ( ). 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 ( ). 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 , 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 ( ). GSEA also showed that GBM samples with high HSPA7 expression were enriched in the MES subtype compared to GBM samples with low HSPA7 expression ( ). In addition, we found that HSPA7 expression was higher in MES GSCs than in PN GSCs ( ). Moreover, Western blot assays revealed that knockdown of HSPA7 reduced the expression of CD44 ( ). 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.

HSPA7 Is Correlated With Stromal and Carcinogenic Activation Pathways

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 ( ). 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 , 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 ( ), as constructed by Mariathasan et al. (28), positively correlated with fibroblast enrichment ( ), as calculated by the MCP-counter method (38), and positively correlated with typical stromal cell activation-related pathways ( ). 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 ( ). 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.

Figure 4

HSPA7 is correlated with stromal and carcinogenic activation pathways. (A) GSVA enrichment analysis showing the activation status of biological pathways in the HSPA7-high and HSPA7-low groups. A heatmap was used to visualize these biological processes. Yellow represents activated pathways, black represents moderately activated pathways, and blue represents inhibited pathways. (B) Correlations between each known gene signature in the TCGA GBM cohort using Pearson correlation analysis. Negative correlations are marked with blue, and positive correlations are marked with orange. (C) Differences in stromal activation-related pathways between the HSPA7-high and HSPA7-low groups. EMT, epithelial–mesenchymal transition; Pan-FTBRS, panfibroblast TGF-β response signature. (D) HSPA7 expression was positively correlated with fibroblast enrichment, as calculated by MCP-counter. (E) HSPA7 expression was positively correlated with stromal activation-related pathways. (F) A heatmap was used to visualize the expression of stromal activation-related genes. Yellow represents high expression, black represents the median expression, and blue represents low expression. The statistical significance is shown as: **P < 0.01; ****P < 0.0001.

Then, we explored the pathways enriched with genes positively correlated with HSPA7 (Pearson’s r ≥ 0.3, P ≤ 0.05, ) 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 ( ). 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 ( ), 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 ( ). 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 5

The genes that were positively correlated with HSPA7 were enriched in immune response- and stromal-related pathways. (A) Bar graph of enriched terms, colored by p-values, and corresponding network of enriched terms. (B) Summary of enrichment analysis in PaGenBase and (C) summary of enrichment analysis in TRRUST across HSPA7 positively correlated genes. (D) GSEA analyses displayed key immune-regulated pathways enriched in the high (up) and low (down) HSPA7 groups, both in the GO biological process (left) and KEGG datasets (right), and each line is for one pathway.

Verification of the Functions of HSPA7 in Two CGGA Cohorts

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.

M6A-Modified HSPA7 Is Potentially Regulated by the Methyltransferase WTAP

To explore the potential molecular mechanism regulating HSPA7 RNA metabolism, we first screened HSPA7-interacting proteins via the NCBI database ( ). We verified that WTAP (45), a regulatory subunit of the m6A 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 , 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 ( ). 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 ( ). 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 ( ), CGGA GBM cohort 1 ( ) and CGGA GBM cohort 2 ( ). These patterns were the same as the patterns observed for HSPA7 expression. The above findings indicate that m6A 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 ( , ). 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 ( and ). However, the exact mechanism by which HSPA7 regulates the GBM TME requires further clarification.

HSPA7 Holds Promise for Predicting the Therapeutic Response to ICB Therapy

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 [ , 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 ( ). 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 [ , 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 ( ).

Figure 6

HSPA7 holds promise for predicting the therapeutic response to ICB. (A) The Kaplan–Meier survival curves showed that HSPA7 was a prognostic risk factor in the melanoma cohort that received anti-CTLA4 therapy. (B) Expression of HSPA7 in distinct anti-CTLA4 clinical response groups, R, response; NR, no response. (C) The proportion of patients who responded to CTLA4 blockade immunotherapy in the low or high HSPA7 expression groups. (D) The Kaplan–Meier survival curves showed that HSPA7 is a prognostic risk factor in the anti-PD-L1 cohort (IMvigor210). (E) Expression of HSPA7 in distinct anti-PD-L1 clinical response groups, R, response; NR, no response. (F) The proportion of patients who responded to PD-L1 blockade immunotherapy in the low or high HSPA7 expression groups. PD-L1 expression in both ICs (G) and TCs (H) was associated with HSPA7 expression, with the highest PD-L1 expression level in cells showing the highest HSPA7 expression. (ANOVA, p = 0.0001, p = 0.0009, respectively). HSAP7 expression between (I) FBXW7 wild-type (WT) and mutant status and (J) FGFR1 wild-type (WT) and mutant status. (K) The waterfall plot of tumor somatic mutations established by those with high HSPA7 expression (right) and low HSPA7 expression (left). The statistical significance is shown as: ns, P > 0.05; *P < 0.05; ***P < 0.001; ****P < 0.0001.

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 ( ), higher stromal activation, higher oncogenic pathway activation, and lower MMR pathway activation ( ), completely consistent with the functional enrichment patterns in GBM ( , and , ). 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 , 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 ( ). 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 ( ). Collectively, our work strongly indicates that HSPA7 expression contributes to predicting the response to immune checkpoint therapy.

HSPA7 Facilitated Macrophage Infiltration and Could Be a Potential Immunotherapy Target for GBM Patients

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 ( ; 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 ( ; 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 ( ). We also confirmed that LOX and YAP1 were higher in MES GSCs than in PN GSCs ( ). Moreover, LOX and YAP1 expression levels were enhanced upon overexpression of HSPA7 in PN subtype GSCs 8–11 ( ) but decreased significantly upon knockdown of HSPA7 in MES subtype GSCs 267 ( ), 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 ( ). Then, THP1-differentiated macrophages were cultured in conditioned medium collected from GSCs 267 transfected with HSPA7 ASO or NC. As shown in , 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 7

HSPA7 facilitated macrophage infiltration and might be a novel immunotherapy target for GBM patients. (A) The heatmap shows that HSPA7 activated YAP1 signaling. The asterisks indicate a significant statistical p-value calculated using the non-parametric Wilcoxon test (*P < 0.05; **P < 0.01; ***P < 0.001). (B) The protein expression of CD44, YAP1, LOX, and YKL40 in PN subtype GSCs (GSCs 8–11) and MES subtype GSCs (GSCs 20 and 267). (C) qPCR assays verified the overexpression efficiency of HSPA7. (D) Representative Western blot assays showed that HSPA7 promoted the expression of YAP1 and LOX, a macrophage chemoattractant (left), and quantification histogram represented the relative protein expression of LOX and YAP1 (right); data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001. Means were compared with Student’s t-test for two groups and one-way ANOVA for three groups, and the NC group is indicated as the control. (E) Representative transwell migration assays showed that HSPA7 inhibited human THP1-differentiated macrophage migration by exposing them to conditioned medium from GSC 267 cells transfected with NC, HSPA7 ASO1, or ASO2 for 24 h. Original magnification, ×100 (left), and quantification histogram represented relative migration of THP-1 macrophages (right); data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001. Means were compared with one-way ANOVA, and the NC group is indicated as the control. (F) Western blot assays showed that HSPA7 inhibited the expression of SPP1 in THP1-differentiated macrophages by exposing them to conditioned medium from GSC 267 cells transfected with NC, HSPA7 ASO1, or ASO2 for 48 h. (G) IF staining in a human GBM tissue microarray showed that the expression of LOX, YAP1, PD-L1, CD68, and SPP1 was higher in the HSPA7 high group than in the low group. Histogram representing statistical proportion data of positive area; data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001, means were compared with Student’s t-test. (H) GBOs at 2 weeks were cocultured with conditioned medium from GSC 267 cells, transfected with NC, HSPA7 ASO1, or ASO2 at 48 h, and PD1 antibody (5 μM) for 5 days as indicated. IF staining for Ki67 and CD44 in GBO sections showed that knockdown of HSPA7 enhanced the effect of anti-PD1 therapy, original magnification, ×630 (scale bars: 20 μm). Histogram representing statistical proportion data of positive area; data are presented as the mean ± SD, n = 3. *p < 0.05, **p < 0.01, ***p < 0.001. Means were compared with one-way ANOVA, and the NC group is indicated as the control.

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 ( ). We further explored the association between the HSPA7 expression level and TMB. As shown in , 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 ( ), which was negatively correlated with immunosuppression, stromal activation, and oncogene pathway activation ( ). 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 ( ). 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 ( ), 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 ( ) 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 , 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.

Characterization of HSPA7 Across 33 Cancer Types

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 , 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 - . 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 ( ) and stromal score ( ) in all 33 cancer types. Moreover, HSPA7 was significantly positively correlated with immune checkpoint expression ( ) and the immune response ( ) 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.

Discussion

GBM, one of the most aggressive brain tumors, currently has no effective and sufficient therapies due to its intratumoral heterogeneity and molecular plasticity. ICB therapy is being actively pursued as a promising treatment option for GBM, but very few patients respond to ICB therapy. Thus, identifying markers regulating the brain TME that are prominent players in ICB therapy for cancer could reveal promising new targets for therapeutic intervention. LncRNAs and m6A modifications are emerging as indispensable regulators of the TME. However, the overall TME infiltration characteristics mediated by m6A-modified lncRNAs have not been comprehensively recognized. Therefore, it is worth obtaining comprehensive knowledge of the cellular TME infiltration characteristics mediated by m6A-regulated lncRNAs.

Here, based on our m6A-seq data, we revealed highly distinct lncRNA m6A methylation modification patterns between GBM and normal brain tissues. In addition, we identified immune-stromal-m6A-related HSPA7 as a novel prognostic risk factor in GBM patients; this gene plays a crucial role in immunophenotype determination, stromal activation, and oncogenic pathway activation and has a robust capacity to predict the ICB response. Furthermore, in this study, WTAP, a methyltransferase mediating m6A modification, was identified as a potential regulator of HSPA7 expression. Further analysis showed that WTAP, like HSPA7, can also regulate immunophenotype determination and stromal activation-related pathways ( , ). This finding again demonstrated that m6A modification is highly important in shaping the TME landscape. Many studies have reported that stromal cell activation in the TME can suppress immune infiltration or facilitate an immunosuppressive response in the TME, mediating therapeutic resistance to ICB (28, 58). Our data also revealed a markedly negative correlation between HSPA7 expression and TMB, a marker of the response to ICB therapy, in a TCGA GBM cohort ( ). In addition, another study showed that PTEN mutations were associated with immunosuppressive expression signatures in non-responders to anti-PD-1 immunotherapy in GBM (2). Chen et al. (52) found that PTEN deficiency in GBM increases the infiltration of SPP1+ macrophages, which can interact with fibroblasts and vascular endothelial cells, inducing angiogenesis, EMT, and some other stromal activation-related pathways in colon cancer and are characterized by expression of the pattern recognition receptor MARCO (59) via a YAP1–LOX-β1 integrin–PYK2 axis. Intriguingly, antibody-mediated depletion of MARCO inhibits cancer progression and metastasis, enhancing ICB efficacy (60), suggesting that SPP1+ macrophages play a non-negligible role in the immunorepressive response and immunotherapeutic resistance. 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) ( ); we also confirmed that HSPA7 could facilitate the macrophage infiltration into the GBM TME via YAP1–LOX axis in vitro, as well as into our local GBM tissue sections ( ). Furthermore, we found that HSPA7 can interact with YAP1 (61), which can induce the secretion of CCL2/CSF1 to recruit monocytes (62), suggesting that HSPA7 may be a target that synergistically regulates SPP1+ macrophages, which then induce stromal activation in the TME. Moreover, we found that HSPA7 and YAP are also contained in GBM extracellular exosomes (63), a fundamental regulator of TME cells. Numerous research reports have indicated that YAP is a hub in the network of signals exchanged within the TME (64), thus regulating the immune response and stromal activation. However, the specific mechanism of HSPA7 needs to be proven by a large number of experiments. We also confirmed the predictive value of HSPA7 in a cohort of melanoma patients who received anti-CTLA4 therapy and an IMvigor210 cohort of patients treated with anti-PD-L1 therapy. A statistically significant difference in HSPA7 expression was found between non-responders and responders. MMR deficiency has recently emerged as a beneficial indicator of the response to PD-1 blockade in patients with cancer (65, 66), whereas we found that HSPA7 was negatively correlated with MMR-related pathway activity ( ), which in turn was negatively correlated with stromal activation and oncogenic pathway activation ( and ). This pattern suggests that MMR deficiency is an unfavorable factor for the response to ICB therapy. Other studies have shown the same results (28, 55). By comprehensive analysis of data for 10,294 samples, a recent study showed that MMR-deficient gliomas are 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 (54), reconfirming the ability of HSPA7 to predict the response to immunotherapy. We also explored the function of HSPA7 in various cancer types. As shown in - , HSPA7 expression had prognostic significance in most cancers; in addition, it can regulate the expression of immune checkpoint genes and the activity of immune response pathways and was positively correlated with the immune score and stromal score in the majority of tumor types. HSPA7 was found to be a risk factor in the cohort of melanoma patients who received anti-CTLA4 therapy ( ) but was found to be a beneficial prognostic factor in the TCGA SKCM cohort (23) ( , ), exhibiting no obvious correlation with immune checkpoint expression ( ) or the immune response ( ). The mechanisms underlying this discrepancy need future exploration.

In summary, via integrated analysis of our own m6A-seq data and public clinical data, we found highly distinct lncRNA m6A methylation modification patterns between GBM and normal brain tissues and determined that HSPA7 could be used to evaluate the prognosis of GBM patients, as well as the corresponding cellular TME infiltration and activation characteristics of individual patients, and could predict the clinical efficacy of anti-PD-1/-PD-L1 immunotherapy in patients. More importantly, this study offers several novel insights regarding cancer immunotherapy; for example, targeting m6A regulators to change the m6A modification patterns of key immune-regulating targets may reverse unfavorable cellular TME infiltration characteristics. This study contributes to future exploration of novel drug combination strategies or novel immunotherapeutic agents.

Data Availability Statement

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.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethical Committee on Scientific Research of Shandong University Qilu Hospital (approval number: KYLL-2018-324). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

GL and HX supervised the project. RRZ designed the research and performed all experiments. BYL completed the basic experiment part. ZH performed statistical analysis with the R language. SJZ, QDG, PZ, WQ, SBW and ZHC were responsible for clinical sample collection and subsequent sample delivery. XG, YHQ, ZWP helped to revise the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the National Natural Science Foundation of China (Nos. 81874083; 82072776; 82072775; 81702468; 81802966; 81902540; 81874082; 81472353), Natural Science Foundation of Shandong Province of China (Nos. ZR2019BH057; ZR2020QH174), Key clinical Research project of Clinical Research Center of Shandong University (2020SDUCRCA011) and Taishan Pandeng Scholar Program of Shandong Province (No. tspd20210322).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.