A five‐CpG signature of microRNA methylation in non‐G‐CIMP glioblastoma

Summary Aims DNA methylation has been found to regulate microRNAs (miRNAs) expression, but the prognostic value of miRNA‐related DNA methylation aberration remained largely elusive in cancers including glioblastomas (GBMs). This study aimed to investigate the clinical and biological feature of miRNA methylation in GBMs of non‐glioma‐CpG island methylator phenotype (non‐G‐CIMP). Methods Prognostic miRNA methylation loci were analyzed, with TCGA and Rennes cohort as training sets, and independent datasets of GBMs and low‐grade gliomas (LGGs) were obtained as validation sets. Different statistical and bioinformatic analysis and experimental validations were performed to clinically and biologically characterize the signature. Results We identified and validated a risk score based on methylation status of five miRNA‐associated CpGs which could predict survival of GBM patients in a series of training and validation sets. This signature was independent of age and O‐6‐methylguanine‐DNA methyltransferase (MGMT) promoter methylation status. The risk subgroup was associated with angiogenesis and accordingly differential responses to bevacizumab‐contained therapy. MiRNA target analysis and in vitro experiments further confirmed the accuracy of this signature. Conclusion The five‐CpG signature of miRNA methylation was biologically relevant and was of potential prognostic and predictive value for GBMs. It might be of help for improving individualized treatment.


| INTRODUC TI ON
Glioblastomas (GBMs) are the most common and devastating subtypes of primary central nervous system tumors. 1 Unfortunately, despite the multimodal treatment of surgical resection, radiotherapy, and chemotherapy, the reported median survivals of GBM patients were only 16-19 months. [1][2][3] Cancer-specific DNA methylation changes play important roles in cancer development and progression. The best-known epigenetic abnormality in cancers is promoter-specific CpG island (CGI) hypermethylation of tumor suppressor genes which consequently cause transcriptional silencing. 4 Altered DNA methylation affected the expressions of not only protein-coding genes but also noncoding RNAs (ncRNAs). 5 Among those ncRNAs, microRNAs (miRNAs), the 20-22 nucleotides small ncRNAs, have been demonstrated to have multiple roles in the pathogenesis of cancers. 6 It has been reported that miRNAs could be regulated by DNA methylation and abnormal methylation in miRNAs was closely correlated with cancer progression. 6,7 However, the biological and clinical implications of miRNA methylation abnormality were largely unstudied in cancers including GBMs. Glioma-CpG island methylator phenotype (G-CIMP) represents a distinct subgroup of glioma which is featured by genome-wide hypermethylated CGIs and favorable prognosis. 8 The G-CIMP+ tumors have already been broadly studied, while the relevance features of non-G-CIMP GBMs remain largely unclear.
In this study, we analyzed miRNA methylation data of non-G-CIMP GBMs from The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and Rennes cohort 9 to reveal the relationship between miRNA methylation and GBM survival. Bioinformatic methods and in vitro experiments were used to validate our results.

| GBM datasets
Rennes cohort of 77 newly diagnosed non-G-CIMP GBMs with clinical and genome-wide DNA methylation microarray data by Infinium HumanMethylation450k BeadChip (Illumina Inc, San Diego, CA, USA) was obtained from the ArrayExpress under the accession number "E-MTAB-4969." 9 All patients received standard adjuvant treatment of radiotherapy (RT) and concurrent temozolomide (TMZ). Public DNA methylation datasets of non-G-CIMP GBM samples were also downloaded from The Cancer Genome Atlas (TCGA) data portal, 10 and Gene Expression Omnibus (GEO) under the accession number "GSE60274." 11 (Detailed clinical data of and relative CpG information are listed in the Supporting Information S1) We also obtained clinical and DNA methylation data of LGGs from TCGA 12 and GSE48462. 13 Among the heterogeneous datasets, only those with age over 18 years old and a molecular diagnosis of non-G-CIMP tumors were included in this study. For survival analysis, patients with a follow-up data >1 month were included, in order to reduce the bias caused by noncancer death. 10 In addition, nontumor brain tissues were obtained from apparently healthy individuals or chronic epilepsy patients with pathological evidence of other neurological or psychiatric diseases in each dataset. The G-CIMP status was determined by K-means (k = 3) clustering on the 1503 probes reported by Noushmehr et al 14 MGMT (O-6-methylguanine-DNA methyltransferase) promoter methylation status was determined by a logistic regression model using two CpGs, that is, cg12434587 and cg12981137. 15 Batch effects from different datasets and platforms were adjusted by a nonparametric empirical Bayes approach (ber package). 16 Methylation level of each integrated CpGs was summarized as M-value. 17

