An Immune Microenvironment-Associated Gene Signature Predicts Outcomes and Therapeutic Response in Pediatric Medulloblastoma

doi:10.21203/rs.3.rs-2723037/v3

Download PDF

Research Article

An Immune Microenvironment-Associated Gene Signature Predicts Outcomes and Therapeutic Response in Pediatric Medulloblastoma

https://doi.org/10.21203/rs.3.rs-2723037/v3

This work is licensed under a CC BY 4.0 License

Version 3

posted

You are reading this latest preprint version

Purpose: The tumor microenvironment (TME) is a pivotal factor in the progression and treatment response of cancers, including medulloblastoma (MB), the most common malignant brain tumor in children. This study investigates TME-specific gene signatures to create a prognostic model for MB patient outcomes.

Methods: We sequenced 240 MB patient samples at Beijing Tiantan Hospital with RNASeq and analyzed TME components. Through co-expression network analysis and regression models, we identified an eight-gene signature that correlated with TME elements. This signature was tested against an independent dataset (GSE85217) for validation.

Results: Our study presents the TME-associated risk score (TMErisk), which incorporates the eight-gene (CEBPB, OLFML2B, GGTA1, GZMA, TCIM, OLFML3, NAT1, and CD1C) signature. Higher TMErisk scores were associated with decreased overall survival and lower immune cell infiltration, immune checkpoint activity, and human leukocyte antigen expression. There was also a notable negative correlation between TMErisk scores and both TMB and IPS, varying across MB molecular subtypes. Moreover, the TME-risk was inversely related to the tumor mRNAsi, implying the influence of the TME on tumor stemness.

Conclusions:Our findings identify a TME-specific eight-gene prognostic model that may serve as a predictive biomarker for MB patient outcomes and responses to immunotherapy. This gene signature model offers a supplementary tool to current WHO molecular subtypes and provides a potential target for future TME-focused MB treatment strategies. However, the distinct TME landscapes across MB subtypes pose new challenges for further research.

Keywords: gene signature, immune cells, medulloblastoma, tumor microenvironment, prognostic model, tumor, stromal cells

Oncology

Bioinformatics

gene signature

immune cells

medulloblastoma

tumor microenvironment

prognostic model

tumor

stromal cells

Medulloblastoma (MB), the predominant malignancy of the central nervous system in the pediatric population, segregates into four principal molecular subgroups—WNT, SHH, Group3, and Group4 (Louis et al., 2016; Northcott et al., 2017). These subgroups exhibit unique clinicopathological features, transcriptional signatures, and tumoral microenvironments (Bockmayr et al., 2018). The heterogeneity within tumors complicates the prognosis for MB patients, necessitating personalized assessments.

The tumor microenvironment (TME)—the milieu where tumors originate, grow, and spread—encompasses a complex array of non-malignant stromal elements and cancerous cells. Acknowledged for its role in oncogenesis, the TME influences cancer immunophenotypes and patient prognoses (Chen et al., 2015; Fridman et al., 2012). Advances in immunotherapy highlight the potential of leveraging the TME for therapeutic interventions and prognostic assessments.

Accurate characterization of the TME's intricate constituents, crucial for prognostication through tumor immunity, presents a formidable challenge. The ESTIMATE algorithm emerges as a potent tool for analyzing immune and stromal cell populations (Yoshihara et al., 2013), thus furnishing insights into TME dynamics. Leveraging a comprehensive database of MB bulk RNA sequences from Beijing Tiantan Hospital Capital Medical University, this investigation has unearthed gene signatures indicative of TME status and prognostic outcomes in MB cases. Subsequently, researchers have developed a TME-related risk score model—TMErisk—to forecast survival rates in MB patients.

Sample sources and analysis

The bulk mRNA expression data of 322 MB patients from Beijing Tiantan Hospital were used as the training cohort (240 patients with complete follow-up data available). Data on sex, age, molecular subgroup, survival status, and survival time were retrospectively collected. All participants or their parents signed an informed consent form for our study. The validation gene expression dataset was obtained from GEO (GSE85217), and clinical data were obtained from published articles (Cavalli et al., 2017) (Table S1).

Our study conformed to the tenets of the Declaration of Helsinki and was approved by the Ethics Committee of Beijing Tiantan Hospital, Capital Medical University. (KY 2019-098-01)

DNA and RNA sample collection and whole genome sequencing

For whole genome sequencing (WGS), the tumor samples of all participants were collected, and genomic samples were extracted using a FlexiGene DNA Kit (Qiagen, Venlo, The Netherlands) and TRIzol RNA Kit. Genomic samples were stored in EDTA tubes at -20 degrees Celsius according to the manufacturer’s protocol. The specimens were sent to the BGI for total RNA sequencing using the BGISEQ-50 platform (BGI, Shenzhen, China). Bioinformatics analysis was performed using RNA-seq (FPKM normalized) and clinical phenotypic data from MB patients. WGS was performed on samples that met high-throughput sequencing standards, and the reads were sequenced on the MGISEQ-2000 platform (DNBSEQ-G400). The size of the reads generated during paired-end sequencing was 150 bp, and the average depth was 60 ×, while the median depth was 44 ×. Reads were mapped to the HG38 reference genome, and variants were called using the Burrows–Wheeler aligner and annotated with an in-house annotation pipeline.

TMErisk scores and pan-cancer analysis based on gene signature

