Systematic transcriptomic analysis of childhood medulloblastoma identifies N6- methyladenosine-dependent lncRNA signatures associated with molecular subtype, immune cell infiltration, and prognosis
Kandarp Joshi, Menglang Yuan, Keisuke Katsushima, Olivier Saulnier, Animesh Ray, Ernest Amankwah, Stacie Stapleton, George Jallo, Michael D. Taylor, Charles G. Eberhart, Ranjan J. Perera
TL;DR
Researchers found that specific RNA modifications and long non-coding RNAs are linked to survival and immune responses in childhood brain tumors, offering new diagnostic and treatment possibilities.
Contribution
A novel m6A-associated lncRNA signature (M6LSig) is identified for medulloblastoma prognosis and immune cell infiltration.
Findings
A five-gene m6A-associated lncRNA signature (M6LSig) is significantly linked to overall survival in medulloblastoma patients.
A machine learning model using m6A-associated lncRNAs classifies medulloblastoma subgroups with over 90% accuracy.
Knockdown of METTL3 and METTL14 reduced tumor cell proliferation and upregulated M6LSig genes in a group 3 cell line.
Abstract
Medulloblastoma, the most common malignant pediatric brain tumor, is classified into four main molecular subgroups, but group 3 and group 4 tumors are difficult to subclassify and have a poor prognosis. Rapid point-of-care diagnostic and prognostic assays are needed to improve medulloblastoma risk stratification and management. N6-methyladenosine (m6A) is a common RNA modification and long non-coding RNAs (lncRNAs) play a central role in tumor progression, but their impact on gene expression and associated clinical outcomes in medulloblastoma are unknown. Here we analyzed 469 medulloblastoma tumor transcriptomes to identify lncRNAs co-expressed with m6A regulators. Using LASSO-Cox analysis, we identified a five-gene m6A-associated lncRNA signature (M6LSig) significantly associated with overall survival, which was combined in a prognostic clinical nomogram. Using expression of the 67…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsRNA modifications and cancer · Cancer-related molecular mechanisms research · Circular RNAs in diseases
Introduction
Medulloblastoma (MB) is a tumor of the developing cerebellum usually seen in children. MB is broadly divided into four main molecular subgroups: WNT signaling-activated (WNT), sonic hedgehog signaling-activated (SHH) and non-WNT/non-SHH tumors split into two main subtypes: group 3 (Grp3), and group 4 (Grp4) [1, 2], although other subtypes have been identified [3]. Tumors in the WNT subgroup primarily harbor mutations in WNT signaling pathway genes such as CTNNB1, while SHH subgroup tumors harbor mutations in the SHH pathway such as PTCH1, SUFU, and SMO [4]. However, a lack of distinguishing mutations or copy number alterations (CNA) between Grp3 and Grp4 MB make them difficult to separate and manage. Recent attempts at subgroup classification have attempted to use epigenetic markers such as DNA methylation [4–6], while gene expression-based subgroup classification is challenging, as Grp3/Grp4 tumors exist as a continuum rather than as distinct subgroups [5]. While WNT MB is curable and certain non-TP53 mutant subtypes of SHH patients also have a good prognosis, Grp3 and Grp4 tumors tend to be highly aggressive and progress to relapse or metastasize [6]. Recent single-cell RNA sequencing (scRNA-seq) studies have revealed distinct cells of origin for Grp3 and Grp4 tumors [7, 8]. Furthermore, treatment of patients with high-risk Grp3 MB with chemoradiotherapy and carboplatin improves the 5-year event-free survival by 19% [9]. Therefore, Grp3 and Grp4 subclassification according to cell type and molecular mechanism is vital for risk stratification, patient care, and improving management and outcomes. However, there is currently no cost-effective nor rapid clinical classification method.
Long non-coding RNAs (lncRNA) and circular RNAs (circRNA) are regulatory non-coding RNAs (ncRNAs) that play important roles in MB gene regulation and tumor progression [10–12]. Several lncRNAs and circRNAs have been reported to have oncogenic roles and act as potential targets for therapy [13–17]. However, this field is still in its infancy, with many new candidates still being actively researched and discovered. Adding to the complexity of RNA regulation and function, chemical modifications of RNA bases such as with N6-methyladenosine (m6A), pseudo-uridine, and inosine are prevalent RNA modifications that add a layer of regulatory information on RNA. These modifications occur co-transcriptionally and contribute to many aspects of the RNA life cycle including stability, translation, and transport. m6A dynamics is tightly regulated by group of methyltransferase proteins (‘writers’) and demethylases (‘erasers’), with another group of proteins able to recognize modified RNA (‘readers’) to exert specific functions [18, 19]. These group of writers, erasers, and readers can act as m6A regulator genes that play important roles in m6A dynamics and, in turn, gene regulation [20]. Modifications of ncRNAs, including lncRNAs, are notable, as they directly affect RNA stability and regulation, but the role of m6A in pediatric cancers including MB has not been explored in detail [21]. Existing data suggest that m6A has prognostic value and affects immune cell infiltration in brain cancers [22–27], and a well-characterized lncRNA GAS5 has been shown to have an m6A-dependent role in tumor suppression in colorectal cancer (CRC). Methylated GAS5 transcripts are degraded by the m6A reader YTHDF3, which increases expression of YAP and in turn drives tumor progression [28]. Cui et al. showed that downregulation of m6A writers METTL3/METTL14 or m6A erasers ALKBH5/FTO reduces m6A on oncogenic ADAM19 mRNA, increasing its expression and driving self-renewal of glioblastoma stem cell, leading to aggressive tumor behavior and metastasis [29]. Dong et al. showed that ALKBH5 was responsible for m6A depletion on lncRNA NEAT1 in glioblastoma cells and increased transcript stability, allowing the transcriptional repressor SFPQ to relocate from the CXCL8 promotor to nuclear paraspeckle bodies, leading to upregulation of IL-8 in the tumor microenvironment and tumor-associated macrophage (TAM) recruitment. The immunosuppressive effect of TAMs allowed tumor progression and was partially rescued by IL-8 overexpression in ALKBH5-deficient GBM cells [30]. NEAT expression was associated with immune checkpoint blocker therapy in glioblastoma (GBM) patients and increased NEAT1 expression in M1 polarized macrophages led to secretion of inflammatory cytokines [31]. LncRNAs have previously been used to subclassify MB, but not with respect to m6A modifications [32]. We therefore reasoned that the diagnostic and prognostic accuracy of biomarkers for MB could be increased by integrating two different layers of regulatory information, lncRNAs and their m6A modification.
To achieve this, we analyzed bulk transcriptome data from over 400 MBs to identify m6A-associated lncRNA signatures. We constructed a co-expression network and identified modules showing representation of m6A regulators identified from a previous study [33]. Co-expressed lncRNAs in these modules associated with prognosis were selected as signature genes, which were used to calculate a risk score which, combined with other clinical parameters, predicted overall survival (OS) for risk stratification. Furthermore, using a machine learning model, the m6A-lncRNA signature was used for subgroup classification. Finally, we correlated the expression of each m6A-associated lncRNA and its signature-derived risk score with immune cell infiltration in the tumor microenvironment. Patient stratification into high and low risk groups based on risk score revealed significant differential expression of immune checkpoint and ligand genes. Collectively, we identify a new role for m6A-modified lncRNAs in immune-related functions and novel candidates for further functional characterization and testing as immune modulatory agents in MB.
Methods
Data collection and pre-processing
Bulk transcriptome data and associated metadata for three discovery cohorts of MB patients were analyzed. Raw RNA-seq counts were obtained from OpenPBTA (n = 119) [34] and Williamson et al. (n = 331) [5]. Alignment files in BAM format were retrieved from the St. Jude Cloud (n = 89) [35], with gene counts generated using stringtie [36]. Additionally, independent validation cohorts GSE164677 (n = 31) and the MAGIC cohort (n = 507) (unpublished) were used to validate the risk score predicting OS. Validation cohorts (GSE181293, GSE164677, and GSE85217; n = 839) were used for subgroup classification accuracy testing. A summary of the discovery and validation cohorts can be found in Supplementary Table 1. Samples with missing age and sex metadata were excluded from the analysis, and genes represented in all three cohorts were selected for analysis. Data were normalized using the vst function in DESeq2 [37], and subgroup clustering of the samples was visualized using a principal component analysis (PCA) plot. In order to perform cross-platform data integration and validation of the subgroup classification signature in RNA-seq and microarray-derived datasets, ComBat batch correction from the sva package [38] was performed on normalized gene expression profiles of discovery and validation cohorts.
Construction of co-expression networks and m6A-associated lncRNA gene selection
A co-expression network was constructed using the weighted gene co-expression network analysis (WGCNA) package in R [39]. An appropriate soft threshold was identified using power analysis for network construction. The constructed unsigned co-expression network contained a minimum module size of 50 genes and a maximum block size of 20,000. m6A-associated modules were identified based on representation of at least one out of 22 m6A regulatory genes within the module. Soft threshold selection to identify the power for network construction was performed with pickSoftThreshold, and an unsigned network of genes was constructed using blockwiseModules. LncRNA genes co-expressed with m6A regulator genes within m6A-associated modules were defined as m6A-associated lncRNA genes.
Survival analysis and gene selection for calculating the M6LSig risk score
LASSO-Cox regression-based feature selection was performed using the glmnet package to identify m6A-associated lncRNA genes associated with OS based on gene expression using the rms package in R, with age and sex as covariates. We selected genes with P-values < 0.05 as m6A-associated lncRNA signature (M6LSig) genes significantly associated with OS. A patient’s risk score was calculated based on M6LSig expression using the following formula:
where xk is the expression of the k^th^ gene, and coef (xk) is the coefficient of the Cox regression of the k^th^ gene.
M6LSig risk score-based prognostic model for medulloblastoma
Risk scores were calculated for all samples, and the median risk score was used to divide the samples into low- and high-risk groups based on lower and higher risk scores, respectively. Associations between risk score, age, and sex and OS were evaluated using a Cox proportional hazards (cph) model in the rms R package. Using the cph model generated using risk scores, age, and sex as variables and OS time and status, a points-based nomogram was developed to assess OS. Nomogram accuracy was assessed with calibration curves and receiver operating characteristic (ROC) curves in independent validation cohorts.
Machine learning model to classify medulloblastoma samples by lncRNA gene expression
Boruta package was used to identify the minimum set of genes capable of accurately subgrouping MBs. Briefly, the Boruta training function was run for 100 iterations to identify the gene set required for subgroup classification based on gene expression of lncRNAs co-expressed with m6A regulator genes in WGCNA modules.
Samples from the discovery cohort were divided into training and test sets in an 80:20 ratio. The training set was used to generate a model for subgroup classification. Eight algorithms were generated using the caret package in R for subgroup prediction based on M6LSig gene expression: glmnet, linear discriminate analysis (LDA), support vector machine (SVM), XGBoost, gradient boosting (GB), random forest (RF), K-nearest neighbors (KNN), and C50. We used the test set to evaluate model performance, precision, recall, and F1 scores, and multiclass ROC curves were calculated to assess the classification accuracy of each model. The performance of the top 5 models was evaluated using multi-class ROC curves using independent validation set not used for training the classifier models. The best performing model was chosen based on the highest multi-class accuracy.
Functional characterization of differentially expressed genes between low- and high-risk groups
Significant differentially-expressed genes (DEGs; adjusted P-value < 0.05) were identified between high- and low-risk groups using DESeq2 in R. We further ranked the expressing genes based on their test statistics for gene ontology and KEGG pathway gene set enrichment analysis (GSEA) using the clusterprofiler package. The significance of differential expression of 19 immune checkpoint-related genes and ligands was also examined.
Association between M6LSig and immune cell infiltration in the tumor microenvironment
Immune cell infiltration in the tumor microenvironment (TME) was identified using the cell type deconvolution method CIBERSORTx [40], with the LM22 dataset of 22 immune cell type gene expression patterns used as a reference set. Normalized gene expression, transcripts per million (TPM) for samples in the discovery cohort and LM22 reference datasets were used as references for immune cell deconvolution. Significant cell type proportions (P < 0.05) were selected, and correlations between individual M6LSig gene expression and risk scores with immune cell type proportions were calculated. Using R, the correlations between individual gene expression to cell type proportions and risk score were visualized. The complete analysis workflow is shown in Fig. 1.
Cell lines and cell culture
The human MB cell line D425-Med was purchased from Sigma-Aldrich (St. Louis, MO), which was STR profiled and Mycoplasma tested. D425 cells were cultured in DMEM/F12 with 10% FBS and 1% penicillin/streptomycin. Cells were grown in a humidified incubator at 37°C in 5% CO_2_, and the culture medium was replaced every 3 to 4 days. The cells were gently trypsinized (0.05%, Gibco, Thermo Fisher Scientific, Waltham, MA) for subculture.
Quantitative real-time PCR
Total RNA was purified as per the manufacturer’s instructions using the Direct-zol RNA Miniprep kit (Zymo Research, Irvine, CA). RNA yields were measured with a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific). High-Capacity cDNA Reverse Transcription Kits (Applied Biosystems, Waltham, MA) were used to reverse transcribe 500 ng RNA, with quantitative PCR performed using SYBR Green Master Mix assays (Applied Biosystems). The primer sequences are listed in Supplementary Table 2.
Cell proliferation assay
D425 cells were harvested in the logarithmic phase and cultured in 96-well plates for varying times for the cell proliferation assay. 20 μl MTS solution (Promega, Madison, WI) was added to each well, and the absorbance (optical density value) was measured at 490 nm on the EnVision 2105 microplate reader (PerkinElmer, Waltham, MA) after 2 hours.
Statistical analysis
Risk score was calculated according to Eq. 1. Samples were divided into high and low risk group based on the risk score value being higher or lower than the median risk score value. Differential expression analysis between high and low risk group samples was performed using DESeq2. Adjusted p-value cut-off 0.05 was used to identify significant differentially expressed genes between high risk and low risk patient groups. Cell proliferation and qRT-PCR data is represented as mean ± SD. Significance of differences between the groups were calculated using standard student’s t-test. P-value cut-off of 0.05 was used for determining statistical significance.
Results
Identification of an m6A regulator and lncRNA gene co-expression network
The discovery cohort consisted of 539 samples (OpenPBTA; Williamson et al. and St. Jude cohorts) and 18,872 genes. Samples with missing clinical metadata, such as age and sex, were excluded. The PCA plot shows sample clustering by subgroup for 469 samples from RNA-seq datasets after normalization (Fig. 2a). Gene expression profiles from 469 samples were used to build WGCNA co-expression modules. An unsigned co-expression gene network with 23 modules was constructed (Fig. 2b). Modules with at least one m6A regulator gene were identified as m6A-associated modules; 74 lncRNA genes (gencode v44) were significantly co-expressed with these m6A regulatory genes (Supplementary Table 3a-b).
Construction of the m6A-lncRNA gene signature
m6A-associated lncRNA genes were associated with patient OS. In total, 268 samples had survival metadata available, including age and sex covariates. Survival time was censored at 10 years. LASSO-Cox regression analysis was performed to identify m6A-associated lncRNAs associated with survival by selecting an appropriate lambda from the regression analysis (Supplementary Fig. 1a-b), and five lncRNAs were significantly associated with OS (Fig. 3a). An m6A-lncRNA signature risk score was calculated for each patient based on the expression of these five lncRNAs, and samples were divided into those above (high risk) or below (low risk) the median of the risk score. Using Cox multivariable regression, we found significant associations between risk score and OS, with low-risk scores showing better survival (hazard ratio (HR) = 0.42; 95%CI 0.27–0.66; P < 0.001). Similarly, patient age (HR 0.48, 95%CI 0.29–0.79; P = 0.004) was associated with OS and, while sex had an HR of 0.69 (95%CI 0.42–1.13; p = 0.146) for OS, this was not quite significant (Fig. 3b). However, given the marginal significance and previous reports of a significant group 3 bias for male vs female patients (2:1) and a significantly poorer five-year survival for males [41], we retained sex in subsequent prognosis calculations.
Discovery datasets were divided into training and test sets (90:10, respectively). Cox proportional hazards was used to model the training data and calculate the probability of 1-, 3-, and 5-year survival based on five m6A-lncRNA signature risk scores, age, and sex. The model was used to construct a nomogram (Fig. 4a), which was evaluated using calibration curves (Fig. 4b). Calibration curves provide information about agreement between observed versus expected probabilities of the estimates, where a linear relationship with a slope of the curve close to 1 indicates a high level of accuracy in predicting the OS of patients with MB. The best survival predictions were obtained using risk score, age, and sex for 1- and 3-year survival. Independent RNA-seq datasets (GSE164677 and MAGIC cohort) were used to validate the association of the five gene signature-derived risk score’s association with OS (Fig. 3c–d). Test data were used to generate ROC plots, and an area under curve (AUC) > 0.9 indicated a high accuracy of the five gene signature derived risk score in predicting OS (Fig. 4c–d).
Machine learning models for subgroup classification based on m6A-associated lncRNA gene expression profiles
Normalized and batch-corrected gene expression from the discovery cohorts (n = 469) were next used to train machine learning-based subgroup classification models. Out of 74 lncRNAs, 67 were identified as important features for subgroup classification using Boruta-based feature selection. Samples were randomly divided into training and test sets (80:20 ratio, respectively), and eight ML methods were assessed for accuracy in subgroup classification (Fig. 5a). LDA was the best-performing ML method, with a median accuracy of 97.3% (Table 1). Using test and independent validation cohorts (GSE164677, GSE181293, GSE164677, and GSE85217), the multiclass AUCs for the top five performing models were calculated, together with overall subgroup classification accuracy, per-group classification accuracy, and F-score. Table 2 shows the F-score for each model for individual subgroup classification accuracy and overall accuracy. Glmnet and XGBoost were the two best classifiers. XGBoost performed slightly better for Grp3 classification (Fig. 5b and Supplementary Fig. 2) and also performed the best for Grp3 vs Grp 4 classification (Fig. 5b).
Functional and immune characterization of the m6A-assocaited lncRNA signature and risk score
Differential expression analysis of genes between high- and low-risk patient groups identified 10,676 DEGs, represented as a volcano plot in Fig. 6a. Functional characterization of DEGs identified key biological processes involved in ribonucleoprotein complex biogenesis, axon development and axonogenesis, and regulation of nervous system development, amongst others (Fig. 6b). Immune cell proportions within the TME were next calculated by gene expression-based cell type deconvolution using CIBERSORTx and the LM22 dataset, which represents the expression of 547 genes in 22 immune cell types derived from scRNA-seq data. Cell type proportions per sample were calculated for 499 samples profiled by RNA-seq, which were then correlated with the expression of individual genes in the m6A-lncRNA signature and the risk score. All lncRNAs significantly correlated with at least one immune cell type in the TME, as determined by cell type deconvolution analysis (Fig. 7a). The differential expression of 19 immune checkpoint (IC) and ligand genes were also examined in high- and low-risk patients (Fig. 7b), and six IC genes were significantly differentially expressed (Fig. 7c).
METTL3 and METTL14 knockdown affects cell proliferation of Grp3 MB cells and expression of M6LSig genes
Finally, we knocked down METTL3 and METTL14, key m6A writers, using siRNAs in D425 Grp3 MB cells. There was an average 62.5% and 70% reduction in proliferation of D425 cells after 72 hours of METTL3 and METTL14 knockdown (Fig. 8a). qRT-PCR of the five genes in the signature showed that three lncRNAs (MATN1-AS1, LINC01963, and GAS5) were significantly upregulated in at least one of the knockdown condition samples (siMETTL3 or siMETTL14) compared with control (siNC) (Fig. 8b). RAET1E-AS1 and LINC02145 were not expressed in this cell line, as confirmed in RNA-seq data of D425 cells in the CCMA database [42] (expression 0 or < 1 CPM).
Discussion
Medulloblastoma (MB) is one of the most prevalent childhood brain tumors, with the highest incidence in children aged 10–14 years [43]. Children aged under three years have a higher risk of death from the disease, especially those with Grp3 or Grp4 tumors [6]. The recurrence rate for MB is high at 30% [44]. Different MB subgroups have different etiologies and cells of origin [4], and the resulting tumors progress differently, requiring personalized management. Subgroup classification requires histopathological analysis, but the Grp3 and Grp4 subgroups are challenging to differentiate, with costly and time-consuming DNA methylation-based classification currently the gold standard for Grp3 and Grp4 classification. There is therefore a clinical need for rapid point-of-care diagnostics for accurate subgroup classification and prognostication, especially in low-resource settings.
Our group previously reported a significant role for ncRNAs in MB subgroup classification and therapy [15–17, 32], and previous studies have reported the dynamics and dysregulation of m6A regulators in MB patients [21, 45]. RNA methylation is vital for regulating lncRNA-related disease progression and is an important therapeutic target. Nevertheless, m6A-based lncRNA regulation of tumor progression and clinical outcomes in MB have largely remained unexplored.
Here we systematically analyzed bulk transcriptomes from MB patients and integrated clinical traits and metadata to identify an m6A-associated lncRNA gene signature (M6LSig) with diagnostic and prognostic potential. LncRNAs co-expressed with m6A regulators are potentially regulated in an m6A-dependent manner. Our analysis revealed a five-gene signature significantly associated with prognosis in the context of age and sex. Due to intragroup heterogeneity and significantly different prognostic outcomes of subtypes within each subgroup, additional covariate of subgroup classification was not included in the model. The risk score calculated based on the expression values of these five genes was significantly associated with OS. High concordance of observed versus expected probabilities on calibration curves indicated that the model generated for predicting OS was highly accurate for one-, three-, and five-year survival. To aid clinical implementation, we generated a nomogram requiring the patient’s risk score, age, and sex for OS prediction. Eight ML-based classification models were also evaluated to determine the most precise algorithm for subgroup classification based on the 67-genes signature. XGBoost outperformed the other approaches, with > 90% accuracy for subgroup classification. These models have direct clinical application and are adaptable to diverse gene analysis platforms. Compared to other high throughput techniques dependent on quantifying bulk transcriptomes or DNA methylation profiles from patients, our small panel of genes are easier to scale and implement in clinical settings, even with limited resources availability. Although the roles of the lncRNAs forming this signature are not currently understood in MB, they have been implicated in many other cancers. For example, GAS5 has been reported to have prognostic value in CRC, where GAS5 regulates the transport and decay of YAP, an oncogene responsible for tumor progression in the disease. YTHDF3 can recognize and degrade m6A-containing GAS5 transcripts and hence negatively impact YAP degradation [28]. Further work is now needed to explore the functions and regulation of these lncRNAs.
Immune cell infiltration in the TME significantly impacts tumor progression, patient survival, therapy responses, and metastasis, including in MB, where CD8^+^ T cell infiltrates are associated with prognosis of MB patients. However, most MB tumors are considered ‘cold tumors’, i.e., the immune response is suppressed or inactive. Hence, immune modulation is a potentially important strategy for treating MB tumors, and novel targets for immune modulation hold important therapeutic and commercial potential. We explored correlations between expression of the five m6A-associated lncRNA genes and the abundance of 22 immune cell types within the MB TME, and all five genes showed significant positive or negative correlations with at least one of the 15 immune cell types (Fig. 6A). Correlating risk scores with immune cell type abundance in the TME identified a negative correlation with follicular helper T cells and eosinophils and a positive correlation with naive CD4^+^ T cells. Hence, the role of m6A-associated lncRNAs in regulating T cell and eosinophil infiltration and disease outcome is a novel therapeutic avenue for further exploration.
Clinically-applicable risk scores are important for personalized medicine. The result shown in Fig. 7b reports the dynamic expression patterns of the immune checkpoint and ligand genes in high- and low-risk patients, classified by their assigned subgroups. The plot shows the heterogeneity in the expression pattern of the immune checkpoint and ligand genes in the risk score-based stratified population of patients. Risk stratification by risk score is distinct from the subgroup classification-based risk stratification and that expression of immune checkpoint and ligand genes in samples with the same assigned subgroups are not similar. This indicates that an immune therapy approach and understanding the tumor microenvironment might involve additional risk factors than subgroup classification alone. Patients with a high-risk score showed significantly increased expression of 4–1BB, which is a positive stimulator of an anti-oncolytic immune response [46, 47]. Previous reports have shown positive outcomes for patients stimulated with anti-4–1BB antibodies, and combined low-dose radiation and anti-4–1BB antibodies have also been explored as a potential alternative therapy in MB patients [47, 48]. Taken together, our newly developed risk score has many clinical and personalized medicine applications in MB treatment and care.
The m6A writers METTL3 and METTL14 were significantly associated with proliferation of Grp3 cells in vitro. Zhang et al. similarly showed that METTL3 knockdown in SHH MB cell lines (DAOY and ONS-76) significantly impacted cell proliferation [45]. Hence, m6A regulatory genes are important for tumor progression and growth. Our results also indicate that manipulation of part of the m6A machinery has a significant impact on the regulation of many genes including three of the five lncRNAs forming part of the M6LSig in Grp3 MB cells, suggesting that these lncRNAs may be functional. Taken together, these results highlight a key role for m6A regulators in MB tumor growth and potentially influencing expression of prognostic lncRNAs in an m6A-dependent manner.
The results of this study highlight a crucial role for m6A-dependent lncRNAs in MB prognosis and immune responses. Further exploration and functional characterization of these identified genes and pathways is now required to determine the role of m6A in regulating immune cell infiltration. Our nomogram provides the basis for a useful and practical tool that can be rapidly deployed in clinical settings for translational applications.
Conclusions
We performed systematic in-silico analysis of publicly available bulk transcriptomes generated from the tumors of MB patients to identify a five gene m6A-associated lncRNA signature potentially linked to tumor progression and patient survival. Risk score derived from the five-gene signature has prognostic prediction ability. Another 67 gene lncRNA signature, which also contains the five prognostic lncRNAs, was identified with capability to perform subgroup classification with high accuracy. Our results indicate that m6A depletion by m6A writers (METTL3 and METTL14) knockdown affected the expression of a subset of these signature lncRNA genes in Grp3 MB cell line. Depletion of global m6A from transcriptome by m6A writers (METTL3 and METTL14) knockdown negatively affected the Grp3 MB cell line proliferation. Patients stratified by their risk score into high and low risk groups showed different gene expression patterns and biological processes. The expression of immune checkpoint genes and ligands was significantly different between high and low risk group patients. Expression of five lncRNA prognostic genes was also found to be significantly correlated with abundance of a subset of immune cells within tumor microenvironment. Taken together, these findings indicate the use of risk score-based stratification as potential marker for personalized immune therapy development application.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, Hawkins C, Ng HK, Pfister SM, Reifenberger G : The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro Oncol 2021, 23(8):1231–1251.34185076 10.1093/neuonc/noab 106PMC 8328013 · doi ↗ · pubmed ↗
- 2Funakoshi Y, Sugihara Y, Uneda A, Nakashima T, Suzuki H: Recent advances in the molecular understanding of medulloblastoma. Cancer Sci 2023, 114(3):741–749.10.1111/cas.15691 PMC 998607536520034 · doi ↗ · pubmed ↗
- 3Sharma T, Schwalbe EC, Williamson D, Sill M, Hovestadt V, Mynarek M, Rutkowski S, Robinson GW, Gajjar A, Cavalli F : Second-generation molecular subgrouping of medulloblastoma: an international meta-analysis of Group 3 and Group 4 subtypes. Acta Neuropathol 2019, 138(2):309–326.31076851 10.1007/s 00401-019-02020-0PMC 6660496 · doi ↗ · pubmed ↗
- 4Hovestadt V, Ayrault O, Swartling FJ, Robinson GW, Pfister SM, Northcott PA: Medulloblastomics revisited: biological and clinical insights from thousands of patients. Nat Rev Cancer 2020, 20(1):42–56.31819232 10.1038/s 41568-019-0223-8PMC 9113832 · doi ↗ · pubmed ↗
- 5Williamson D, Schwalbe EC, Hicks D, Aldinger KA, Lindsey JC, Crosier S, Richardson S, Goddard J, Hill RM, Castle J : Medulloblastoma group 3 and 4 tumors comprise a clinically and biologically significant expression continuum reflecting human cerebellar development. Cell Rep 2022, 40(5):111162.35926460 10.1016/j.celrep.2022.111162 PMC 9638015 · doi ↗ · pubmed ↗
- 6Ray S, Chaturvedi NK, Bhakat KK, Rizzino A, Mahapatra S: Subgroup-Specific Diagnostic, Prognostic, and Predictive Markers Influencing Pediatric Medulloblastoma Treatment. Diagnostics (Basel) 2021, 12(1).10.3390/diagnostics 12010061 PMC 877496735054230 · doi ↗ · pubmed ↗
- 7Luo Z, Xia M, Shi W, Zhao C, Wang J, Xin D, Dong X, Xiong Y, Zhang F, Berry K : Human fetal cerebellar cell atlas informs medulloblastoma origin and oncogenesis. Nature 2022, 612(7941):787–794.36450980 10.1038/s 41586-022-05487-2 · doi ↗ · pubmed ↗
- 8Hovestadt V, Smith KS, Bihannic L, Filbin MG, Shaw ML, Baumgartner A, De Witt JC, Groves A, Mayr L, Weisman HR : Resolving medulloblastoma cellular architecture by single-cell genomics. Nature 2019, 572(7767):74–79.31341285 10.1038/s 41586-019-1434-6PMC 6754173 · doi ↗ · pubmed ↗