| Construction and validation of a miRNA methylation-based risk score model
CpG probes were filtered by removing those targeting the X and Y chromosomes, containing a single nucleotide polymorphism (SNP) within five base pairs of the targeted CpG. We then selected probes F I G U R E 1 Schematic diagram of the probe selection workflow for the study annotated with miRNAs (n = 2448) for this study. The discovery phase was performed within TCGA and Rennes cohort (training sets).
Univariate Cox regression analysis with permutation test was performed within each training set. Potential prognostic CpGs with consistent survival correlation (permutation P < 0.2) in each training set were subjected to multivariate Cox regression analysis within the combined training set (TCGA and Rennes collectively), and those with a P value < 0.05 were finally selected for risk score modeling. The risk score formula was constructed by integrating the M-values of all selected CpGs which were weighted by their multivariate Cox regression coefficients after adjusted by patient age and MGMT promoter methylation status. 18,19 Patients were then classified into high-risk or low-risk groups with the cutoff point as the median risk score from the combined training set. The validation phase was performed in GSE60274 and datasets of LGGs and in particular those with wide-type IDH.

| Gene Set Enrichment Analysis (GSEA)
GSEA was performed to evaluate the functional gene expression profiles between the risk subgroups on reported gene sets from Molecular Signature Database (MSigDB), with nominal P value ≤ 0.05 for significance. 20

| MiRNA target gene prediction and pathway analysis
The online databases TargetScan  were considered to be significant.

| Pyrosequencing
Pyrosequencing was performed by Pyromark Q96 ID platform and

| Cell proliferation assay
Cells with different treatments were implanted in 96-well plates at 5 × 10 3 per well. At indicated time points, CCK-8 kit (Yeasen, Shanghai, China) was assayed for cell viability measurement.

| Wound-healing assay
Cell motility was assessed by wound-healing assay as described previously. 26 A scratch wound was generated by a 200 μL pipette tip on the confluent cell monolayers in 6-well plates. The spread of the wound closure was observed after 48 hours of the scratch.

