This article is mentioned in:

Introduction

Glioblastoma (GBM) is the most common and aggressive primary brain tumor, accounting for ~30% of primary intracranial tumors worldwide (1). Because of its highly invasive growth and heterogeneous nature, the average survival time of patients with GBM is ~1 year (2). The standard GBM treatment of surgical resection followed by radiotherapy and postoperative chemotherapy has been improved dramatically, but the disease prognosis remains poor (3). While considerable progress has been made in the past decade in the understanding of the pathology of GBM (4), the underlying pathogenic mechanism of this tumor remains poorly understood.

Long non-coding RNAs (lncRNAs) have recently attracted considerable attention and become an important area of research. These RNAs are longer than 200 nucleotides and have no protein-coding potential (5). They can play complex and critical roles in tumor initiation and progression (6). For instance, the AGAP2-AS1 expression level was found to be elevated in GBM, and it functions as an oncogenic lncRNA to modulate GBM cell proliferation and apoptosis, suggesting that AGAP2-AS1 is a potential therapeutic target for GBM (7). Furthermore, the lncRNA CASP5 is upregulated in GBM tissues and promotes the malignant phenotypes of GBM (8). These discoveries suggested that lncRNAs could be excellent prognostic biomarkers and potential therapeutic targets for GBM. Although recent studies have identified some lncRNAs that exert regulatory activities during the development of GBM, these have mainly focused on the expression pattern of lncRNAs (9,10). At present, the regulatory mechanism of the vast majority of lncRNAs in GBM, especially DNA methylation, remains unclear.

DNA methylation, as an essential epigenetic modification, is involved in a variety of biological processes and mediates the dysregulation of gene expression (11). This modification can form a molecular basis for the silencing of tumor suppressors and the activation of oncogenes (12). In general, the hypermethylation of the gene promoter can downregulate or even silence gene expression, while the hypomethylation of the gene promoter tends to activate gene expression. For example, a previous study showed that the gene expression level of microRNA (miRNA/miR)-205 decreased in GBM tissues compared to controls, and lower expression was significantly associated with promoter hypermethylation (13). Tabu et al (14) found that promoter hypomethylation was an important determinant of CD133 overexpression in GBM, and this epigenetic event may be associated with the development of brain tumor-initiating cells expressing CD133. Recent studies have only described the aberrant methylation of some specific genes in GBM (15,16). However, the association between lncRNA methylation events and transcriptional changes at a global scale in GBM remains unknown.

With the improvement of high-throughput sequencing technology, large-scale Illumina Infinium Human Methylation 450 BeadChips (HM450k; Illumina, Inc.) and RNA-sequencing (RNA-seq) data have been applied for the analysis of cancer (17). The present study used a reannotation strategy to construct the DNA methylation profile of lncRNAs in GBM. The lncRNAs whose expression might be regulated by aberrant promoter methylation were determined according to two criteria, and the collective lncRNAs obtained in both scenarios were used for further analysis, including 314 hypermethylated lncRNAs (UhyperLncs) and 668 hypomethylated lncRNAs (UhypoLncs). Then, a methylation-mediated lncRNA regulatory network (MLRN) and functional analysis were used to elucidate the regulatory mechanism of lncRNAs and predict the functions of aberrantly methylated lncRNAs. Specially, it was found that four lncRNAs may have a good diagnostic and prognostic function. Finally, through the construction of drug-target association networks, the present study provided potential therapeutic targets and small-molecule drugs for GBM treatment. The present study enhanced the understanding of the regulatory mechanism of lncRNAs through DNA methylation and provided potential cancer biomarkers for the diagnosis and treatment of GBM.

Materials and methods

Data sources

The molecular data used in the present study were collected from various platforms. The DNA methylation data (level three) were generated on the Infinium HM450k platform (18). The HM450k data of tumor samples of GBM were downloaded from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov) and the data of normal samples, which contained 58 normal glial cell samples [GSE41826 (19)], were downloaded from the Gene Expression Omnibus (GEO) (20). The normal data were used in previous research and generated at the same organization via the same pipeline (21). The RNA-seq data of GBM from the Illumina RNAseqV2 platform (Illumina, Inc.) were downloaded from TCGA, which contained 156 tumor samples and five normal samples. To ensure the quality of the research, 50 tumor samples were specifically selected as the experimental dataset, among which DNA methylation data and RNA-seq data were both available. The remaining samples were used as a validation dataset (90 DNA methylation tumor samples and 106 RNA-seq tumor samples), where the normal samples in the experimental dataset were also applied to the validation dataset.

