Clinical and Molecular Features in Medulloblastomas Subtypes in Children in a Cohort in Taiwan

Simple Summary Medulloblastoma (MB) was classified into four subgroups: WNT, SHH, group 3, and group 4. In 2017, 12 subtypes within 4 subgroups and 8 subtypes within non-WNT/non-SHH subgroups according to the heterogenous features were announced. In this study, we aimed to identify the heterogeneity of molecular features for discovering subtype specific factors linked to diagnosis and prognosis. We retrieved 70 MBs to perform RNA sequencing and a DNA methylation array. Integrated with clinical annotations, we classified 12 subtypes of pediatric MBs. We found that M2 macrophages were enriched in SHH β, which correlated with good outcomes of SHH MBs. The high infiltration of M2 macrophages may be an indicator of a favorable prognosis and therapeutic target for SHH MBs. Furthermore, C11orf95-RELA fusion was observed to be associated with recurrence and a poor prognosis. These results will contribute to the establishment of a molecular diagnosis linked to prognostic factors of relevance for MBs. Abstract Medulloblastoma (MB) was classified into four molecular subgroups: WNT, SHH, group 3, and group 4. In 2017, 12 subtypes within 4 subgroups and 8 subtypes within non-WNT/non-SHH subgroups according to the differences of clinical features and biology were announced. In this study, we aimed to identify the heterogeneity of molecular features for discovering subtype specific factors linked to diagnosis and prognosis. We retrieved 70 MBs in children to perform RNA sequencing and a DNA methylation array in Taiwan. Integrated with clinical annotations, we achieved classification of 12 subtypes of pediatric MBs in our cohort series with reference to the other reported series. We analyzed the correlation of cell type enrichment in SHH MBs and found that M2 macrophages were enriched in SHH β, which related to good outcomes of SHH MBs. The high infiltration of M2 macrophages may be an indicator of a favorable prognosis and therapeutic target for SHH MBs. Furthermore, C11orf95-RELA fusion was observed to be associated with recurrence and a poor prognosis. These results will contribute to the establishment of a molecular diagnosis linked to prognostic indicators of relevance and help to promote molecular-based risk stratified treatment for MBs in children.


Introduction
Medulloblastoma (MB) is a common malignant brain tumor in children. Demographics, clinical information, and molecular data are significantly predictive factors for survival. According to the 2016 WHO classification, four molecular subgroups: WNT, SHH, group 3 (G3), and group 4 (G4) are included in MBs [1]. The subtypes within the molecular subgroups are defined as 12 subtypes [2]. G3 and G4 are merged as non-WNT/non-SHH MBs and comprised of eight subtypes by Northcott et al. [3]. The diversity of clinical features, demographics, and genetic and cytogenetic aberrations exists in MB subtypes. Two subtypes are included in the WNT subgroup: α and β, which exhibit favorable outcomes. WNT α mainly exists in children and presents with monosomy 6. Four subtypes are included in the SHH subgroup: α, β, γ, and δ, with different age distributions. SHH α presents in children and has the following features: TP53 mutations; focal amplifications in MYCN, GLI2, and YAP1; and broad loss in 9q, 10q, 17p. SHH β presents in infants and is associated with a high metastatic rate. SHH β presents the worst outcomes, which is associated with focal PTEN deletion. SHH γ presents in infants and is enriched histologically by MBEN, which indicates favorable outcomes. SHH δ mainly presents in adults and shows a favorable outcome as SHH γ.
Recently, two independent studies have announced various subtypes in non-WNT/non-SHH MBs. Cavalli and colleagues identified the G3 (α, β, γ) and G4 (α, β, γ) subtypes [2]. Northcott and colleagues identified eight subtypes (I to VIII) in non-WNT/non-SHH MBs, which were recruited in the 2021 WHO CNS5 classification [3,4]. Usually, subtype II to IV belong to G3 and subtype V to VIII belong to G4 [5]. Subtype I represents the least common subtype, whereas subtype VIII is the most common and only consists in G4 [3,6]. Generally, no chromosome aberrations are found in subtype I, while i17q are enriched in subtype VIII [6]. MYC amplification is enriched in subtype II and III and is associated with poor outcomes (5-year survival: 50% in subtype II, 43% in subtype III) [5]. Subtype VII is associated with a favorable 5-year survival (85%) [6].
Gene expression and DNA methylation profiles are the current standard for MB subgrouping and subtyping. Recently, the similarity network fusion (SNF) method for clustering was proposed [7]. By integrating gene expression and DNA methylation data, MB subgroups can divide into various subtypes [2]. In the previous study, we collected childhood MBs to identify a molecular-clinical correlation and defined an adjusted Heidelberg risk stratification scheme for treatment protocol guidelines in multiple centers in Taiwan [8]. Different MB subtypes need to be classified based on molecular and clinical heterogeneity for establishing molecular diagnostic and prognostic markers.
In this study, we retrieved 70 childhood MBs to perform RNA sequencing (RNA-Seq) and a DNA methylation array to perform subtype clustering in Taiwan. Integrated with clinical annotations, we achieved classification of 12 subtypes of pediatric MBs in our cohort series with reference to the other reported series. We characterized high infiltration of M2 macrophages in SHH β, which may be an indicator of a favorable prognosis and a therapeutic target for SHH MBs. Furthermore, C11orf95-RELA fusion was observed and associated with recurrence and poor outcomes. These results will contribute to the establishment of a molecular diagnosis linked to prognostic factors of relevance and further help to promote molecular-based, risk-stratified treatment for MBs in children.