| Identification of prognostic miRNA methylation loci from the training sets
The included cohorts and the workflow of probe selection were schematically presented in Figure 1, and patient characteristics were summarized in Table 1. By employing a multistep selection criterion, we identified a five-CpG panel of miRNA methylation that showed consistent prognostic significance in both training sets (Table 2).
Among the panel, two CpGs (eg, cg05744073 and cg08244382) were hypermethylated and one CpG (eg, cg13767001) was hypomethylated in GBMs, while the other two were not differentially methylated in GBMs ( negative correlation with OS, while two CpGs (eg, cg24082174 and cg13767001) with positive correlation (Table 2).
Accordingly, the risk score model was constructed as follows:

| Validation of the five-GpG miRNA methylation signature for prognostication
To validate the prognostic performance of the 5-CpG miRNA methylation signature, we applied it to the independent validation set of GSE60274. With the prespecified cutoff, patients were classified into a low-risk group (n = 28) and a high-risk group (n = 31).
Consistent with the training sets, low-risk patients were associated with longer OS than high-risk ones (log-rank P = 0.013; Figure 2B).
We also observed a significant correlation between progressionfree survival (PFS) and the miRNA methylation-based risk groups in Rennes cohort ( Figure 2C).
In addition, we applied the GBM-derived signature to independent validation cohorts of IDH wild-type LGGs. The miRNA methylation signature failed to yield significant OS differences between the risk subgroups within both grade III and II gliomas using their median risk scores as cutoffs, respectively, which supported the signature as a GBM-specific prognostic model ( Figure 2D).

| Molecular and clinical correlation of the 5-CpG miRNA methylation signature
Correlation with current established molecular features showed that the 5-CpG signature appeared not to be significantly correlated with TCGA gene expression subtypes, DNA methylation clusters, and MGMT promoter methylation status ( Figure 3A). Also, the signature seemed not to be correlated with treatments, gender, and age ( Figure 3A). GSEA on gene expression data of TCGA samples showed that the high-risk tumors were enriched with prooncogenic gene sets such as ErbB signaling pathway (P < 0.0001), MAPK signaling pathway (P = 0.029), pro-angiogenic gene sets such as hypoxia (P = 0.035), and VEGF pathway (P = 0.029), which might biologically explain the inferior survival of those high-risk tumors ( Figure 3B).

| High-risk patients appeared to be beneficial for bevacizumab therapy
As reported by GSEA, high-risk tumors seemed to be featured by upregulation of various pro-angiogenic gene sets ( Figure 3B).
Accordingly, we tested the potential survival benefits conferred by the anti-angiogenic agent bevacizumab as combined therapy within each risk subgroup. In Rennes cohort with available second-line therapies, we found that the addition of bevacizumab did confer significant OS benefits in high-risk tumors, but was associated with similar OS to bevacizumab-free therapy ( Figure 3C).

| The 5-CpG signature was an independent prognostic factor in non-G-CIMP GBMs
Within Rennes cohort with RT/TMZ, univariate Cox regression analysis showed that patient age, MGMT promoter methylation The survival correlation of the five-CpG signature within current GBM classification. A, The five-CpG signature predicted overall survival (OS) in both MGMT promoter methylated and unmethylated patients treated with both radiotherapy (RT) and temozolomide (TMZ). B, It was also correlated with different OS in subgroups of ≤60 or >60 y. C, The correlation between five-CpG signature and different prognoses was significant in proneural and neural subtypes and marginally significant in the classical and mesenchymal subtypes status, and our miRNA methylation signature were significantly associated with OS (Table 3). Multivariate Cox regression analysis further demonstrated the prognostic independence of the abovementioned factors (Table 3). Multivariate Cox regression model within the combined cohorts of TCGA, GSE60274, and Rennes cohort not only confirmed the prognostic independence of our miRNA methylation signature but also supported its treatment independence (Table 3).

| The prognostic value of the miRNA methylation signature with respect to current GBM classification
We also tested the prognostic interrelationship of the 5-CpG signature with known prognostic factors within available patients from the combined training and validation sets. We found that it could consistently predict OS within the subtypes of unmethylated F I G U R E 5 Target prediction results of signature associated miRNAs  Figure 4A), and the subgroups of ≤60 or>60 years ( Figure 4B). Regarding the TCGA expression subtypes, the signature was significantly associated with OS in the proneural and neural subtypes, and also yielded marginally significant OS difference in the classical and mesenchymal subtypes ( Figure 4C).

| Target gene prediction of the 5-CpG signaturerelated miRNAs
Targetscan, miRanda, and miRDB databases were used to predict the target genes of miR-132, miR-127, miR-433, miR-1284, miR-1248, and miR-759. To ensure the specificity and sensitivity of our prediction, we kept the identical targets predicted in all three databases without setting additional criterion. Totally 1578 target genes left for further functional analysis ( Figure 5). Among them, 139 genes were candidate targets for two miRNAs, seven genes were common targets of three miRNAs, and one was target of four miRNAs.

| Biological characteristics of predicted target genes
Predicted target genes were further analyzed using PANTHER GO-slim tools on MF, BP, and CC ( Table 4). The most enriched MF terms were MAP kinase activity, ubiquitin-like protein transferase activity, DNA binding, and RNA polymerase II transcription factor activity ( Figure 6A). The most enriched BP terms were regulation of transcription by RNA polymerase II, transcription by RNA polymerase II, and cellular protein modification process ( Figure 6B). The most enriched CC terms were nuclear chromatin and extracellular space ( Figure 6C). PANTHER pathway analysis showed EGF receptor F I G U R E 6 Bioinformatic analysis of predicted target genes. A, PANTHER GO-Slim biological process. B, PANTHER GO-Slim molecular function. C, PANTHER GO-Slim cellular component. D, PANTHER pathway enrichment. E, KEGG pathway enrichment analysis, relative genes were shown as well signaling pathway, cadherin signaling pathway, FGF signaling pathway, CCKR signaling pathway, PDGF signaling pathway, and Wnt signaling pathway were the most enriched pathways ( Figure 6D).
Then, we utilized ClueGO to make a KEGG pathway enrichment analysis ( Figure 6E, Table 5). The most enrichment terms were adherens junction, cell cycle, TGF-beta signaling pathway, ErbB signaling pathway, axon guidance, renal cell carcinoma, oocyte meiosis, and cellular senescence.