The human comprehensive gene annotation data were derived from GENCODE (release 19) (22). The experimental interactions between lncRNAs and miRNAs were collected from the starBase v2.0 database (23) and DIANA-LncBase v2 database (24). Human miRNAs and their targets were downloaded from starBase v2.0 and miRTarBase (release 7.0) (25). Both databases store manually curated collections of experimentally supported miRNA targets. The relationships between miRNAs and diseases were downloaded from the Human MicroRNA Disease Database (HMDD) (26). Comprehensive information about small molecule effects on miRNA expression was collected from SM2miR (27). Approved and experimentally validated drug target information were downloaded from DrugBank (v5.1.2) (28) and PharmGKB (29). Clinical data in XML format were downloaded from TCGA data portal for survival analysis.

Data normalization and construction of lncRNA methylation profile

The lncRNA and mRNA expression values were recalculated by TMM normalization and voom transformation for the expression data in raw read count format (30,31).

To evaluate the methylation values of a given probe, the methylation level of each probe was measured as a β-value, which is calculated as the ratio of methylated signal to the sum of the methylated and unmethylated signals. The range of β-values is from 0 (unmethylated) to 1 (completely methylated). To estimate the quality of the probe, the number of probes with missing values in all tumor samples was calculated. In total, 89,512 probes were removed, and the remaining missing values were filled in using the k-nearest neighbors method (32) with the knnImputation function in the DMwR package (https://mirrors.tuna.tsinghua.edu.cn/CRAN/src/contrib/DMwR_0.4.1.tar.gz). Finally, 392,867 CpG sites were obtained and used for constructing the lncRNA methylation profile.

Human genome annotation data were filtered to extract the lncRNA promoter information. Since the regulatory mechanism of lncRNA transcription is similar to the regulation of coding genes, 10 kb upstream from the transcriptional start site (TSS) was used for a relatively comprehensive range of lncRNA promoters (17). The 392,867 probes were mapped to the lncRNA promoter regions and only the probe closest to the TSS was used to determine the DNA methylation level of each lncRNA promoter (17,33).

Obtaining aberrant methylation-mediated lncRNAs

The limma package (34) was used to identify aberrantly methylated lncRNAs and differentially expressed lncRNAs between the tumor and normal samples, based on DNA methylation and RNA-seq data. The P-values were corrected using the Benjamini-Hochberg method (35), and only the lncRNAs with a corrected P≤0.05 were considered significant. To evaluate the correlation between methylation levels and expression levels of aberrantly methylated lncRNAs, the Pearson correlation coefficient (PCC) was calculated for each lncRNA between the methylation value and the corresponding expression value.

To comprehensively analyze lncRNAs whose expression levels were regulated by corresponding aberrant promoter methylation, the lncRNAs were selected for further study if they satisfied one of the following criteria: i) Hypermethylated lncRNAs and lncRNAs with significantly downregulated expression were overlapped as hypermethylated-underexpressed lncRNAs; similarly, hypomethylated lncRNAs and lncRNAs with significantly upregulated expression were overlapped as hypomethylated-highly expressed lncRNAs; and ii) the lncRNAs with PCC <-0.2 and P<0.05 were retained. Additionally, Student's t-tests and one-way ANOVA with Bonferroni's correction as the post hoc test were used to compare the validation dataset with the normal samples to obtain P-values; P<0.05 was considered to indicate a statistically significant difference.

Regulatory network construction and visualization

Experimentally supported interaction pairs of lncRNAs-miRNAs and miRNAs-targets were collected and integrated from multiple databases. A total of 32,452 non-redundant lncRNA-miRNA interactions and 729,567 miRNA-mRNA pairs were retained for further analysis. A candidate methylation-mediated lncRNA binary tuple (MlncBT) was defined as a lncRNA-mRNA interaction pair found to interact with the same miRNA in which the lncRNA was aberrantly methylated. In total, 3,531,816 potential MlncBTs were identified from lncRNA-miRNA and miRNA-mRNA interactions. To identify competing MlncBTs, the hypergeometric distribution was calculated to evaluate the significance of the shared miRNAs between each lncRNA and mRNA. The P-values were adjusted by false discovery rate (FDR), and MlncBTs with FDR≤0.05 were considered significant. In total, 916,628 candidate MlncBTs were retained for further identification.