Patient Cohort
There were 70 MB cases collected from Taipei Veterans General Hospital (Taipei VGH) and Taipei Medical University Hospital (TMUH), retrieved between 1989-2019,

Retrieve of Clinical Data
The retrieved clinical data included age, sex, metastasis status, histological variant, follow-up, and death. The centers of the tumor locations were defined as midline of the fourth ventricle (Midline/4thV) and cortex-centered, including cerebellar vermis (CV), cerebellar hemisphere (CH), and cerebellar pontine angle (CPA) location tumors. We defined the status of metastasis at diagnosis as M0-1 and M2-3 according to Chang's operative staging system [9].

Gene Expression Profiles by RNA-Seq
RNA-Seq was performed as described in the previous study [8]. Briefly, RNA-Seq was run in a Nextseq 500 sequencing instrument (Illumina) for paired-end reads. Gene expression tables were extracted by Kallisto [10] and the tximport [11] package in the R environment. The RNA-Seq data of 70 MB cases are available in the Gene Expression Omnibus (GSE143940 and GSE158413).

Applying RNA-Seq to Identify Mutations
There were 73 clinically relevant mutations selected for mutations that were detected in this cohort series. These selected mutations were linked to DNA damage response (DDR), MB genesis, a genetic predisposition for MB, the MAPK and PI3K/mTOR pathways, and pediatric cancer predisposition syndromes [3,8,[12][13][14][15][16]. RNA-Seq raw data were aligned using HISAT2 [17], followed by variant calling using the HaplotypeCaller tool in GATK. Variants were annotated using ANNOVAR [18] based on COSMIC database [19], and all variants in IGV with alignment level were visualized [20].

Immune Cell Deconvolution
Cell type deconvolutions were estimated as described in the previous study [21]. Briefly, the scores of 64 cell types in 5 major cell populations were computed with the gene expression data set normalized to TPM level of 489 cell population specific markers with xCell [22]. The scores of 34 immune cell types were compared between MB subtypes. The resulting scores are presented in arbitrary units.

DNA Methylation Array Profiling
The DNA methylation array was performed as described in the previous study [21]. Raw data files were read and preprocessed using the capabilities of Minfi [23] and the ChAMP [24] package in the R environment.

Applying DNA Methylation Profiles to Identify Copy Number Variations
The genetic status of chromosomes or selected genes was deciphered from the methylation array data. Selected copy number variations were identified from array data by using the conumee package in the R environment, as previously described [2,25,26]. The log2 ratio of chromosomes or genes more than 0.2 was defined as gain (amplification), and that of less than −0.2 was defined as loss (deletion).

Similarity Network Fusion (SNF) Analysis for WNT and SHH Subtype Clustering
The SNF method was performed in the cohort series as described in the previous study [21]. Briefly, subtype clustering was performed by the SNFtool package in the R environment based on the top 1% of the most differentially expressed common genes (n = 216) and probes (n = 3211) from a previous study [2]. The parameters of SNF were referred to the previous study [21].

Random Forest (RF) for Non-WNT/Non-SHH Subtype Clustering
The subtyping of non-WNT/non-SHH MBs was based on a web-based classifier of MB G3/4 subgroups (https://www.molecularneuropathology.org/mnp, accessed on 14 August 2021), as described in the previous study [21]. Briefly, Illumina Infinium Methyla-tionEPIC array raw signal IDAT-files were uploaded and normalized by a two-factor linear model on log2 transition to the web-based classifier with adjustment for frozen derivatives and patient gender. The most differential 50,000 CpG loci were implemented to calculate a RF score between 0 and 1 with multinomial logistic regression for non-WNT/non-SHH subtypes prediction [27].

Survival Analysis
Overall survival (OS) analysis was based on the date of first tumor surgery (diagnosis date), last follow-up, and death. OS analysis based on the scores of various cell types or the expression of genes was performed by the Kaplan-Meier method by using the surv_cutpoint function within the survminer package in the R environment. The differences of survivals were assessed using the log-rank test. The association between categorized variables was determined by the Kruskal-Wallis test. A p value < 0.05 was considered statistical significance.

Subsection of Molecular Subgroups Were Identified by Integrative Gene Expression and DNA Methylation Profiles
We retrieved 70 pediatric MBs to perform RNA-Seq and 66 of this cohort for the DNA methylation array. By clustering analysis, three established molecular subgroups were identified: WNT (n = 8, 11.4%), SHH (n = 24, 34.3%), and non-WNT/non-SHH (n = 38, 54.3%) (Figure 1a,b). There were 20 cases (28.6%) with metastasis at diagnosis in this cohort. SHH presented the highest recurrent rate (n = 11, 45.8%), which correlated with the worst prognosis (5-year overall survival (OS) after recurrence: 12.8%). The male/female ratio was 1 in all MBs; however, the ratio was 0. We combined gene expression and the DNA methylation profile to perform subtype clustering. In this cohort series, MBs were classified into WNT (α, β), SHH (α, β, γ), and non-WNT/non-SHH (II to VIII) ( Figure S1). The characteristics including gender, age, histological variants, tumor location, metastasis status, survival, cytogenetic, genetic aberrations, and immune cell enrichment of subtypes were identified (Table S1). We further compared demographics and clinical annotations of SHH and non-WNT/non-SHH subtypes in our and SickKids cohorts (Tables S2 and S3).

Discussion
In a genomic study of a large cohort, four distinct molecular subgroups of MB were identified and reported: WNT, SHH, G3, and G4 [30]. These subgroups presented specific demographics, histology, metastatic status, and prognosis [31]. By integration of gene expression and DNA methylation profiles, various subtypes with distinct demographic and clinical features were identified [2]. G3 and G4 subgroups exhibit similarities in molecular and biological profiling and are formally defined as non-WNT/non-SHH MB, which comprise eight subtypes [4]. In this cohort, 70 MBs in children and infants were retrieved and subjected to RNA-Seq and DNA methylation array analysis. SNF method and random forest scores were applied for subtype clustering to refine genetic and cytogenetic landscape within subtypes. The demographic, clinical annotations, molecular, and immune features of MB subtypes in our cohorts are summarized in Table 1.
Upon the clinical results, no metastasis presented in WNT, which was associated with very good outcomes (Figure 1g-i). Females mainly exhibited WNT in our cohort (Figures 1c and 2c), however, the male/female ratio was approximately 0.8 in other studies [2,6]. Classic histology mainly existed in WNT in our cohort and the other study (Figures 1e and 2e) [2]. We observed the CTNNB1 mutation in all WNT (Figure 2h), and monosomy 6 in all WNT α, as previously described (Figure 2g) [2,32,33]. DDX3X mutation existed in approximately half of WNT α, which was consistent with a previous study (Figure 2h) [32]. Interestingly, only one male WNT patient presented with recurrence.
SHH β typically occurs in children and is associated with better outcomes in our cohort (Figure 3d,h,i, and Table S1). However, SHH β occurred in infants and was associated with worse outcomes in the other studies [2,6]. A clinical trial study reported that iSHH-II (equivalent to SHH α and SHH γ) had improved survival with reduced intensity therapy compared to iSHH-I (equivalent to SHH-β) [34]. Another clinical trial study enrolled infants with DNMB/MBEN histology, which is the majority in SHH γ, but it was closed prematurely due to an excess of relapses [35]. The treatment strategy remains a key factor to affect the prognosis in SHH MBs.
TP53 mutation in the R273H point was associated with poor follow-up in one SHH α in this cohort according to the previous study [8]. The TP53 mutation was enriched in SHH α, which was associated with worse outcomes [2,6]. The R273H mutation in TP53 can develop highly metastatic tumors in mice models [36,37]. MYCN and GLI2 were frequently co-amplified in our cohort and other cohorts [38]. MYCN and GLI2 amplification are risk factors for SHH MBs [28,39]. Consequently, treatment with SMO inhibitors by targeting MYCN and GLI2 in the SHH pathway might have efficacy for SHH MBs, which exhibited these two co-amplifications [40][41][42][43]. No PTEN deletion was identified in SHH β, whereas the deletion existed in SHH α and γ in our cohort (Figure 3d,k). SHH β presented with focal PTEN deletion, which is associated with high metastatic rates and worse survival [2].   In this study, we found that M2 macrophages and their associated genes: CCL2, CD68, CD163, CD204, CD206, CD209, CSF1R, and Dectin-1, were enriched in SHH β (Figures 3n and S3-5a). The enrichment of M2 macrophages and their associated gene expression were correlated with favorable outcomes of SHH MB in our and the SickKids cohort ( Figure 3o, Figure S3b, Figure S4b and Figure S5b). SHH MBs have significant immune signatures of T cells, fibroblasts, and macrophages [44]. During tumorigenesis, hypoxia induces angiogenesis and recruits immune cells, such as macrophages, to initiate a pro-or anti-tumor response in the tumor microenvironment (TME) [45]. The increased M2 macrophages in SHH MBs might be due to the increased CCL2, a neuroinflammatory cytokine, which could recruit and promote M2 macrophage polarization [29,46]. The infiltration of tumor-associated macrophages (TAMs) and increased expression of their associated genes, CD163 and CSF1R, were significantly observed in SHH MBs [47]. An increase in M1 macrophages was reported to correlate with good outcomes of MB [48]. Furthermore, macrophage reduction in TME were correlated with poorer outcomes of SHH MBs, and TAMs might be involved in inhibiting tumor growth in SHH MBs [46]. However, another study demonstrated that recruitment of M1 macrophages was correlated with poor outcomes of SHH MBs [49]. The above contrary studies illustrate the inconclusive and incomplete roles of TAMs to promote or suppress tumor growth in TME.
Among non-WNT/non-SHH MBs, subtype IV occurred in younger patients with a median age of 3 years and was associated with favorable outcomes in our cohort and the other study (Figure 4c,g,h, and Table S1) [6]. Subtype VIII occurred in older children with a median age of 10 years and enriched i17q in our cohort and the other study (Figure 4c,i, and Table S1) [6]. MYC amplification was reported as a risk factor for G3, which was mainly enriched in subtype II, which was associated with poor outcomes in our cohort and other studies (Figure 4g,h,j) [28,50]. MYCN amplifications were predominantly found in subtype V, followed by subtype VI, which was associated with poor outcomes (Figures 4g,h,j and S6) [6,50]. However, MYCN amplification was not found in the only subtype V patient (Figure 4j). CDK6 amplification, which was predominantly found in G4 MBs highly enriched in subtype VI in our cohort [3,51] (Figure 4j).
In this study, we found NK and NKT cells enriched in subtype II (Figure 4m). NK cells can migrate to TME and exhibit cytolytic activity to kill tumor cells directly without specific immunization. NK cells were found to exist in MBs in the previous studies [44,[52][53][54]. It was reported that NK cells can suppress SHH MB tumor growth in a Daoy xenografted mouse model [53]. NKT cells can only recognize glycolipids or lipid antigens presented by CD1d, which is a monomorphic class I HLA molecule. The expression of CD1d was reported in GBM and SHH MBs and could be a potential target for NKT cell immunotherapy [55,56]. On the other hand, some studies reported that MB can suppress NK cell attacks with TGF-β, which is an immune suppressive strategy used by tumor cells [57][58][59][60]. Therefore, subtype II MBs may secrete TGF-β to fight against the cytotoxicity of NK or NKT cells in TME and promote tumor progression. The roles of NK or NKT cells in the tumorigenesis of non-WNT/non-SHH tumors need further study.

Conclusions
In conclusion, we highlighted genomic diversities in MB subtypes in a cohort series in Taiwan. We combined two platforms: gene expression and DNA methylation profiles, for MB subtype clustering. Genetic aberrations and prognosis within subtypes were identified. We found high enrichment of M2 macrophages, and their associated genes may be an indicator of a favorable prognosis in SHH MBs. TAMs might be a therapeutic target to improve the prognosis of SHH MBs. These results will contribute to the establishment of a nationwide molecular diagnosis linked to a prognostic indicator of relevance for MBs in children.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14215419/s1, Figure S1: The clustering and relationship between subgroups and subtypes of MBs in this cohort; Figure S2: Fusion gene distributions in MB subgroups; Figure S3: (a) The scores of the infiltrating M2 macrophage cells in SHH MB subtypes in the SickKids cohort. p value calculated by Kruskal-Wallis test. (b) The overall survival based on high or low M2 macrophage cell infiltration in SHH MBs in the SickKids cohort. p value calculated by log-rank test; Figure S4 Figure S7: (a) The scores of the infiltrating NK and NKT cells in non-WNT/non-SHH MB subtypes in the SickKids cohort. p value calculated by Kruskal-Wallis test. The overall survival based on high or low NK or NKT cell infiltration in non-WNT/non-SHH MBs in our (b) and SickKids cohort (b). p value calculated by log-rank test; Table S1: Demography and clinical data of molecular subtypes in our cohort series of 70 childhood medulloblastomas (MBs) in Taiwan; Table S2: Comparison of demographic and clinical annotations of SHH MB subtypes in our and SickKids cohorts; Table S3: Comparison of demographic and clinical annotations of non-WNT/non-SHH MB subtypes in our and SickKids cohorts.