Systematic transcriptomic analysis of childhood medulloblastoma identifies N6- methyladenosine-dependent lncRNA signatures associated with molecular subtype, immune cell infiltration, and prognosis

doi:10.21203/rs.3.rs-4810070/v1

Download PDF

Research Article

Systematic transcriptomic analysis of childhood medulloblastoma identifies N6- methyladenosine-dependent lncRNA signatures associated with molecular subtype, immune cell infiltration, and prognosis

https://doi.org/10.21203/rs.3.rs-4810070/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 28 Aug, 2024

Read the published version in Acta Neuropathologica Communications →

You are reading this latest preprint version

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 m6A-associated lncRNAs, a subgroup classification model was generated using the XGBoost machine learning algorithm, which had a classification accuracy > 90%, including for group 3 and 4 samples. All M6LSig genes were significantly correlated with at least one immune cell type abundance in the tumor microenvironment, and the risk score was positively correlated with CD4⁺ naïve T cell abundance and negatively correlated with follicular helper T cells and eosinophils. Knockdown of key m6A writer genes METTL3 and METTL14 in a group 3 medulloblastoma cell line (D425-Med) decreased cell proliferation and upregulated many M6LSig genes identified in our in silico analysis, suggesting that the signature genes are functional in medulloblastoma. This study highlights a crucial role for m6A-dependent lncRNAs in medulloblastoma prognosis and immune responses and provides the foundation for practical clinical tools that can be rapidly deployed in clinical settings.

Epitranscriptomics

immunity

long noncoding RNA

medulloblastoma

N6-methyladenosine

prognosis

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.

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:

$$\:\:M6LSig\:risk\_score={\sum\:}_{k=1}^{n}coef\:\left({x}_{k}\right)*{x}_{k}$$

where $\:{x}_{k}$ is the expression of the k^th gene, and $\:coef\:\left({x}_{k}\right)$ 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₂, 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.

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

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.

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.

MB Medulloblastoma

WNT WNT signaling-activated

SHH Sonic hedgehog signaling-activated

Grp3 Group 3

Grp4 Group 4

CNA Copy number alteractions

scRNA-seq Single cell RNA sequencing

lncRNA Long noncoding RNA

circRNA circular RNA

ncRNA noncoding RNA

m6A N6-methyladenosine

CRC colorectal cancer

TAM tumor-associated macrophage

GBM Glioblastoma

PCA Principal component analysis

WGCNA Weighted gene co-expression network analysis

M6LSig m6A-associated lncRNA signature

cph Cox proportional hazards

ROC Receiver operating characteristic

LDA Linear discriminate analysis

SVM Support vector machine

GB Gradient boosting

RF Random Forest

KNN K-nearest neighbors

DEG Differentially expressed genes

GSEA Gene set enrichment analysis

TME Tumor microenvironment

TPM Transcripts per million

HR Hazards ratio

AUC Area under curve

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

Data sharing is not applicable to this article as no new datasets were generated during the current study. The code used for analysis and visualization is made available at https://github.com/kandarpRJ/MB_epitranscriptomics/tree/main/m6A_lncRNA_project

Competing interests

The authors declare no competing interests.

Funding

Schamroth Project funded by Ian’s Friends Foundation to RJP and GJ. The Hough Foundation made a grant to RJP and GJ. Mahaffey foundation grant to SS. JHACH CBDI pilot research grant to KJ. This study was supported by P30 CA006973 (JHU SKCCC) to RJP & CGE, NCI 5P30CA030199 (SBP) RW-R, R01NS124668-01A1 to RJP, and CPRIT Scholar to MDT.

Author’s contributions

Conceptualization: K.J. & R.J.P.; Data analysis: K.J.; Project Administration: R.J.P.; Resources: R.J.P., S.S., G.J., C.G.E., M.D.T.; Supervision: R.J.P., E.A., A.R.; Experiments & Validation: M.Y., K.K., O.S.; Visualization: K.J., O.S., M.Y., K.K; Writing—Original Draft Preparation: K.J., R.J.P.; Writing—Review and Editing: R.J.P., K.J. All authors read and approved the final manuscript.

Acknowledgements