Previous studies have indicated that increased lncRNA expression can enhance corresponding coding gene expression (36,37). The PCC of each MlncBT was calculated based on the lncRNA and mRNA expression profiles. PCC>0 and FDR≤0.05 were used as thresholds to screen out 28,721 MlncBTs comprising 262 lncRNAs and 9,441 mRNAs.

Finally, Cytoscape software (v3.7.0) (38) was used to visualize the regulatory network and divide the network into the hypermethylated lncRNA-mediated network and the hypomethylated lncRNA-mediated network.

Functional prediction of lncRNAs with different methylation patterns

Functional annotation of lncRNAs based on a co-expression network has been verified to be effective and accurate in previous studies (39,40). The mRNAs co-expressed with UhyperLncs and UhypoLncs in the regulatory network were used to perform functional enrichment analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to identify the significantly enriched biological processes and pathways, using the R package clusterProfiler (41). Only the GO terms and pathways with corrected P≤0.05 were retained to assess the potential functions of aberrantly methylated lncRNAs in GBM.

Identification of MlncBTs associated with GBM prognosis

To identify the clinical effect of MlncBTs, the patients were randomly divided into a training set and a test set, based on all the expression profile data (the sample sizes were the same in both groups). A Cox proportional hazards regression model was fitted to evaluate the association between the expression profile of each gene and patient survival in GBM. The prognostic index (PI) was adopted to classify the risk groups, as follows: PI=∑i=1nβiXi, where n was the number of survival correlated genes, βi was the Cox regression coefficient for genei and Xi was the expression level of genei in a corresponding patient. The median PI was used as a cut-off to divide patients in the training set into high- and low-risk groups. This PI model and cut-off point were also applied to the test set to divide the patients into high- and low-risk groups. Kaplan-Meier survival analysis and the log-rank test (P≤0.05) were performed to estimate the survival difference between the two patient groups. Moreover, receiver operating characteristic (ROC) analysis was performed to see whether the lncRNAs could be used as cancer biomarkers for early diagnosis of GBM.

Prediction of small molecule drugs for GBM treatment

To improve the accuracy of predicting small molecule drugs and targets, high-competing binary sub-networks were derived from the regulatory network by applying a PCC threshold >0.5. It was deemed that the sub-networks revealed a more stable regulatory relationship. Since the perturbation of miRNA expression could influence the expression level of many lncRNAs and mRNAs (42), the lncRNA-miRNA-mRNA triple sub-networks were further constructed, with the binary sub-networks as the background. A total of 231 GBM-related miRNAs were filtered from HMDD and 2,583 non-redundant miRNA-small molecule associations from SM2miR were integrated, with 1,084 associations of small molecules that could downregulate miRNA expression and 1,499 associations of small molecules that could upregulate miRNA expression. It was ensured that each miRNA in the sub-networks was among the GBM-related miRNAs that were screened from HMDD. The drug-target association network targeting miRNAs based on the hypermethylated lncRNAs was constructed by matching 1,084 associations with the hypermethylation-mediated lncRNA sub-network. Similarly, the drug-target association network targeting miRNAs based on the hypomethylated lncRNAs was constructed by matching 1,499 associations with the hypomethylation-mediated lncRNA sub-network. Furthermore, 14,256 gene-drug associations were integrated from DrugBank and PharmGKB. These gene-drug associations were used as seeds to match the drug-target networks targeting miRNAs to construct drug-target networks targeting mRNAs. Finally, drug-target association networks were constructed and illustrated using Cytoscape software (v3.7.0).

D-lnc (43) curates experimentally validated regulatory effects of drugs on lncRNA expression and contains 4,960 lncRNA-drug regulatory relationships for Homo sapiens. These lncRNA-drug regulatory relationships were integrated and matched with the associations in the drug-target networks that had been constructed in the present study, in order to verify the predicted drugs and targets. All analyses were performed using R 3.5.1 software (44).

Results

Characterization of global differences in DNA methylation and expression

To construct the DNA methylation profile of lncRNAs in GBM, a computational strategy was adopted to reannotate data from Infinium 450 k arrays into human lncRNA associated promoter regions. In the present study, 44,523 probes were located in 7,200 lncRNA promoter regions, in which each lncRNA had at least one probe mapped to the corresponding promoter region. Although each lncRNA had several probes mapping to the corresponding promoter region, only the probes closest to each TSS were retained to determine the DNA methylation status of the lncRNA promoters.

