All microarray expression data, including RNA-seq and miRNA-seq, and their correlated clinical data were download from The Cancer Genome Atlas (20) (TCGA) database (LGG and GBM datasets). After invalid information was discarded, 670 RNA-seq data of glioma patients with clinical information were finally obtained. RNA-seq data of level 3 HTSeq-FPKM were transformed into TPM (transcripts per million reads) format for subsequent analysis. The clinical data which were “Unavailable” or “Unknown” were considered as missing values. The prognosis data, such as overall survival (OS), progression-free interval (PFI) and disease-specific survival (DSS), were obtained from the study of Liu (21). The RNA-seq data of normal brain samples from standardized Genotype-Tissue Expression (GTEx) (22) were considered as normal controls and were obtained from the Xena (23) database (https://xenabrowser.net/).

The data were divided into high and low REST expression groups, and the expression profile data of the top 30% (0–30%) and the bottom 30% (70–100%) were extracted for subsequent analysis. The DESeq2 package (24) was used to identify the DEGs of glioma. HTSeq counts were used for the analysis of DE-mRNAs and DE-lncRNAs, while miRNA-seq counts were used for the analysis of DE-miRNAs. |logFC| > 1.5 and adjusted p < 0.05 was set as the cut-off criteria for DEGs, and the results were visualized with volcano plot and heatmap.

The Cluster Profiler package (25) was used for functional enrichment analysis. To explore the possible functions, annotations of the REST-related DEGs, GO enrichment analysis and KEGG pathway enrichment analysis were conducted and visualized.

Differences in REST protein expression between glioma and normal brain tissues were examined by immunohistochemistry (IHC) images from the HPA database (26). The staining intensity evaluation was determined from the HPA (https://www.proteinatlas.org).

Starbase (27) (http://starbase.sysu.edu.cn/) was used to explore the regulatory interaction between REST and miRNA or lncRNA, and the lncRNA-miRNA-REST ceRNA network was constructed with Cytoscape software (28).

Human cell lines U251, T98G, and 293T were purchased from ATCC and cultured in DMEM medium (Invitrogen, Thermo Fisher Scientific, USA). Culture media contained 10% fetal bovine serum (FBS) (Gibco) plus penicillin Streptomycin (Gibco), GlutaMAX(Gibco), and MEM non-Essential Amino Acids (MEM-NEAA) (Gibco). The above cells were all incubated in an incubator at 37°C with 5% CO2.

The mimics of miRNAs and siRNA were obtained from GenePharma Co. Ltd (Shanghai, China). The primers of all the genes were synthesized from TsingKe Biological Technology (Beijing, China). All the sequences were displayed in Supplementary Table 1. Transfection of miRNAs and siRNA was performed using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) when the cell density reached approximately 50%-70%(usually 24 h after inoculation) following instruction of manufacture, transfection of plasmid was performed using Hieff TransTM Liposomal Transfection Reagent (YEASEN, Shanghai, China), according to the manufacturer's instructions. The concentrations of the siRNA and miRNA used in transfection was 5 pmol/well for 96-well-plate (40 pmol/well for 12-well-plate and 100 pmol/well for 6-well-plate, respectively) following the manufacturer's instructions (GenePharma Co. Ltd, Shanghai, China).

The methods of RNA extraction and real-time RT-PCR analysis were performed as previously described (29). Briefly, total RNA was extracted and converted into complementary DNA (cDNA). Quantitative PCR (qPCR) was performed using QuantiNova SYBR Green PCR Kit (QIAGEN, Hilden, Germany). Quantifications were normalized and calculated by the 2^−ΔΔCt method.

Western blot was performed as previously (29). Briefly, total protein was quantified using bicinchoninic acid (BCA) method and separated using 10% SDS-PAGE. Then the gels were transferred to PVDF membranes. After blocked with 5% skim milk in TBST, the membranes were incubated with the primary antibodies: GAPDH (1:6000; Proteintech, USA) and REST (1:2000; Proteintech, USA) overnight at 4°C. After washing with TBST, the membranes were reacted with the secondary antibody. Finally, enhanced chemiluminescence (ECL) kit (MD Millipore, Germany) was used to measure the signals in the Bio-Rad ChemiDocXRS Imaging system (Bio-Rad, Hercules, CA, USA) and calculated with Image J software.

The proliferative capacity of glioma cells was determined using the CCK8 kit (DOJINDO, Japan). After cell transfection, log phase grown cells were implanted in 96-well-plates with 3,000 cells per well. At 0, 24, 48, and 72 h post-transfection, 10 μl CCK-8 reaction solution was added to each well. The plates were then placed back in the incubator for 1 h. Finally, the absorbance of each well was measured at 450 nm using a multifunctional enzyme marker.

After cell transfection, cells were digested and resuspended with serum-free DMEM culture medium, inoculated into 6-well-plates at a cell count of 150,000 per well, and cell transfection was performed at 24 h. When the cells were grown to fusion (usually 48 h after inoculation), a 200 ul pipette tip was used to make a vertical “wound” in the cell layer. This was followed by incubation with DMEM containing 2% FBS to eliminate the effect of proliferation on the migratory phenotype. Scratch pictures were collected at 0, 24, 48, and 96 h (until the scratch healed in the experimental or control group). Image J software was applied to calculate the scratch healing area. Scratch healing ratio of X h = (0 h scratch area–X h scratch area)/0 h scratch area.

The methods of luciferase assay were performed as previously described (29). Direct target sites were predicted based on TargetScan (http://www.targetscan.org/vert_72/) and starBase v2.0 (http://starbase.sysu.edu.cn/) (27). The luciferase plasmid was constructed, and the 3′-UTR sequences of REST and lncRNAs were cloned into psiCHECK-2 promoter vector (Promega) by Ruibio Biotech (Beijing, China).

All statistical analyses were performed with R (v3.6.2) or GraphPad Prism 8.0 (GraphPad Inc., California, CA, USA). Wilcoxon rank sum test was used to analyze the expression of glioma and normal samples. Receiver operating characteristic (ROC) and the area under curve (AUC) were used to test the diagnostic value of REST expression for glioma. Kruskal–Wallis test and Wilcoxon signed-rank test were used to analyze the correlation between REST expression and clinical characteristics. Spearman rank correlation was used to analyze the correlation between REST expression and hub genes in ceRNA network. Cox regression analysis or Kaplan-Meier analysis were used to evaluate the correlation between REST expression and clinical prognosis. In Cox regression analysis, variables with p < 0.1 in univariate Cox regression were included in multivariate Cox regression analysis. Spearman correlation test was used to analyze the correlation of gene expression. P < 0.05 were considered statistically significant.

The flow chart of our study is shown in Figure 1. To analyze the difference in REST expression between gliomas and normal tissues, we used the TCGA database combined with GTEX data to analyze the gene expression level of REST in pan-cancer and adjacent tissues. The results showed that the expression of REST was significantly up-regulated in most tumors compared to normal tissues, such as esophageal carcinoma, acute myeloid leukemia (LAML) and testicular germ cell tumors (TGCT) (all p < 0.001, Figure 2A). Furthermore, we observed that REST expression level in glioma (LGG + GBM) was significantly higher than that in normal samples (Figure 2B). The ROC curve was used to analyze the diagnostic value of REST. The AUC was 0.9, which indicates that REST is a potential diagnostic biomarker (Figure 2C). In addition, we explored the protein expression of REST in glioma and normal samples and examined the immunohistochemical results of REST protein using the HPA (Supplementary Figure 1). The representative results showed that REST was not expressed in normal brain tissue, whereas high expression of REST was observed in glioma tissue (Figure 2D).

The Kruskal–Wallis test and Wilcoxon signed-rank test were used to analyze the correlation between REST expression and clinical characteristics of glioma. We found that higher REST expression level positively correlated with higher WHO grade (p < 0.001, Figure 3A), IDH status-WT (p < 0.001, Figure 3B), non-1p/19q codeletion (p < 0.001, Figure 3C), histological type-GBM (p < 0.001, Figure 3D), tumor status -with tumor (p < 0.001, Figure 3E) and primary therapy outcome-PD (p < 0.001, Figure 3F). Consistent results were obtained from chi-square test (Table 1). Logistic regression analysis also suggested similar results. The results showed that REST expression was significantly correlated with 1p/19q codeletion (non-codeletion) (OR = 28.73, p < 0.001), IDH status (WT) (OR = 4.60, p < 0.001), WHO grade (G3–4) (OR = 3.84, p < 0.001), histological type (GBM) (OR = 3.44, p < 0.001), tumor status (with tumor) (OR = 2.5, p < 0.001), primary therapy outcome (SD + PD) (OR = 2.09, p < 0.001) and Karnofsky score (≤80) (OR = 2.01, p < 0.001, Table 2). These results suggest that REST expression is closely related to worse clinical characteristics.

To explore the relationship between REST expression and prognosis of glioma patients, we analyzed the correlation of REST expression with OS, PFI and DSS. The results showed that higher REST expression was associated with worse OS (HR = 2.69; 95% CI, 2.06–3.52; Figure 3G), worse PFI (HR = 2.33; 95% CI, 1.85–2.94; p < 0.001; Figure 3H) and worse DSS (HR = 2.97; 95% CI, 2.22–3.95; p < 0.001; Figure 3I). Univariate Cox regression analysis of OS showed that WHO grade (G3–4), IDH status (WT), 1p/19q codeletion (non-codeletion), histological type (GBM), tumor status (with tumor), primary therapy outcome (SD + PD), age (>45), Karnofsky score (≤80) and high REST expression level was associated with poor OS (all p < 0.001, Table 3). Further multivariate Cox regression analysis suggested that high REST expression level were an marginally independent prognostic marker of OS.

To further explore the biological function of REST in glioma, we grouped the expression profiles according to the level of REST expression. We analyzed differences in whole-gene expression profile between samples with REST expression in the top 30% (0–30%) and those in the bottom 30% (70–100%), and found that 2,336 up-regulated genes and 1,599 down-regulated genes were significantly correlated with REST expression (|logFC| > 1.5 and adjusted p < 0.05, Figure 4A). The top 20 up-regulated and the top 20 down-regulated DE-mRNAs are shown in a heat map (Figure 4B). In addition, we conducted GO and KEGG enrichment analysis of the REST-associated DEGs. The results of GO analysis showed that changes in biological processes (BP) of REST-associated DEGs were mainly enriched in extracellular matrix (ECM) organization, extracellular structure organization and neurotransmitter transport (Figure 4C, Table 4), while changes in cellular components (CC) were mainly enriched in the collagen-containing ECM, synaptic membrane and pre-synapse (Figure 4D, Table 4). Changes in molecular function (MF) were mainly enriched in passive transmembrane transporter activity, channel activity and ECM structural constituents (Figure 4E, Table 4). KEGG pathway analysis revealed that the REST-associated DEGs were mainly enriched in neuroactive ligand-receptor interaction, ECM-receptor interaction, and nicotine addiction (Figure 4F, Table 4).

To study the lncRNA and miRNA associated with REST expression, we analyzed the DEGs of the lncRNA profile and miRNA profile based on the REST high-low expression group. In the lncRNA DEGs, there were 508 up-regulated lncRNAs and 544 down-regulated lncRNAs (|logFC| > 1.5 and adjusted p < 0.05, Figure 5A). In the miRNA DEGs, there were 36 up-regulated miRNAs and 27 down-regulated miRNAs (|logFC| > 1.5 and adjusted p < 0.05, Figure 5B). The top 20 up-regulated and the top 20 down-regulated DE-lncRNAs and DE-miRNAs are shown with a heat map (Figures 5C,D).

In recent years, the ceRNA hypothesis has suggested a new mechanism for RNA regulatory interactions (30). miRNA can silence genes by binding to the 3′UTR of mRNA, while ceRNA can competitively sponge miRNA to up-regulate the expression of its downstream genes. ceRNA and mRNA expression are positively correlated, while miRNA and mRNA expression are negatively correlated. To explore the REST-associated ceRNA regulatory network, we constructed a lncRNA-miRNA-REST network. Starbase (27) was used to predict miRNA that may target REST, which led to the discovery of 165 miRNAs; of which, 4 miRNAs were in the list of REST-associated down-regulated miRNAs (Figure 5E). These 4 miRNAs were used to predict 321 target lncRNAs; among which, 8 lncRNAs were in the list of REST-associated up-regulated lncRNAs (Figure 5F). We then constructed a lncRNA-miRNA-REST ceRNA network including predicted miRNAs and lncRNAs (Figure 5G). The correlation between miRNAs or lncRNAs obtained from the above analysis and REST expression was analyzed individually, and the results showed that the expression of NR2F2-AS1 (R = 0.307, p < 0.001), HOTAIRM1 (R = 0.484, p < 0.001), CYTOR (R = 0.371, p < 0.001), MIR4435-2HG (R = 0.420, p < 0.001), HAND2-AS1 (R = 0.430, p < 0.001), PVT1 (R = 0.463, p < 0.001), AC009502.1 (R = 0.285, p <0.001), and AC0098864.1 (R = 0.119, p = 0.002) were positively correlated with REST expression, while the expression of hsa-miR-137 (R = −0.459, p < 0.001), hsa-miR-129-1-3p (R = −0.282, p < 0.001), hsa-miR-876-5p (R = −0.261, p < 0.001), and hsa-miR-873-5p (R = −0.223, p < 0.001) were negatively correlated with REST expression (Figure 6).

In order to determine the effectiveness of predicted miRNAs and its regulatory relationship with REST, we first overexpressed 4 miRNAs with miRNA-mimics and observed their own expression efficiency. The results showed that the level of miRNAs-mimics were significantly increased in U251 and T98G cell lines (Figures 7A,B). CCK8 assay was used to investigate the role of miRNAs in glioma cell growth and proliferation. The results showed that miR-129 and miR-137 can significantly inhibited the proliferation in glioma cell lines (Figures 7C,D). Wound healing assay was performed to evaluate the migration ability of miRNAs, and the results showed that all four miRNAs could significantly inhibit the migration of glioma cell lines (Figures 7E,F). Then, the silencing efficiency of REST mRNA (Figures 7G,H) and protein (Figure 7I) were determined by qRT-PCR and Western Blot. The results showed that only miR-129 and miR-137 can effectively regulate the expression levels of REST, at both mRNA and protein levels. Therefore, we selected miR-129 and miR-137, which can effectively regulate REST expression and have tumor inhibitory effects, for our further research. Then luciferase assays were performed to evaluate the direct binding of miR-129 or miR-137 with REST. The result showed that both miR-129 and miR-137 can directly bind to REST 3'UTR (Figures 7J,K).

Further, LncRNA NR2F2-AS1 and HOTAIRM1, which were predicted to be associated with miR-129 or miR-137, were selected for further experimental verification. First, we transfected U251 cells with siRNAs of LncRNAs to knockdown the expression of LncRNAs (Figure 8A). Then, qPCR and Western Blot were performed to evaluate the regulatory effect of REST. Results showed that both two lncRNAs can reduce the expression levels of REST mRNA (Figure 8B) and protein (Figure 8C). Next, luciferase assay was performed to evaluate the direct binding of miRNA and lncRNA. It was found that miR-129 can directly bind to NR2F2AS1 (Figure 8D), and miR-137 can directly bind to HOTAIRM1 (Figure 8E). Meanwhile, in TCGA-GBM, it can be found that the expression levels of NR2F2AS1 and HOTAIRM1 were significantly higher than that of normal tissues (Figures 8F,G), and the high expression levels were significantly correlated with poor prognosis in TCGA-GBM database (Figures 8H,I). Finally, the mechanistic scheme of lncRNA-miRNA-REST in glioma is shown in Figure 8J.