Not applicable

Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, Hawkins C, Ng HK, Pfister SM, Reifenberger G et al: The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro Oncol 2021, 23(8):1231–1251.
Funakoshi Y, Sugihara Y, Uneda A, Nakashima T, Suzuki H: Recent advances in the molecular understanding of medulloblastoma. Cancer Sci 2023, 114(3):741–749.
Sharma T, Schwalbe EC, Williamson D, Sill M, Hovestadt V, Mynarek M, Rutkowski S, Robinson GW, Gajjar A, Cavalli F et al: Second-generation molecular subgrouping of medulloblastoma: an international meta-analysis of Group 3 and Group 4 subtypes. Acta Neuropathol 2019, 138(2):309–326.
Hovestadt 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.
Williamson D, Schwalbe EC, Hicks D, Aldinger KA, Lindsey JC, Crosier S, Richardson S, Goddard J, Hill RM, Castle J et al: Medulloblastoma group 3 and 4 tumors comprise a clinically and biologically significant expression continuum reflecting human cerebellar development. Cell Rep 2022, 40(5):111162.
Ray 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).
Luo Z, Xia M, Shi W, Zhao C, Wang J, Xin D, Dong X, Xiong Y, Zhang F, Berry K et al: Human fetal cerebellar cell atlas informs medulloblastoma origin and oncogenesis. Nature 2022, 612(7941):787–794.
Hovestadt V, Smith KS, Bihannic L, Filbin MG, Shaw ML, Baumgartner A, DeWitt JC, Groves A, Mayr L, Weisman HR et al: Resolving medulloblastoma cellular architecture by single-cell genomics. Nature 2019, 572(7767):74–79.
Leary SES, Packer RJ, Li Y, Billups CA, Smith KS, Jaju A, Heier L, Burger P, Walsh K, Han Y et al: Efficacy of Carboplatin and Isotretinoin in Children With High-risk Medulloblastoma: A Randomized Clinical Trial From the Children's Oncology Group. JAMA Oncol 2021, 7(9):1313–1321.
Martinez de Estibariz I, Jakjimovska A, Illarregi U, Martin-Guerrero I, Gutierrez-Camino A, Lopez-Lopez E, Bilbao-Aldaiturriaga N: The Role of the Dysregulation of Long Non-Coding and Circular RNA Expression in Medulloblastoma: A Systematic Review. Cancers (Basel) 2023, 15(19).
Joshi P, Katsushima K, Zhou R, Meoded A, Stapleton S, Jallo G, Raabe E, Eberhart CG, Perera RJ: The therapeutic and diagnostic potential of regulatory noncoding RNAs in medulloblastoma. Neurooncol Adv 2019, 1(1):vdz023.
Katsushima K, Joshi K, Perera RJ: Diagnostic and therapeutic potential of circular RNA in brain tumors. Neurooncol Adv 2023, 5(1):vdad063.
Qin N, Paisana E, Picard D, Leprivier G, Langini M, Custodia C, Cascao R, Conrad C, Peitzsch M, Stefanski A et al: The long non-coding RNA OTX2-AS1 promotes tumor growth and predicts response to BCL-2 inhibition in medulloblastoma. J Neurooncol 2023, 165(2):329–342.
Kesherwani V, Shukla M, Coulter DW, Sharp JG, Joshi SS, Chaturvedi NK: Long non-coding RNA profiling of pediatric Medulloblastoma. BMC Med Genomics 2020, 13(1):87.
Katsushima K, Lee B, Kunhiraman H, Zhong C, Murad R, Yin J, Liu B, Garancher A, Gonzalez-Gomez I, Monforte HL et al: The long noncoding RNA lnc-HLX-2-7 is oncogenic in Group 3 medulloblastomas. Neuro Oncol 2021, 23(4):572–585.
Lee B, Katsushima K, Pokhrel R, Yuan M, Stapleton S, Jallo G, Wechsler-Reya RJ, Eberhart CG, Ray A, Perera RJ: The long non-coding RNA SPRIGHTLY and its binding partner PTBP1 regulate exon 5 skipping of SMYD3 transcripts in group 4 medulloblastomas. Neurooncol Adv 2022, 4(1):vdac120.
Katsushima K, Pokhrel R, Mahmud I, Yuan M, Murad R, Baral P, Zhou R, Chapagain P, Garrett T, Stapleton S et al: The oncogenic circular RNA circ_63706 is a potential therapeutic target in sonic hedgehog-subtype childhood medulloblastomas. Acta Neuropathol Commun 2023, 11(1):38.
Fang Z, Mei W, Qu C, Lu J, Shang L, Cao F, Li F: Role of m6A writers, erasers and readers in cancer. Exp Hematol Oncol 2022, 11(1):45.
Chen L, He Y, Zhu J, Zhao S, Qi S, Chen X, Zhang H, Ni Z, Zhou Y, Chen G et al: The roles and mechanism of m(6)A RNA methylation regulators in cancer immunity. Biomed Pharmacother 2023, 163:114839.
Yang Y, Hsu PJ, Chen YS, Yang YG: Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res 2018, 28(6):616–624.
Chen Y, Miao L, Lin H, Zhuo Z, He J: The role of m6A modification in pediatric cancer. Biochim Biophys Acta Rev Cancer 2022, 1877(2):188691.
Wu WJ, Xiao F, Xiong Y, Sun GF, Guo Y, Zhou X, Hu GW, Huang K, Guo H: N6-methyladenosine (m6A)-connected lncRNAs are linked to survival and immune infiltration in glioma patients. Biosci Rep 2023, 43(5).
Li R, Chen H, Li C, Qi Y, Zhao K, Wang J, You C, Huang H: The prognostic value and immune landscaps of m6A/m5C-related lncRNAs signature in the low grade glioma. BMC Bioinformatics 2023, 24(1):274.
Sharma A, Wang Y, Ge F, Chen P, Dakal TC, Carro MS, Schmidt-Wolf IGH, Maciaczyk J: Systematic integration of m6A regulators and autophagy-related genes in combination with long non-coding RNAs predicts survival in glioblastoma multiforme. Sci Rep 2023, 13(1):17232.
Chen Z, Zhang W, Yan Z, Zhang M: Comprehensive analyses indicated the association between m6A related long non-coding RNAs and various pathways in glioma. Cancer Med 2023, 12(1):760–788.
Cao X, Geng Q, Fan D, Wang Q, Wang X, Zhang M, Zhao L, Jiao Y, Deng T, Liu H et al: m(6)A methylation: a process reshaping the tumour immune microenvironment and regulating immune evasion. Mol Cancer 2023, 22(1):42.
Xie P, Yan H, Gao Y, Li X, Zhou DB, Liu ZQ: Construction of m6A-Related lncRNA Prognostic Signature Model and Immunomodulatory Effect in Glioblastoma Multiforme. Front Oncol 2022, 12:920926.
Ni W, Yao S, Zhou Y, Liu Y, Huang P, Zhou A, Liu J, Che L, Li J: Long noncoding RNA GAS5 inhibits progression of colorectal cancer by interacting with and triggering YAP phosphorylation and degradation and is negatively regulated by the m(6)A reader YTHDF3. Mol Cancer 2019, 18(1):143.
Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, Sun G, Lu Z, Huang Y, Yang CG et al: m(6)A RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem Cells. Cell Rep 2017, 18(11):2622–2634.
Dong F, Qin X, Wang B, Li Q, Hu J, Cheng X, Guo D, Cheng F, Fang C, Tan Y et al: ALKBH5 Facilitates Hypoxia-Induced Paraspeckle Assembly and IL8 Secretion to Generate an Immunosuppressive Tumor Microenvironment. Cancer Res 2021, 81(23):5876–5888.
Toker J, Iorgulescu JB, Ling AL, Villa GR, Gadet J, Parida L, Getz G, Wu CJ, Reardon DA, Chiocca EA et al: Clinical Importance of the lncRNA NEAT1 in Cancer Patients Treated with Immune Checkpoint Inhibitors. Clin Cancer Res 2023, 29(12):2226–2238.
Joshi P, Jallo G, Perera RJ: In silico analysis of long non-coding RNAs in medulloblastoma and its subgroups. Neurobiol Dis 2020, 141:104873.
Joshi K, Wang DO, Gururajan A: The m6A-methylome in major depression: A bioinformatic analysis of publicly available datasets. Psychiatry Research Communications 2022, 2(4):100089.
Shapiro JA, Gaonkar KS, Spielman SJ, Savonen CL, Bethell CJ, Jin R, Rathi KS, Zhu Y, Egolf LE, Farrow BK et al: OpenPBTA: The Open Pediatric Brain Tumor Atlas. Cell Genom 2023, 3(7):100340.
McLeod C, Gout AM, Zhou X, Thrasher A, Rahbarinia D, Brady SW, Macias M, Birch K, Finkelstein D, Sunny J et al: St. Jude Cloud: A Pediatric Cancer Genomic Data-Sharing Ecosystem. Cancer Discov 2021, 11(5):1082–1099.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL: StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol 2015, 33(3):290–295.
Love MI, Huber W, Anders S: Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014, 15(12):550.
Leek JT: svaseq: removing batch effects and other unwanted noise from sequencing data. Nucleic Acids Res 2014, 42(21):e161.
Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008, 9:559.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA: Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol 2018, 1711:243–259.
Khan M, Wong K, Jardel D, Broggio J, Stiller C, McCabe M: MEDULLOBLASTOMA INCIDENCE AND SURVIVAL - A POPULATION BASED STUDY. Neuro-Oncology 2018, 20(suppl_5):v350-v351.
Sun CX, Daniel P, Bradshaw G, Shi H, Loi M, Chew N, Parackal S, Tsui V, Liang Y, Koptyra M et al: Generation and multi-dimensional profiling of a childhood cancer cell line atlas defines new therapeutic opportunities. Cancer Cell 2023, 41(4):660–677 e667.
Khanna V, Achey RL, Ostrom QT, Block-Beach H, Kruchko C, Barnholtz-Sloan JS, de Blank PM: Incidence and survival trends for medulloblastomas in the United States from 2001 to 2013. J Neurooncol 2017, 135(3):433–441.
Hill RM, Richardson S, Schwalbe EC, Hicks D, Lindsey JC, Crosier S, Rafiee G, Grabovska Y, Wharton SB, Jacques TS et al: Time, pattern, and outcome of medulloblastoma relapse and their association with tumour biology at diagnosis and therapy: a multicentre cohort study. The Lancet Child & Adolescent Health 2020, 4(12):865–874.
Zhang ZW, Teng X, Zhao F, Ma C, Zhang J, Xiao LF, Wang Y, Chang M, Tian Y, Li C et al: METTL3 regulates m(6)A methylation of PTCH1 and GLI2 in Sonic hedgehog signaling to promote tumor progression in SHH-medulloblastoma. Cell Rep 2022, 41(4):111530.
Singh R, Kim YH, Lee SJ, Eom HS, Choi BK: 4-1BB immunotherapy: advances and hurdles. Exp Mol Med 2024, 56(1):32–39.
Stoll A, Bruns H, Fuchs M, Völkl S, Nimmerjahn F, Kunz M, Peipp M, Mackensen A, Mougiakakos D: CD137 (4-1BB) stimulation leads to metabolic and functional reprogramming of human monocytes/macrophages enhancing their tumoricidal activity. Leukemia 2021, 35(12):3482–3496.
Cachia D, Eskandari R, McDonald DG, Infinger LK, Vandergrift WA, Varma AK, Patel SJ, Zukas AM, Lindhorst SM, Das A: Low-dose radiation followed by on-target inhibition of Galectin-3 in combination with anti-4-1BB monoclonal antibody regulates immune responses in Group 3 and Group 4 medulloblastoma mouse model (P11-13.001). Neurology 2023, 100(17_supplement_2):3857.