In some cancer types, previous studies have found that genome-wide DNA hypomethylation may be observed early in tumorigenesis, impacting genome stability and contributing to cellular transformation (12,45). After preprocessing the profiles, the present study focused on lncRNAs with significantly aberrant methylation between GBM and normal samples. In total, 5,567 aberrantly methylated lncRNAs were identified, including 1,214 hypermethylated lncRNAs and 4,353 hypomethylated lncRNAs. Meanwhile, hierarchical clustering analysis of these aberrantly methylated lncRNAs was performed according to the DNA methylation level (Fig. 1A). It was noted that these lncRNAs markedly differed between the tumor and normal tissues, and the number of hypomethyl-ated lncRNAs was much greater than that of hypermethylated lncRNAs. These results indicated that lncRNAs exhibit a more hypomethylated pattern during the occurrence and development of GBM. This global lncRNA hypomethylation may cause oncogene activation and genomic instability, and initiate tumorigenesis (46,47). Moreover, 2,567 differentially expressed lncRNAs were obtained in GBM compared to corresponding normal samples, in which 713 lncRNAs were upregulated, and 1,854 lncRNAs were downregulated. The expression levels of the differentially expressed lncRNAs are shown in a heatmap (Fig. 1B).

Figure 1 Unsupervised hierarchical clustering analysis of lncRNAs in GBM. Heatmaps of (A) aberrantly methylated and (B) differentially expressed lncRNAs in GBM relative to the normal control. On the x-axis, brown represents the GBM samples and dark blue represents the normal controls. The y-axis represents the lncRNAs. lncRNA, long non-coding RNA; GBM, glioblastoma.

Exploring aberrant methylation regulation patterns of lncRNAs in GBM

Since aberrant promoter methylation silences tumor suppressor genes or activates oncogenes (48), it was considered that aberrant promoter methylation of lncRNAs might be an important epigenetic regulator of lncRNA expression in GBM. To comprehensively analyze the regulatory effects of lncRNA aberrant promoter methylation on their expression in GBM, two criteria with unique biological and statistical significance were considered (Fig. 2A). For the first criterion, 213 hypermethylated-underexpressed lncRNAs were obtained by overlapping the hypermethylated lncRNAs and significantly underexpressed lncRNAs. Similarly, 288 hypomethylated-highly expressed lncRNAs were identified by overlapping the hypomethylated lncRNAs and significantly highly expressed lncRNAs. For the second criterion, 564 significantly negatively correlated lncRNAs were screened out by calculating the PCC of each aberrantly methylated lncRNA between methylation and expression level. The collective lncRNAs obtained in both scenarios was used for further analysis (Table SI), including 314 hypermethylated lncRNAs and 668 hypomethylated lncRNAs. Notably, 83 lncRNAs satisfied both criteria, among which 45 were hypermethylated and 38 were hypomethylated (Fig. S1). It was considered that the aberrant promoter methylation of these 83 lncRNAs was more likely to be an epigenetic regulator of their expression.

Figure 2 Regulatory patterns of lncRNA aberrant methylation in GBM. (A) The number of lncRNAs in different categories, where sncLnc represents lncRNAs with a significantly negative correlation between methylation level and expression level in GBM. (B) The pie chart shows the proportions of UhyperLncs and UhypoLncs relative to their locations with respect to protein-coding genes. The boxplots show comparisons of DNA methylation and expression of (C) UhyperLncs and (D) UhypoLncs. (E) Boxplot of the comparison of methylation levels between UhyperLncs, UhypoLncs and nsamLncs. (F) Boxplot of the comparison of expression levels between UhyperLncs, UhypoLncs and nsamLncs. lncRNA, long non-coding RNA; GBM, glioblastoma; UhyperLncs, hypermethylated lncRNAs; UhypoLncs, hypomethylated lncRNAs; nsamLncs, non-significantly aberrantly methylated lncRNAs.

A previous report suggested that different biotypes of lncRNAs perform distinct functions (49). Therefore, the two categories of lncRNAs were subdivided based on their location with respect to protein-coding genes to determine their types (Fig. 2B). The results showed that the majority of these lncRNAs were from antisense to protein-coding loci (antisense, 49.39%) and intergenic regions (lincRNAs, 40.53%).