| MiR-1284 suppressed glioma cell proliferation and migration
To further validate the functional relevance of this miRNA methylation signature, we selected miR-1284 for in vitro experiments.
Pyrosequencing of cg20382675 showed that the miR-1284-associated CpG was consistently associated with high DNA methylation status in glioma cells, that is, U87MG, U251, T98MG, and SHG44 ( Figure 7A). Accordingly, the expression levels of miR-1284 were comparable in those glioma cells ( Figure 7B). However, after treated with 5-Aza-dC, we found that the expressions of miR-1284 were significantly decreased in U251 and U87MG, indicating a positive impact of DNA methylation on miRNA expression ( Figure 7C). By transferring miR-1284 mimics into U251, we established a miR-1284-overexpressed U251 model, which was validated by qPCR ( Figure 7D). CCK-8 analysis showed that over-expression of miR-1284 reduced cell viability of U251 ( Figure 7E). Flow cytometry analysis showed that overexpression of miR-1284 was also associated with lower frequency of tumor cells in S and G2 phase ( Figure 7F), and higher frequency of apoptotic cells ( Figure 7G). Finally, wound-healing assay showed that migration was inhibited by over-expression of miR-1284 ( Figure 7H).

| DISCUSS IONS
This study investigated genome-wide DNA methylation microarray data of miRNA-associated CpGs to explore the clinical value of TA B L E 5 KEGG pathway enrichment analysis of predicted target genes  <0.001  21  ACP1, ACTN4, BAIAP2, CREBBP, CTNND1, EP300, INSR,  MAPK1, MAPK3, MET, NLK, SMAD4, SNAI1, SNAI2,  SORBS1, SRC, SSX2IP, TGFBR1,  The renal cell carcinoma pathway includes HIF-α pathway and strongly correlates with VEGF and PDGF production. This is consistent with the above GSEA results on gene expression data of TCGA and reminds the effects of these signature relative miRNAs.
To further explore the biological relevance of the 5-CpG signature, we selected one miRNA (miR-1284) for functional experiments. MiR-1284 has been reported to have tumor-inhibiting roles in lung, gastric and ovarian cancers, [30][31][32]  Among the other panel-associated miRNAs, some have been reported to be implicated in glioma biology. MiR-132 was upregulated in GBMs and was a potential indicator of poor prognosis. 35,36 MiR-127 and miR-433 are both derived from an overlapping gene locus and colocalized within a cancer-associated genomic region. 37 MiR-127 was reported to promote GBM cell migration and invasion by targeting tumor-suppresser gene SEPT7. 38 MiR-433 was reported to be commonly dysregulated in GBMs and suppressed glioma cell proliferation, migration, invasion, and enhanced sensitivity to TMZ therapy. 39,40 Regarding miR-759 and miR-1248, no biological or clinical evidences have been reported in cancers so far.
Our study has several limitations. First, as a retrospective study, the identification and validation of the signature were based on open source databases which had already been uploaded before. The follow-up information of these researches could not be considered in our study. Also, clinical information such as drug data and recurrent therapy of some cases was not detailed enough, which made it hard to make more subtle analysis. Second, bias caused by the differences among these selected trial platforms should be considered even with compensatory statistical measure. More proof should be collected before conducting further trials. Third, we only performed in vitro study on one miRNA, more in vitro and in vivo studies are needed, especially those on GBM angiogenesis.
In conclusion, by analyzing genome-wide DNA methylation microarray data of miRNAs-associated CpGs, we presented the initial report on the prognostic relevance of aberrant DNA methylation in miRNA regions in GBMs. The identification of the biologically and clinically relevant miRNA methylation signature may represent a promising approach for optimizing prognostication of GBMs, and be of potential value for improving individualized treatment and antiangiogenic therapy in particular.

ACK N OWLED G M ENTS
This work was supported by grants from the National Natural Science Foundation of China (No: 81471266; 81671302; and 81672909).

CO N FLI C T O F I NTE R E S T
No potential conflicts of interest were disclosed.