The TME of MB was evaluated by immune and stromal scores, and the "estimate" R package provided details of the algorithm (Yoshihara et al., 2013). The correlation between the TME and TME scores was evaluated using the Spearman method. The "maxstat" R package was used to calculate the optimal cut-off point for continuous variables (Table S2).

Differentially expressed genes (DEGs) in the high- and low-risk groups were screened using the Wilcoxon method. Data with an average gene expression > 0 were filtered. Genes with P-values (false discovery rate, FDR) adjusted by the Benjamin–Hochberg method < 0.05 and Log2 Fold-change (FC) > 1 were defined as DEGs. Weighted gene co-expression network analysis (WGCNA) was used to select gene modules associated with MB immunity and stromal cell scores, and correlation coefficients greater than 0.5 were considered to be strongly correlated with TME. The gene expression matrix was transformed to log2 (FPKM+1). The DEGs and WGCNA gene intersections were shown in a Venn diagram. The least absolute shrinkage and selection operator (LASSO) algorithm was used to screen continuous variables (Tibshirani, 1997). The results of LASSO were transformed into binary variables, and the final TMErisk prognosis model was screened based on the Akaike information criterion (AIC), using Cox stepwise backward regression. The MB risk-scoring formula was calculated as follows: TMErisk = ∑ Expi*Coefi. A flowchart of the TMErisk score development process is presented in Figure 1. UCSCXenaShiny software was used to analyze the pan-cancer data of gene signatures in the Cancer Genome Atlas (TGCA).

Analysis of immune microenvironment and MB TMErisk

The correlation between the MB immune microenvironment and TMErisk includes two parts: (1) correlation and differences between immune checkpoint genes and human leukocyte antigen (HLA) family genes in TMErisk. (2) Differences in immune cell and stromal cell infiltration among the different risk groups. We used four algorithms to analyze immune cell infiltration: CIBERSORT (Newman et al., 2015), xCell (Aran et al., 2017), ssGSEA (Barbie et al., 2009) and MCPcounter (Becht et al., 2016). "ClusterProfiler" packages were used for gene enrichment analysis and enrichment in the different risk groups (Yu et al., 2012). Cell function and pathways were annotated using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG).

Gene expression-based stemness index, immune-phenoscore, and drug-sensitivity prediction