To further dissect the promoter methylation patterns of lncRNAs associated with lncRNA expression, DNA methylation levels and expression levels between tumor and normal control samples in the validation dataset were used to evaluate statistical differences based on the aforementioned lncRNAs. As expected, it was found that UhyperLncs had higher methylation levels and lower expression levels in the tumor samples (Fig. 2C). On the contrary, the overall methylation levels of UhypoLncs in tumors tended to decrease, but the overall expression levels rose compared with normal control samples (Fig. 2D). These results indicated that aberrant promoter methylation of lncRNAs could affect their expression, and these lncRNAs may be cancer biomarkers for diagnosis and treatment in GBM. Furthermore, the DNA methylation levels and expression levels between UhyperLncs, UhypoLncs and non-significantly aberrantly methylated lncRNAs (nsamLncs, Fig. 2E and F) were observed; the overall methylation level of UhyperLncs was higher, but the overall expression level of UhyperLncs was lower compared with UhypoLncs. Notably, it was found that the nsamLncs, among the three groups, had the highest methylation level but the lowest expression level.

Construction of the MLRN and functional annotation

Previous studies have demonstrated that aberrant methylation of lncRNA promoters in tumors can lead to silencing or activation of lncRNA expression, and that dysregulation of lncRNA expression can regulate the expression of mRNAs by sharing common miRNA-binding sites with mRNAs (37). To further analyze how lncRNA dysregulation contributes to GBM progression through methylation-mediated epigenetic regulation, a comprehensive MLRN was constructed. The MLRN contained 262 lncRNAs, 9,441 mRNAs, and 28,721 MlncBTs. The MLRN was divided into a highly methylation-mediated lncRNA regulatory network (high MLRN) and a low methylation-mediated lncRNA regulatory network (low MLRN) based on the hypermethylated lncRNAs and hypomethylated lncRNAs that were singled out (Fig. 3A and B). There were 76 UhyperLncs, 4,110 mRNAs and 8,879 MlncBTs in the high MLRN and 186 UhypoLncs, 8,178 mRNAs, 19,842 MlncBTs in the low MLRN.

Figure 3 Construction of the methylation-mediated lncRNA regulatory network. (A) Regulatory network mediated by the hypermethylated lncRNAs; (B) regulatory network mediated by the hypomethylated lncRNAs. The node degree is indicated by the node size. lncRNA, long non-coding RNA.

Although tumorigenesis is a complicated dynamic process, a recent study demonstrated that the dysregulation of lncRNA methylation plays critical and complex roles during the development and progression of tumors (21). These aberrantly methylated lncRNAs mediate diverse biological functions, such as metabolism, cell apoptosis, angiogenesis, or many other cancer-related functions (50,51).

As the functional prediction of lncRNAs is hampered by the shortage of annotated information, functional annotation analysis of lncRNAs has frequently been conducted based on the guilt by association principle (52). To dissect whether different methylation patterns of lncRNAs in GBM correspond to distinct biological functions, the present study annotated lncRNAs with significantly enriched functional terms among the protein-coding genes that were co-expressed with the lncRNAs (53). A sum of 555 GO terms and 101 KEGG pathways were identified to be associated with the UhyperLncs, and 1,006 enriched GO terms and 66 pathways were obtained for the UhypoLncs. The results indicated that GO terms enriched by the UhyperLncs and the UhypoLncs involved cell proliferation, cell adhesion, cell apoptosis, cellular biological processes and methylation of related proteins, which are closely associated with the development of tumors (54) (Fig. 4A and B). Also, cell migration-related processes (such as 'nuclear migration' and 'positive regulation of fibroblast migration') were enriched in the UhyperLncs, while more biological transcription-related processes (such as 'RNA splicing' and 'regulation of DNA repair') were enriched in the UhypoLncs. For the KEGG pathway analysis, it was found that UhyperLncs and UhypoLncs were enriched in many cancer-related KEGG pathways (Fig. 4C and D). The UhyperLncs were mainly enriched in 'ErbB signaling pathway', 'FoxO signaling pathway', 'Wnt signaling pathway' and 'MAPK signaling pathway'. The UhypoLncs were mainly enriched in 'cell cycle', 'Alzheimer disease', 'colorectal cancer' and 'renal cell carcinoma'.

Figure 4 Aberrant methylation patterns of lncRNAs revealing different biological functions. Significantly enriched GO terms of (A) hypermethylated lncRNAs and (B) hypomethylated lncRNAs by co-expressed protein-coding genes. Significantly enriched KEGG pathways of (C) hypermethylated lncRNAs and (D) hypomethylated lncRNAs by co-expressed protein-coding genes. lncRNA, long non-coding RNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of four small nucleolar RNA host gene (SNHG) family lncRNAs associated with good prognosis