Tables 1 to 3 are available in the Supplementary Files section.

No competing interests reported.

Download PDF

Journal Publication

published 28 Aug, 2024

Read the published version in Acta Neuropathologica Communications →

Editorial decision: Accepted
08 Aug, 2024
Editor assigned by journal
30 Jul, 2024
Submission checks completed at journal
30 Jul, 2024
First submitted to journal
26 Jul, 2024

You are reading this latest preprint version

Systematic transcriptomic analysis of childhood medulloblastoma identifies N6- methyladenosine-dependent lncRNA signatures associated with molecular subtype, immune cell infiltration, and prognosis

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Methods

Data collection and pre-processing

Construction of co-expression networks and m6A-associated lncRNA gene selection

Survival analysis and gene selection for calculating the M6LSig risk score

M6LSig risk score-based prognostic model for medulloblastoma

Machine learning model to classify medulloblastoma samples by lncRNA gene expression

Functional characterization of differentially expressed genes between low- and high-risk groups

Association between M6LSig and immune cell infiltration in the tumor microenvironment

Cell lines and cell culture

Quantitative real-time PCR

Cell proliferation assay

Statistical analysis

Results

Identification of an m6A regulator and lncRNA gene co-expression network

Construction of the m6A-lncRNA gene signature

Machine learning models for subgroup classification based on m6A-associated lncRNA gene expression profiles

Functional and immune characterization of the m6A-assocaited lncRNA signature and risk score

Discussion

Conclusions

Abbreviations

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1