Malta et al. used one-class logistic regression (OCLR) machine-learning algorithms to train a stemness signature using gene expression in pluripotent stem cells from the Progenitor Cell Biology Consortium dataset (Malta et al., 2018). We implemented OCLR using the "gelnet" package (https://www.rdocumentation.org/packages/gelnet/versions/1.2.1) to determine the gene expression-based stemness index (mRNAsi), and the results were normalized to 0–1. We used the GDSC2 dataset after RMA standardization and log transformation as the training set, and build the prediction model with "oncoPredict" package (Maeser et al., 2021). The sensitivity of all MB patients to 198 drugs was calculated. To assess the reactivity of MB patients, we calculated the immune-phenoscore (IPS) of the sample using the R script provided by the Cancer Immunome Atlas (https://tcia.at) (Charoentong et al., 2017; Hakimi et al., 2016).

Tumor mutations in MB patients

The quality of the sequencing data obtained from tumor tissue samples was tested, and a database was built for samples that met the conditions. WGS was used to analyze the tumor mutation burden (TMB) and mutation status of the eight gene signatures. Mutation data files containing MB clinical phenotype samples were obtained from the BGI. Mutation data files were analyzed using the "maftools" package.

Other statistical analysis

The Kaplan–Meier survival curve was used to show the overall survival time of MB patients with different risks, and the log-rank test was used to compare K-M curves between groups. The area under the curve was calculated to evaluate the predictive efficacy of the gene signature prognostic model. We performed 1:1 propensity score matching according to the risk groups to balance the baseline levels between the groups. All P values were analyzed using double-tailed tests, and P < 0.05 was defined as statistically significant. R software 4.2.0 (https://www.r-project.org) was used for all statistical analyses.

Clinical characteristics and tumor microenvironment of MB patients

In our comprehensive analysis of medulloblastoma cases, immune and stromal scores were quantitatively assessed for 322 patients at Beijing Tiantan Hospital. The range of stromal scores spanned a broad spectrum from -2043.4 to 1152.5, whereas immune scores varied between -2022 to 1977. A significant elevation in immune scores characterized Group3 compared to other molecular subgroups, including WNT (P = 0.010), SHH (P < 0.001), and Group4 (P = 0.004), as depicted in Figure S1A. Investigations into demographic variables revealed no notable differences in immune scores attributable to gender (Figure S1B) or across four distinct age brackets (Figure S1C). Nevertheless, stromal scores distinguished Group3 and the SHH subgroup with statistical significance (P = 0.001, Figure S1D), while gender and age-related disparities remained insignificant (Figures S1E and S1F).

Risk classification, derived from immune and stromal cell metrics, enabled the construction of gene expression heatmaps for each subgroup (Figure S1G and S1H). Survival analysis, employing the Kaplan–Meier method, disclosed a favorable prognostic implication for patients with heightened immune scores (log-rank P = 0.027, Figure S1I). Conversely, stromal scores did not manifest a significant prognostic differentiation (log-rank P = 0.14, Figure S1J). The samples under study exhibited over 80% tumor purity, with a pronounced inverse correlation between tumor purity and both immune and stromal scores (Figures S2A and S2B).

Development of TMErisk score model

In our study, we sought to establish a robust prognostic model that hinges on the unique characteristics of the tumor microenvironment (TME) in medulloblastoma (MB). The development of the TMErisk score model involved analyzing differentially expressed genes (DEGs) among groups stratified by immune and stromal scores within the Beijing Tiantan Hospital cohort (Figures 2A,B). Over 80% tumor purity was observed in our sample cohort, and a significant negative correlation was found between both immune and stromal scores and tumor purity (Figures S2A, B).

After selecting a soft threshold of 16 (Figure S2C), we applied weighted gene co-expression network analysis (WGCNA) to identify relevant gene modules. The analysis yielded 27 gene modules (Figure 2C), with dark red and yellow modules showing strong correlations with immune and stromal scores, respectively (Figure 2D). We then pinpointed 281 intersecting genes between the DEGs and WGCNA genes (Figure 2E).

Subsequent least absolute shrinkage and selection operator (LASSO) regression analysis refined this list to 12 pivotal genes (CEBPB, HCST, CSF1R, ADAMTS15, OLFML2B, GGTA1 (also named GGTA1P), RHOD, GZMA, and TCIM (also named C8ORF4), OLFML3, NAT1, and CD1C). Following backward stepwise regression, we arrived at an optimal eight-gene signature that underpins the TMErisk score. This prognostic model, expressed as a linear combination of these gene signatures weighted by their respective coefficients, demonstrated potential as a predictive tool for MB patient survival. The final prognostic model was TMErisk score = 1.08 × CEBPB + 0.65 × OLFML2B − 1.16 × GGTA1 − 0.92 × GZMA + 0.88 × TMIM - 1.45 × OLFML3 + 1.24 × NAT1 - 1.27 × CD1C according to the minimum Akaike’s Information Criterion (AIC).

Prognostic assessment and external validation of TMErisk

We subjected the TMErisk score model to rigorous prognostic assessment and validation using both the development cohort from Beijing Tiantan Hospital and an external dataset (GSE85217). The univariate and multivariate Cox regression analyses, visualized through forest plots (Figure 3A), affirmed the model's predictive strength. Patients with elevated TMErisk scores consistently correlated with poorer outcomes in both cohorts (Figures 3B, C).

The relationship between molecular subgroups and risk stratification was further elucidated via Sankey plots, highlighting the distribution of patients across varying risk levels based on TMErisk scores (Figure 3D). The concordance statistic (C-statistic) was notably higher for the TMErisk score compared to molecular subgroups alone in the development cohort, indicating superior prognostic accuracy (Figure 3E).

Additionally, the model's capacity to predict patient outcomes was supported by the time-dependent area under the receiver operating characteristic (ROC) curve (AUC) values. These values suggested that the TMErisk scores could effectively forecast the prognoses of MB patients in both the Beijing Tiantan Hospital dataset and the GSE85217 dataset (Figure 3F). The proportional hazards assumption underlying Cox regression was satisfactorily met, as confirmed by the Schoenfeld residuals (Figure S3).

The AUC values for 1-year, 3-year, and 5-year overall survival (OS) further validated the risk scoring model's moderate predictive efficiency in the development group (Figure S4A). This efficiency was mirrored in the validation group's ROC curve analysis for the same time points (Figure S4B). Calibration plots for the development and validation cohorts exhibited close alignment between predicted outcomes and actual observations, reinforcing the model's robustness (Figures S4C, D).

Differential analysis of clinical characteristics and risk scores indicated that the risk score for Group3 was significantly elevated compared to other groups (Figures S4E–G). Survival analyses for high- and low-risk MB patients from both Beijing Tiantan Hospital and the external validation cohort (GSE85217) are depicted in Figures S5A–D. Propensity score matching (PSM) was employed to adjust for covariate differences between high-risk and low-risk groups, confirming the consistency of results pre- and post-PSM (Figures S5E and F; Table S3).

Prognosis of gene signatures in MB and pan-cancer

High expression of GGTA1 and OLFML3 resulted in better OS in MB patients in the development and validation groups (Figures S6A, B). However, high expression of TCIM was associated with worse OS. For GZMA and CD1C, the OS results of the development and validation groups were opposite. The results of Cox proportional hazards model analyses were consistent with the OS results (Figures S7A, B).

Pan-cancer analysis also suggested that GGTA1 expression had a protective effect in low-grade gliomas (LGG) (Figure S7C). As in MB patients, high TCIM expression was also a risk factor for LGG (Figure S7C).

TMErisk score is associated with MB immune signature

The results of Gene Set Enrichment Analysis (GSEA) in the high- and low-risk groups are shown using volcano diagrams and are annotated based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (Figures 4A, B; Tables S4, S5). Muscle system process, rRNA metabolic process, and RNA splicing epidermis were enriched in the high-risk group, while epidermis development and the external side of the plasma membrane were enriched in the low-risk group.

KEGG enrichment analysis identified ribosomes, hypertrophic cardiomyopathy, dilated cardiomyopathy, cardiac muscle contraction, and oxidative phosphorylation as enriched in the high-risk group. Olfactory transduction and cell adhesion molecules were enriched in the low-risk group. Figure 4 C–F shows the top-five results for the high- and low-risk groups after GO and KEGG annotation.

We also studied the correlation of 19 HLA family genes and 48 immune checkpoint genes with the identified gene signatures and TMErisk scores. The Spearman test showed that expression of seven HLA family genes (HLA-DRB1, HLA-DRA, HLA-DPA1, HLA-DPB1, HLA-DMA, HLA-DMB, HLA-DOA) and 16 immune checkpoint genes (including CD200, IDO1, LAIR1, TIGIT, PDCD1LG2, and SELPLG) were significantly correlated (Figures 5 A–D, S8). In addition, the Wilcoxon test was used to compare the high- and low-risk groups (Figure S9A, B).

Correlation analysis of immune and stromal scores with TMErisk scores showed a significant negative correlation in the SHH group (Figure S10).

Four algorithms: CIBERSORT, MCP-counter, ssGSEA, and xCell, were used to estimate TMErisk, gene signatures, and immune cell infiltration (Figure 6A–D). We investigated the infiltration of 22 types of immune cells in the different molecular subgroups in the CIBERSORT results, and found significant differences for seven types of immune cells (Figure S11). The immune cell differences between the high- and low-risk groups based on the TMErisk score as estimated by the four algorithms were calculated (Figure S12). CIBERSORT results showed that CD8 T cells were significantly increased in the low-risk group as compared to in the high-risk group (Figure S12A) SsGSEA results showed that activated B cells were significantly increased in the low-risk group as compared to in the high-risk group (Figure S12C). The results of pan-cancer analysis showed that eight gene signatures, such as GGTA1, OLFML3, GZMA, and CD1C, were significantly correlated with most immune cells (Figure S13).

Immuno-phenoscore, mRNAsi, and drug sensitivity

Our analysis revealed that the immuno-phenoscore (IPS), which evaluates the likelihood of response to immunotherapy, was significantly elevated in Group3 compared to other MB molecular subgroups, showcasing a negative correlation with the TMErisk score (Figures S17A, C). This relationship emphasizes the potential of TMErisk to reflect immunotherapy responsiveness. The correlation circle diagram, which included IPS, gene expression stemness index (mRNAsi), immune score, stromal score, TMErisk score, and gene signature, indicated a significant association between IPS and gene expression, except for CD1C (Figures S17B, D).

The study also examined the gene expression-based stemness index (mRNAsi) across different molecular subtypes, finding that Group4 possessed significantly higher mRNAsi scores than other groups (Figure S18A). Wilcoxon’s test revealed a higher mRNAsi in the low-risk group compared to the high-risk group, underscoring a significant negative correlation between TMErisk and mRNAsi (Figures S18B, C). Further analysis of the association between gene signatures and mRNAsi identified six genes with significant correlations (Figure S18D).

Moreover, we explored the sensitivity to various therapeutic drugs among different TMErisk-defined risk groups. The results suggested distinct drug response profiles, which may inform therapeutic decision-making (Figures S14, S18E).

Mutation status in the signature genes of MB patients

In our cohort, sequencing and subsequent analysis verified that out of 94 samples suitable for high-throughput sequencing, 91 (96.81%) exhibited mutations in the eight signature genes derived from the TMErisk model. Most notably, GZMA and TCIM mutations were present in 91% of the samples, with OLFML2B mutations identified in 83%. Predominantly, these mutations were of the missense type. Single-nucleotide polymorphism analysis highlighted that C to T transitions constituted the majority of these alterations (Figure S15A, B).

Tumor mutational burden (TMB) has gained recognition as a biomarker that can predict immune response to tumors (Yarchoan et al., 2017). Our findings indicate that in the Group4 molecular subgroup, there was a notable negative correlation between TMB and the TMErisk score (Figure S16), suggesting a potential link between mutation load and the TME.

The present study developed an eight-gene signature prognostic model grounded in the immune microenvironment utilizing 322 MB patients from Beijing Tiantan Hospital and externally validated it employing a public dataset (GSE85217). Our findings suggest the TMErisk score may function as an independent predictor of survival outcomes in MB patients. The TMErisk score and gene signature exhibited significant correlation with HLA family gene expression, immune checkpoint gene expression, immune-related cell infiltration, IPS, and mRNAsi. Additionally, somatic mutations occurred in five of the eight signature genes (Figure S15A). Given the varying drug sensitivities across TMErisk groups, TMErisk may serve as a potential biomarker for predicting immunotherapy response in MB patients.

The ESTIMATE algorithm's efficacy in prognostic exploration has been demonstrated in glioblastoma (Jia et al., 2018), clear cell renal cancer (Xu et al., 2019) and gastric cancer (H. Wang et al., 2019). Unlike other cancers, MB patients with high immune and stromal cell scores based on ESTIMATE have superior prognosis compared to those with low scores. This implies score specificity across cancers. Wang et al. analyzed TME scores of MB data in the GSE85217 dataset, finding the SHH subgroup had the highest tumor immune and stromal scores (Li et al., 2020). However, Group3 exhibited the highest immune and stromal scores in our study. Furthermore, we analyzed the correlation between TME and TMErisk scores. TMErisk scores negatively correlated with stromal scores (r = -0.43, P < 0.001) and immune scores (r = -0.37, P = 0.004).

The ESTIMATE algorithm necessitates bulk RNA expression data to analyze the immune microenvironment, and sequencing technology costs remain challenging. Therefore, we developed a model comprising eight signature genes that could accurately predict prognosis in MB patients. The eight TMErisk model genes also regulate the TME in other tumors. OLFML3 encodes the human olfactory medin-like 3, a stromal cell protein containing an olfactory mediator domain that promotes angiogenesis. Studies show OLFML3 can increase NK T cell recruitment to play an anti-tumor role in colorectal cancer (Stalin et al., 2021). Using an OLFML3 antibody may improve PD-1 immunotherapy sensitivity (Stalin et al., 2021). In a pan-cancer analysis, the related protein OLFML2B was highly expressed across 14 cancers and associated with immune-related pathways like macrophages and tumor-associated fibroblasts (Hu et al., 2022). Additionally, OLFML2B may influence chemotherapeutic drug resistance and serves as a new prognostic marker and potential immunotherapy target (Hu et al., 2022). OLFML2B primarily mediates PI3K-Akt signaling and extracellular matrix–receptor interaction, regulating the tumor immune microenvironment (Hu et al., 2022).

The time-dependent AUC and C-statistics indicated the immune microenvironment-based TMErisk score could predict prognosis in MB patients. The TMErisk score may supplement current molecular subtypes. mRNAsi analysis of MB patients showed higher stemness in Group4, the low-risk group, potentially aiding risk stratification. Drug sensitivity analysis revealed varying sensitivities to some drugs like cytarabine, cediranib, navitoclax, and fludarabine across TMErisk groups. This confirms the TME's role in regulating sensitivity to antitumor drugs, as demonstrated in recent studies (Quail et al., 2017; Wu et al., 2017). Therefore, our immune microenvironment-based gene signature may provide chemotherapy, targeted therapy, and immunotherapy guidance.

New MB therapies focus on identifying distinct MB molecular characteristics and exploring new therapeutic targets (J. Romer et al., 2005; J. T. Romer et al., 2004). However, MB TME studies remain limited in scale (Margol et al., 2015; Murata et al., 2018; Pham et al., 2016; Vermeulen et al., 2018) likely due to the minimal immune infiltration and immunosuppression characterizing the MB TME (Bockmayr et al., 2018; Margol et al., 2015; Vermeulen et al., 2018) posing large-sample research challenges. Yauch et al. found SHH signaling between tumor cells and the immune microenvironment plays a key role in tumor progression, representing a new therapeutic target (Yauch et al., 2008). This evidence implies TME phenotype correlates to immunotherapy responsiveness and prognostic outcomes (S. Wang et al., 2019; Zeng et al., 2019). The TMErisk score negatively correlated with immune checkpoint and HLA family gene expression, indicating immunosuppression in high-risk MB patients. For instance, HLA-DMA (r = -0.19, P = 0.004), HLA-DMB (r = -0.18, P = 0.005), CD48 (r = -0.13, P = 0.040), PDCD1LG2 (r = -0.25, P < 0.001), and TIGIT (r = -0.25, P < 0.001) expression negatively correlated with TME-risk scores (Figure S8).

CIBERSORT and MCPcounter algorithm analysis also revealed patients in the low TMErisk group were more likely to exhibit low CD8 T cell levels (Figures 6A, B; Figures S12A, B). Our immune cell analysis across MB subgroups suggested tumor-associated macrophages (TAMs) M1 were more abundant in the SHH subgroup (Figure S11), consistent with previous studies (Margol et al., 2015; Maximov et al., 2019). Among TAMs, monocytes can activate microglia to produce CCL2, binding CCR2 and promoting leukocyte exudation into the nervous system (Eferl, 2013; Maximov et al., 2019). SHH subgroup tumor cells can promote MB intracranial and hematogenous metastasis through these pathways (Garzia et al., 2018). Additionally, this subtype produces TGF-β, causing CD8 T cell immunosuppression (Gate et al., 2014). However, the mechanisms by which immune cells affect MB patient prognosis require further exploration.

We utilized the OCLR machine learning algorithm to calculate mRNAsi in MB patients. We also calculated mRNAsi, gene signatures, and TMErisk scores across MB molecular subtypes. An inverse TMErisk score-mRNAsi association was observed, consistent with a previous study (Lian et al., 2019). These findings suggest the immune microenvironment may impact the stem cell transcriptome. We derived a prediction model based on eight signature genes whose expression significantly associated with stemness from the immune microenvironment.

The strength of our study is that we developed a TMErisk prognostic model based on a large dataset of MB patients in China and validated it in an external dataset.

Limitations include the model's moderate predictive power in external datasets requiring confirmation in populations. Due to strict tumor tissue library construction and sequencing criteria, only 94 cases underwent whole genome sequencing analysis. Non-qualifying samples were re-extracted and retested for inclusion in future studies. Finally, while our model shows some advantages for prognosis prediction and immunotherapy response prediction in MB patients, the limited immunotherapy-receiving MB patients prevented exploring correlations between TMErisk score and immunotherapy response, necessitating future genetic testing in such patients.

In conclusion, we constructed and validated an eight signature gene TME-based risk-scoring model that can predict OS in MB patients. The model associated with immune checkpoint gene expression, immune cell infiltration, tumor immunogenicity, and mRNAsi. The risk model we developed may supplement current WHO molecular subtypes and has promising application prospects in MB patients. Additionally, this study identified candidate molecular targets for further TME immunotherapy exploration in MB. However, the distinct TME landscapes across MB subtypes pose new research challenges.

Funding

This work was supported by grants from the National Natural Science Foundation of China (82172608, 81902975) and the National Key R&D Program of Science and Technology Innovation 2030 (2022ZD0210100).

Competing Interests

The authors have no relevant financial or non-financial interests to disclose

Authors' contributions

All authors contributed to the study conception and design. DH, ZJ, TJ and XJ were involved in research design. DH and ZJ were responsible for data collection and analysis. ZJ, HL, and DH prepared the manuscript. ZH, ZW, ZY, XC,YJ and DH prepared the chart. TJ, WZ, RL, YH, XQ, CL, JL,SL and HL made revisions to the manuscript. All authors confirm that they contributed to manuscript review and critical revision for important intellectual content, and that they read and approved the final draft for submission. All authors confirm that they contributed to manuscript review and critical revision for important intellectual content, and that they read and approved the final draft for submission. All authors read and approved the final manuscript

Availability of data and material

The original data are available from the corresponding author and GSE85217.

Ethics approval

This study was performed in line with the principles of the Declaration of Helsinki. This study was approved by the Ethics Committee of Beijing Tiantan Hospital, Capital Medical University (KY 2019-098-01).

Consent to participate

All participants and their families signed informed consent forms.

Consent for publication

Not applicable

Acknowledgements

We would like to thank the patients and their families for their trust and participation in this study. We appreciate the support of the China National Gene Bank.

Aran D, Hu Z, Butte AJ (2017) xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 18: 220. 10.1186/s13059-017-1349-1

Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, Fröhling S, Chan EM, Sos ML, Michel K, Mermel C, Silver SJ, Weir BA, Reiling JH, Sheng Q, Gupta PB, Wadlow RC, Le H, Hoersch S, Wittner BS, Ramaswamy S, Livingston DM, Sabatini DM, Meyerson M, Thomas RK, Lander ES, Mesirov JP, Root DE, Gilliland DG, Jacks T, Hahn WC (2009) Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462: 108-112. 10.1038/nature08460

Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, de Reyniès A (2016) Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 17: 218. 10.1186/s13059-016-1070-5

Bockmayr M, Mohme M, Klauschen F, Winkler B, Budczies J, Rutkowski S, Schuller U (2018) Subgroup-specific immune and stromal microenvironment in medulloblastoma. Oncoimmunology 7: e1462430. 10.1080/2162402X.2018.1462430

Cavalli FMG, Remke M, Rampasek L, Peacock J, Shih DJH, Luu B, Garzia L, Torchia J, Nor C, Morrissy AS, Agnihotri S, Thompson YY, Kuzan-Fischer CM, Farooq H, Isaev K, Daniels C, Cho BK, Kim SK, Wang KC, Lee JY, Grajkowska WA, Perek-Polnik M, Vasiljevic A, Faure-Conter C, Jouvet A, Giannini C, Nageswara Rao AA, Li KKW, Ng HK, Eberhart CG, Pollack IF, Hamilton RL, Gillespie GY, Olson JM, Leary S, Weiss WA, Lach B, Chambless LB, Thompson RC, Cooper MK, Vibhakar R, Hauser P, van Veelen MC, Kros JM, French PJ, Ra YS, Kumabe T, López-Aguilar E, Zitterbart K, Sterba J, Finocchiaro G, Massimino M, Van Meir EG, Osuka S, Shofuda T, Klekner A, Zollo M, Leonard JR, Rubin JB, Jabado N, Albrecht S, Mora J, Van Meter TE, Jung S, Moore AS, Hallahan AR, Chan JA, Tirapelli DPC, Carlotti CG, Fouladi M, Pimentel J, Faria CC, Saad AG, Massimi L, Liau LM, Wheeler H, Nakamura H, Elbabaa SK, Perezpeña-Diazconti M, Chico Ponce de León F, Robinson S, Zapotocky M, Lassaletta A, Huang A, Hawkins CE, Tabori U, Bouffet E, Bartels U, Dirks PB, Rutka JT, Bader GD, Reimand J, Goldenberg A, Ramaswamy V, Taylor MD (2017) Intertumoral Heterogeneity within Medulloblastoma Subgroups. Cancer Cell 31: 737-754.e736. 10.1016/j.ccell.2017.05.005

Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z (2017) Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 18: 248-262. 10.1016/j.celrep.2016.12.019

Chen F, Zhuang X, Lin L, Yu P, Wang Y, Shi Y, Hu G, Sun Y (2015) New horizons in tumor microenvironment biology: challenges and opportunities. BMC Med 13: 45. 10.1186/s12916-015-0278-7

Eferl R (2013) CCL2 at the crossroad of cancer metastasis. Jakstat 2: e23816. 10.4161/jkst.23816

Fridman WH, Pagès F, Sautès-Fridman C, Galon J (2012) The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer 12: 298-306. 10.1038/nrc3245

Garzia L, Kijima N, Morrissy AS, De Antonellis P, Guerreiro-Stucklin A, Holgado BL, Wu X, Wang X, Parsons M, Zayne K, Manno A, Kuzan-Fischer C, Nor C, Donovan LK, Liu J, Qin L, Garancher A, Liu KW, Mansouri S, Luu B, Thompson YY, Ramaswamy V, Peacock J, Farooq H, Skowron P, Shih DJH, Li A, Ensan S, Robbins CS, Cybulsky M, Mitra S, Ma Y, Moore R, Mungall A, Cho YJ, Weiss WA, Chan JA, Hawkins CE, Massimino M, Jabado N, Zapotocky M, Sumerauer D, Bouffet E, Dirks P, Tabori U, Sorensen PHB, Brastianos PK, Aldape K, Jones SJM, Marra MA, Woodgett JR, Wechsler-Reya RJ, Fults DW, Taylor MD (2018) A Hematogenous Route for Medulloblastoma Leptomeningeal Metastases. Cell 172: 1050-1062.e1014. 10.1016/j.cell.2018.01.038

Gate D, Danielpour M, Rodriguez J, Jr., Kim GB, Levy R, Bannykh S, Breunig JJ, Kaech SM, Flavell RA, Town T (2014) T-cell TGF-β signaling abrogation restricts medulloblastoma progression. Proc Natl Acad Sci U S A 111: E3458-3466. 10.1073/pnas.1412489111

Hakimi AA, Reznik E, Lee CH, Creighton CJ, Brannon AR, Luna A, Aksoy BA, Liu EM, Shen R, Lee W, Chen Y, Stirdivant SM, Russo P, Chen YB, Tickoo SK, Reuter VE, Cheng EH, Sander C, Hsieh JJ (2016) An Integrated Metabolic Atlas of Clear Cell Renal Cell Carcinoma. Cancer Cell 29: 104-116. 10.1016/j.ccell.2015.12.004

Hu P, Zhang X, Li Y, Xu L, Qiu H (2022) Pan-Cancer Analysis of OLFML2B Expression and Its Association With Prognosis and Immune Infiltration. Front Genet 13: 882794. 10.3389/fgene.2022.882794

Jia D, Li S, Li D, Xue H, Yang D, Liu Y (2018) Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY) 10: 592-605. 10.18632/aging.101415

Li C, Zou H, Xiong Z, Xiong Y, Miyagishima DF, Wanggou S, Li X (2020) Construction and Validation of a 13-Gene Signature for Prognosis Prediction in Medulloblastoma. Front Genet 11: 429. 10.3389/fgene.2020.00429

Lian H, Han YP, Zhang YC, Zhao Y, Yan S, Li QF, Wang BC, Wang JJ, Meng W, Yang J, Wang QH, Mao WW, Ma J (2019) Integrative analysis of gene expression and DNA methylation through one-class logistic regression machine learning identifies stemness features in medulloblastoma. Mol Oncol 13: 2227-2245. 10.1002/1878-0261.12557

Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, Ohgaki H, Wiestler OD, Kleihues P, Ellison DW (2016) The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol 131: 803-820. 10.1007/s00401-016-1545-1

Maeser D, Gruener RF, Huang RS (2021) oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 22. 10.1093/bib/bbab260

Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, Kaminska B, Huelsken J, Omberg L, Gevaert O, Colaprico A, Czerwinska P, Mazurek S, Mishra L, Heyn H, Krasnitz A, Godwin AK, Lazar AJ, Cancer Genome Atlas Research N, Stuart JM, Hoadley KA, Laird PW, Noushmehr H, Wiznerowicz M (2018) Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell 173: 338-354 e315. 10.1016/j.cell.2018.03.034

Margol AS, Robison NJ, Gnanachandran J, Hung LT, Kennedy RJ, Vali M, Dhall G, Finlay JL, Erdreich-Epstein A, Krieger MD, Drissi R, Fouladi M, Gilles FH, Judkins AR, Sposto R, Asgharzadeh S (2015) Tumor-associated macrophages in SHH subgroup of medulloblastomas. Clin Cancer Res 21: 1457-1465. 10.1158/1078-0432.Ccr-14-1144

Maximov V, Chen Z, Wei Y, Robinson MH, Herting CJ, Shanmugam NS, Rudneva VA, Goldsmith KC, MacDonald TJ, Northcott PA, Hambardzumyan D, Kenney AM (2019) Tumour-associated macrophages exhibit anti-tumoural properties in Sonic Hedgehog medulloblastoma. Nat Commun 10: 2410. 10.1038/s41467-019-10458-9

Murata D, Mineharu Y, Arakawa Y, Liu B, Tanji M, Yamaguchi M, Fujimoto KI, Fukui N, Terada Y, Yokogawa R, Yamaguchi M, Minamiguchi S, Miyamoto S (2018) High programmed cell death 1 ligand-1 expression: association with CD8+ T-cell infiltration and poor prognosis in human medulloblastoma. J Neurosurg 128: 710-716. 10.3171/2016.11.Jns16991

Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA (2015) Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 12: 453-457. 10.1038/nmeth.3337

Northcott PA, Buchhalter I, Morrissy AS, Hovestadt V, Weischenfeldt J, Ehrenberger T, Gröbner S, Segura-Wang M, Zichner T, Rudneva VA, Warnatz HJ, Sidiropoulos N, Phillips AH, Schumacher S, Kleinheinz K, Waszak SM, Erkek S, Jones DTW, Worst BC, Kool M, Zapatka M, Jäger N, Chavez L, Hutter B, Bieg M, Paramasivam N, Heinold M, Gu Z, Ishaque N, Jäger-Schmidt C, Imbusch CD, Jugold A, Hübschmann D, Risch T, Amstislavskiy V, Gonzalez FGR, Weber UD, Wolf S, Robinson GW, Zhou X, Wu G, Finkelstein D, Liu Y, Cavalli FMG, Luu B, Ramaswamy V, Wu X, Koster J, Ryzhova M, Cho YJ, Pomeroy SL, Herold-Mende C, Schuhmann M, Ebinger M, Liau LM, Mora J, McLendon RE, Jabado N, Kumabe T, Chuah E, Ma Y, Moore RA, Mungall AJ, Mungall KL, Thiessen N, Tse K, Wong T, Jones SJM, Witt O, Milde T, Von Deimling A, Capper D, Korshunov A, Yaspo ML, Kriwacki R, Gajjar A, Zhang J, Beroukhim R, Fraenkel E, Korbel JO, Brors B, Schlesner M, Eils R, Marra MA, Pfister SM, Taylor MD, Lichter P (2017) The whole-genome landscape of medulloblastoma subtypes. Nature 547: 311-317. 10.1038/nature22973

Pham CD, Flores C, Yang C, Pinheiro EM, Yearley JH, Sayour EJ, Pei Y, Moore C, McLendon RE, Huang J, Sampson JH, Wechsler-Reya R, Mitchell DA (2016) Differential Immune Microenvironments and Response to Immune Checkpoint Blockade among Molecular Subtypes of Murine Medulloblastoma. Clin Cancer Res 22: 582-595. 10.1158/1078-0432.Ccr-15-0713

Quail DF, Joyce JA (2017) The Microenvironmental Landscape of Brain Tumors. Cancer Cell 31: 326-341. 10.1016/j.ccell.2017.02.009

Romer J, Curran T (2005) Targeting medulloblastoma: small-molecule inhibitors of the Sonic Hedgehog pathway as potential cancer therapeutics. Cancer Res 65: 4975-4978. 10.1158/0008-5472.Can-05-0481

Romer JT, Kimura H, Magdaleno S, Sasai K, Fuller C, Baines H, Connelly M, Stewart CF, Gould S, Rubin LL, Curran T (2004) Suppression of the Shh pathway using a small molecule inhibitor eliminates medulloblastoma in Ptc1(+/-)p53(-/-) mice. Cancer Cell 6: 229-240. 10.1016/j.ccr.2004.08.019

Stalin J, Imhof BA, Coquoz O, Jeitziner R, Hammel P, McKee TA, Jemelin S, Poittevin M, Pocard M, Matthes T, Kaci R, Delorenzi M, Ruegg C, Miljkovic-Licina M (2021) Targeting OLFML3 in Colorectal Cancer Suppresses Tumor Growth and Angiogenesis, and Increases the Efficacy of Anti-PD1 Based Immunotherapy. Cancers (Basel) 13. 10.3390/cancers13184625

Tibshirani R (1997) The lasso method for variable selection in the Cox model. Stat Med 16: 385-395. 10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3

Vermeulen JF, Van Hecke W, Adriaansen EJM, Jansen MK, Bouma RG, Villacorta Hidalgo J, Fisch P, Broekhuizen R, Spliet WGM, Kool M, Bovenschen N (2018) Prognostic relevance of tumor-infiltrating lymphocytes and immune checkpoints in pediatric medulloblastoma. Oncoimmunology 7: e1398877. 10.1080/2162402x.2017.1398877

Wang H, Wu X, Chen Y (2019) Stromal-Immune Score-Based Gene Signature: A Prognosis Stratification Tool in Gastric Cancer. Front Oncol 9: 1212. 10.3389/fonc.2019.01212

Wang S, He Z, Wang X, Li H, Liu XS (2019) Antigen presentation and tumor immunogenicity in cancer immunotherapy response prediction. Elife 8. 10.7554/eLife.49020

Wu T, Dai Y (2017) Tumor microenvironment and therapeutic response. Cancer Lett 387: 61-68. 10.1016/j.canlet.2016.01.043

Xu WH, Xu Y, Wang J, Wan FN, Wang HK, Cao DL, Shi GH, Qu YY, Zhang HL, Ye DW (2019) Prognostic value and immune infiltration of novel signatures in clear cell renal cell carcinoma microenvironment. Aging (Albany NY) 11: 6999-7020. 10.18632/aging.102233

Yarchoan M, Hopkins A, Jaffee EM (2017) Tumor Mutational Burden and Response Rate to PD-1 Inhibition. N Engl J Med 377: 2500-2501. 10.1056/NEJMc1713444

Yauch RL, Gould SE, Scales SJ, Tang T, Tian H, Ahn CP, Marshall D, Fu L, Januario T, Kallop D, Nannini-Pepe M, Kotkow K, Marsters JC, Rubin LL, de Sauvage FJ (2008) A paracrine requirement for hedgehog signalling in cancer. Nature 455: 406-410. 10.1038/nature07275

Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, Mills GB, Verhaak RG (2013) Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 4: 2612. 10.1038/ncomms3612

Yu G, Wang LG, Han Y, He QY (2012) clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16: 284-287. 10.1089/omi.2011.0118

Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, Bin J, Liao Y, Rao J, Liao W (2019) Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res 7: 737-750. 10.1158/2326-6066.Cir-18-0436

The authors declare no competing interests.

TableS1.xls
Table S1 Dataset for the Validation Group;
TableS2.xls
Table S2 Cutoff value of Beijing Tiantan Hospital dataset and validation dataset;
TableS3.xls
Table S3 Clinical characteristics before and after propensity score matching;
TableS4.xls
Table S4 The analysis results of Gene Ontology (GO);
TableS5.xls
Table S5 The analysis results of Kyoto Encyclopedia of Genes and Genomes (KEGG).

Download PDF

Version 3

posted

You are reading this latest preprint version

An Immune Microenvironment-Associated Gene Signature Predicts Outcomes and Therapeutic Response in Pediatric Medulloblastoma

Status:

Version 3

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 3