Previous research discovered that hub genes play essential roles in networks, and the top 10–20% of the nodes in networks are usually defined as hubs (36). Since hub nodes in the low MLRN had more degrees compared to those in the high MLRN, it was decided to analyze the hub nodes in the low MLRN. The present study found that SNHG16 was the hub lncRNA with the highest degree (degree, 1,012) in the low MLRN and another lncRNA, SNHG7, was also a hub node (degree, 272) in the low MLRN. Some lncRNAs in the SNHG family have been reported to be involved in glioma as potential oncogenes, and the dysregulation of their expression may promote the growth of glioma. For instance, the expression of SNHG7 was found to be upregulated in GBM tissues, and SNHG7 knockdown markedly suppressed cell proliferation and migration, while inducing cell apoptosis in GBM (54). Another lncRNA, SNHG18, was found to strengthen glioma cell radioresistance, and the expression was observably upregulated in clinical glioma tissues (55).

Notably, the present study found that SNHG16, SNHG7 and two other lncRNAs (SNHG9 and SNHG18) that were not in the MLRN all belonged to the UhypoLncs group, which meant that they were all hypomethylated in GBM and their increased expression in GBM might be affected by the corresponding promoter hypomethylation. It has been reported that SNHG7 expression is significantly increased in hypopharyngeal cancer, and that metformin decreases SNHG7 expression by mediating hypermethylation of the SNHG7 promoter (56), but the methylation status of the other three lncRNAs has not yet been reported, to the best of our knowledge. The present study also validated the methylation states of all probes annotated to these four lncRNA promoters. The results indicated that the methylation level of all CpGs on the SNHG7 promoter (Fig. 5A) and SNHG18 promoter (Fig. 5B) was significantly lower than that in paired normal samples. The methylation level of CpGs on the SNHG9 promoter and SNHG16 promoter was mostly decreased in GBM samples (Fig. S2A and B). Moreover, CpGs annotated to the four lncRNA promoters were more significant when they were closer to the TSS, which revealed that the CpGs closer to TSS were key in determining the methylation status of the lncRNA promoter. The diagnostic value of the four lncRNAs was further appraised to see whether they could be used as cancer biomarkers for early diagnosis of GBM. ROC analysis was performed according to the methylation and expression of the lncRNAs. The overall area under the ROC curve values for the diagnostic potential of the methylation and expression of these lncRNAs in GBM were >0.8 (Fig. 5C and D), which suggested that they could distinguish GBM samples from normal samples and might become diagnostic cancer biomarkers for GBM, especially SNHG16.

Figure 5 Four SNHG family lncRNAs are associated with good prognosis. Boxplots are presented with comparisons of methylation levels between GBM and normal samples of all CpG probes annotated to the (A) SNHG7 promoter and the (B) SNHG18 promoter. The probe sequences closer to the right in the figure were closer to the transcriptional start site. ***P≤0.001, ****P≤0.0001. Receiver operating characteristic analysis of (C) methylation and (D) expression of the four SNHG family lncRNAs was performed for GBM diagnosis. Survival analysis curves of the (E) SNHG16-HEATR1 MlncBT and (F) SNHG16-MUL1 MlncBT in the training set. Survival analysis curves of the (G) SNHG16-HEATR1 MlncBT and (H) SNHG16-MUL1 MlncBT in the test set. Numbers on the x-axis represent the number of living patients at that time. SNHG, small nucleolar RNA host gene; GBM, glioblastoma; AUC, area under the curve; MlncBT, methylation-mediated lncRNA binary tuple.

Furthermore, several studies have reported that SNHG16 is highly expressed in glioma tissues, and SNHG16 can upregulate the expression of some coding genes by interacting with miRNAs (57–59). In the present study, the lncRNA SNHG16 interacted with many mRNAs by combining with the target miRNAs in the low MLRN. To evaluate whether these MlncBTs were prognostic factors for GBM, the expression profiles were combined with clinical annotations and a subset of MlncBTs that significantly correlated with the overall survival of GBM was identified. It was found that two MlncBTs related to SNHG16 had a significant effect on survival [SNHG16-heat repeat containing 1 (HEATR1) and SNHG16-mitochondrial E3 ubiquitin protein ligase 1 (MUL1)]. HEATR1 is an mRNA related to GBM, and its expression in GBM tissues is significantly higher than that in normal samples (60), but the expression of the mRNA MUL1 in GBM has not been reported, to the best of our knowledge. Kaplan-Meier survival analysis of the training set indicated that the two MlncBTs could divide the GBM patients into two different risk groups with high and low PIs (Fig. 5E and F). The GBM patients in the test set could be separated into two groups based on the two MlncBTs (Fig. 5G and H). Taken together, it was concluded that the hypomethylation of the SNHG16 promoter might contribute to the activation of its expression, and then the high expression of SNHG16 may upregulate HEATR1 and MUL1 expression via the miRNA sponge mechanism. At the same time, these findings suggested their potential as novel prognostic signatures in GBM. All of the prognostic MlncBTs are summarized in Table SII.

Prediction of small molecule drugs for GBM treatment

As precision medicine becomes increasingly relevant in healthcare, the field of pharmacogenomics also continues to gain prominence in the clinical setting (29). It was inferred that some potential small molecules and targets that could be used for the treatment of patients with GBM may be identified by constructing drug-target association networks targeting miRNAs (Fig. 6A and B). In the hypermethylated sub-network, these potential small molecules may indirectly promote lncRNA expression by inhibiting the corresponding expression of miRNAs to achieve the purpose of treating GBM (Fig. 6A); for example, lenalidomide (downregulates hsa-miR-449a), emodin (downregulates hsa-miR-137) and marine fungal metabolite 1386A (down-regulates hsa-miR-320a). Previous studies have shown that lenalidomide and emodin may be used for the treatment of GBM (61,62). Besides, experiments have demonstrated that marine fungal metabolite 1386A has strong cytotoxicity towards cancer cells and can alter the miRNA profiles of MCF-7 breast cancer cells (63). It was deduced that 1836A might be a novel potential small-molecule treatment for GBM. In the hypomethylated sub-network, these potential small molecules might indirectly inhibit lncRNA expression by promoting the corresponding expression of miRNAs to treat GBM (Fig. 6B). Temozolomide (TMZ) is widely used to treat GBM, and a previous study demonstrated that some miRNAs, including hsa-miR-137, were significantly upregulated in GBM after treatment with both TMZ and Olea europaea leaf extract (64). Studies have revealed that the lncRNAs SNHG7, SNHG16 and ZFAS1 are upregulated in GBM or glioma tissues compared with non-tumor brain tissues (54,58,65), and cisplatin has shown good performance in patients with recurrent GBM (66). It was inferred that cisplatin might upregulate the expression of these lncRNAs by downregulating the expression of hsa-miR-449a. Notably, the present results indicated that TMZ may upregulate hsa-miR-137 and further down-regulate SNHG16. Furthermore, the present study sought to identify a combined treatment with different targets for use in patients with GBM by constructing drug-target networks targeting mRNAs (Fig. 6C and D). A total of 70 drugs (39 mRNA targets) were acquired in the hypermethylated sub-network (Fig. 6C) and 209 drugs (112 mRNA targets) in the hypomethylated sub-network (Fig. 6D). The partial validation results of the D-lncRNA database are discussed in the Discussion.

Figure 6 Construction of the drug-target association network. Drug-target networks targeting miRNAs based on (A) the hypermethylated lncRNAs and (B) the hypomethylated lncRNAs were constructed. Drug-target networks targeting mRNAs based on (C) the hypermethylated lncRNAs and (D) the hypomethylated lncRNAs were constructed. Blue circles, yellow circles, red circles and green diamonds represent lncRNAs, miRNAs, mRNAs and small-molecule drugs, respectively. lncRNA, long non-coding RNA; miRNA/miR, microRNA.

Discussion

Increasing evidence suggests that lncRNAs play crucial roles in carcinogenesis (6). In recent years, with the development of high-throughput sequencing technology, epigenetic regulation has become a hotspot in biomedical research. DNA methylation is an important pattern of epigenetic regulation, which can be involved in different biological processes by regulating the transcriptional activity of genes, including the occurrence and development of tumors (11).

The present study comprehensively investigated the changes in DNA methylation of lncRNA promoters in GBM. DNA methylation and expression profiles from TCGA and GEO were integrated to study the methylation regulation pattern of lncRNAs. The present study found that many lncRNAs in GBM exhibited epigenetic dysregulation in promoter regions, most of which were hypomethylated. There is evidence that reduced expression of lncRNAs is due to hypermethylation inhibition and increased expression of lncRNAs is due to hypomethylation activation (51). Sang et al (67) preliminarily explored the functions and pathways of coding genes in hepatocellular carcinoma by overlapping differentially methylated mRNAs and differentially expressed mRNAs. Therefore, to establish a strategy for predicting lncRNAs whose expression level might be associated with aberrant promoter methylation, the present study not only obtained hypermethylated-underexpressed lncRNAs and hypomethylated-highly expressed lncRNAs, but also retained lncRNAs with a significant negative correlation between methylation and expression level in GBM samples. In total, 314 hypermethylated lncRNAs and 668 hypomethylated lncRNAs were used for further analysis. A recent study predicted six lncRNAs that might improve the prognosis of GBM by using weighted gene co-expression network analysis, Cox regression and L1-LASSO penalization (68). Notably, one of the six lncRNAs (PRRT3-AS1) was among the hypermethylated lncRNAs identified in the present study. The results suggested that the aberrantly methylated lncRNAs screened in the present study may be biomarkers of GBM prognosis. Moreover, it was worth noting that 83 lncRNAs satisfied both criteria, and it was thought that these lncRNAs would have greater value in the future to experimentally verify their epigenetic regulation.

Functional analysis of lncRNAs has revealed that they can regulate gene expression at the transcriptional, post-transcriptional and epigenetic levels (1), and a competing endogenous relationship is one way in which lncRNAs influence the expression of mRNAs, by binding to the target miRNAs of mRNAs (36). The construction of the MLRN provided a global perspective to study the interactions between mRNAs and differentially methylated lncRNAs. Although it was confirmed that epigenetic modification of lncRNAs may be a modulator of their expression, it was unclear what important roles these lncRNAs may play in GBM pathogenesis. In the present study, both hypermethylated and hypomethylated lncRNAs were involved in many biological processes and pathways related to tumorigenesis and progression in GBM. Notably, it was found that aberrant methylation of either UhyperLncs or UhypoLncs may perturb specific common pathways, such as 'ErbB signaling pathway', 'FoxO signaling pathway', 'mTOR signaling pathway', 'glioma' and 'apoptosis'. As the most enriched pathway, ErbB [epidermal growth factor receptor (EGFR)] amplification and mutations are the most common oncogenic events in GBM, and oncogenic EGFR induces DNA methylation-mediated transcriptional silencing of tumor suppressors (69).

Furthermore, it was found that four SNHG family lncRNAs had an excellent diagnostic effect and clinical prognostic value. Three of the four lncRNAs (SNHG7, SNHG16, SNHG18) have been reported as oncogenic lncRNAs in glioma or GBM, and were highly expressed in tumors (54,55,58). In the present study, the expression of SNHG9 was upregu-lated in GBM, thus it was assumed that SNHG9 may also be a potentially carcinogenic lncRNA. Moreover, it was found that the four lncRNAs were hypomethylated, and one lncRNA (SNHG7) has been reported to improve the survival of patients with hypopharyngeal cancer by inhibiting the hypomethylation of its promoter (56). Furthermore, two MlncBTs (SNHG16-HEATR1 and SNHG16-MUL1) were associated with a good prognosis. Finally, the drug-target association networks were constructed to provide potential small molecule drugs and targets for the precise treatment of GBM. Additionally, it was found that doxorubicin can inhibit GAS5 expression by upregulating hsa-miR-449a in the present results. In the D-lnc database, it was also found that doxorubicin can downregulate GAS5 expression (data not shown).

There is no denying that the present research has some shortcomings which should be addressed. Since GBM is different from general tumors, sample acquisition is a problem, and aberrantly methylated lncRNA research should be replicated in a larger cohort. Moreover, due to technical and time constraints, it was not possible to verify the association between these lncRNA expression changes and promoter methylation events. Meanwhile, the lack of analysis of the high MLRN was a limitation to the present study. Further studies in animal models of GBM and brain tissues from patients with GBM may help to solve these issues.

In conclusion, the present study investigated the regulatory mechanism and functions of DNA methylation of lncRNAs in GBM by integrating multi-omics data. The identified GBM-related or clinically relevant lncRNAs (MlncBTs) might help to improve understanding of the mechanisms of tumorigenesis and progression in GBM, and could be further evaluated for use as cancer biomarkers. Meanwhile, the present study provides an insight into the discovery of potential drug targets for GBM treatment.

Supplementary Data

Funding

This study was supported by The Cooperation Project of Basic Clinical Scientific Research in Capital Medical University [grant nos. 17JL61 and 17JL(TT2X)] and The Scientific Research Common Program of Beijing Municipal Commission of Education (grant no. KM201710025010).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

JJ, DL and LZ conceived and designed the study. JJ and XZ wrote the manuscript. JJ, XZ and YA made substantial contributions to the acquisition of data. JJ, LL and QL analyzed and interpreted the data. LL, DL and LZ reviewed and revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Acknowledgments

Not applicable.

References