Comprehensive Transcriptomic Analysis and Experimental Validation Identify lncRNA HOXA-AS2/miR-184/COL6A2 as the Critical ceRNA Regulation Involved in Low-Grade Glioma Recurrence

Introduction

Glioma is the most common malignant brain tumor in human adults, with a low 5-year overall survival (OS) rate (35%).¹ According to the World Health Organization (WHO) grading system, gliomas are classified from grade I–IV. Grade II and III gliomas are grouped as LGGs, and GBM is regarded as grade IV. GBM, as the most invasive and aggressive glioma, shows the highest mortality despite aggressive treatment.^2,3 Although LGGs have identical morphological phenotypes, the range of patient outcomes is wide.⁴ Because of the aggressive nature of gliomas, a residual tumor after excision is the leading cause of recurrence and disease progression, which is a major reason for the poor prognosis in glioma patients.⁵ Therefore, it is necessary to explore the tumor progression mechanisms and identify potential therapeutic targets in glioma.

Recently, many studies have revealed the distinctions between GBM and LGG by integrating and analyzing multilevel molecular data, including mutation, methylation, and gene expressions. Previous studies have reported that somatic intronic microsatellite loci can differentiate between GBM and LGG.^6–8 Zhang et al revealed distinctions between LGG and GBM by establishing a differential mRNA-lncRNA network using mRNAs and lncRNAs differentially expressed between GBM and LGG.⁹ Therefore, a separate exploration of the mechanisms of LGG and GBM progression will be more appropriate.

Long noncoding RNAs (lncRNAs) and microRNAs are classic noncoding RNAs and have been reported to regulate the development of glioma.^10–12 For example, Li et al found that lncRNA SNHG1 could upregulate FOXP2 and KDM5B by regulating the expression of miR-154-5p and miR-376b-3p, which contributed to the malignant behavior of glioma cells.¹³ Wu et al reported that lncRNA lnc-TALC increased the expression of O⁶-methylguanine-DNA methyltransferase by regulating the c-Met pathway by competitively binding with miR-20b-3p, and this biological process is required for temozolomide resistance and GBM recurrence.¹⁴ To date, a number of public open-access databases have helped uncover more complex protein-coding gene-lncRNA-microRNA regulatory networks. By informatics analysis using TCGA and GEO data, reciprocal regulatory networks among protein-coding genes, lncRNAs and microRNAs have been revealed to participate in cancer progression.^15,16 Lou et al described a miRNA-mRNA regulatory network by using DEPs between GBM samples and normal samples that was involved in tumorigenesis and the development of GBM.¹⁷ Deng et al also revealed hub DEPs by constructing PPI networks of DEPs between primary and recurrent LGG tumor tissues in the TCGA database.¹⁸ Besides, lncRNAs, microRNAs and mRNAs could form complex ceRNA network to promote GBM growth and progression.^19,20 Long et al revealed a complex lncRNA mediated ceRNA network associated with GBM progression by analyzing TCGA and GEO data.²¹ However, the complex regulatory networks and ceRNA networks among protein-coding genes, lncRNAs and microRNAs are still unclear in the recurrence of glioma.

In this study, DEPs, DELs and DEMs were identified, and GO analysis was performed to explore transcriptome variation in the process of LGG or GBM recurrence by cross-database analysis using CGGA and TCGA data. By intersecting the results from analyses of the CGGA and TCGA databases and constructing a PPI network, three hub genes (TIMP1, COL1A1 and COL6A2) were revealed to participate in LGG recurrence. Using the three genes, primary LGGs were clustered into two groups (High_risk and Low_risk), and the High_risk group showed poorer prognosis. Finally, protein-coding gene-lncRNA-microRNA regulatory network analysis and ceRNA network analysis were conducted to reveal hub protein-coding genes, lncRNAs and microRNAs associated with LGG recurrence.

Materials and Methods

Data Collection

We obtained transcriptomic data of GBM and LGG samples from the TCGA (https://portal.gdc.cancer.gov/) and CGGA (http://www.cgga.org.cn/) databases. The details are illustrated in Table 1. We normalized raw counts data from TCGA database by TMM using GDCRNATools package in R software. In the TCGA database, 14 LGG patients had both recurrent and primary samples, and RNA-Seq and miRNA-Seq data were available for a total of 32 samples. In addition, 6 GBM patients had both recurrent and primary samples, and RNA-Seq data were available for a total of 13 samples. RNA-Seq data of grade III glioma cell lines (SW1783, SW1088, HS683 and GOS3) were downloaded from CCLE database (https://portals.broadinstitute.org/ccle). RNA levels of COL6A2, HOXA-AS2 and miR-184 were extracted and log2-transformed.

Table 1 Datasets Used in This Study

Differentially Expressed Protein-Coding Genes, lncRNAs and microRNAs

Differentially expressed RNAs (including DEPs, DELs, and DEMs) were analyzed between primary and recurrent samples of LGG or GBM by using the “limma” package in R software. Because some cases had both primary and recurrent samples in the TCGA database, we identified DEPs, DELs, and DEMs between primary and recurrent tumors in these patients. The DEPs, DELs, and DEMs that met the criteria of P value < 0.05 and |log2 (fold change)| > 1 were selected for further analysis. Volcano plot and heatmap were plotted to show DEPs, DELs, and DEMs by using the “ggplot2” and “pheatmap” package in R software. Venn diagrams were drawn to find overlap using the Bioinformatics & Evolutionary Genomics website (http://bioinformatics.psb.ugent.be/webtools/Venn/).

GO Enrichment Analysis

GO analysis was performed using DEPs with the “clusterProfiler” package in R software. Dotplots were generated to display the most significant GO terms involved in LGG or GBM recurrence (P<0.001). Using the “enrichplot” package in R, circular diagrams were plotted to show the 5 most related GO terms (P<0.001).

Protein–Protein Interaction (PPI) Network Analysis

To better understand the interactive relationships among these enriched genes in GO terms, a PPI network was established by using the STRING database.²² According to the degree calculated using Cytoscape software (version 3.6.1), the top 10 hub genes were identified from the PPI network.

Protein Coding Gene-lncRNA-miRNA Coexpression Network Analysis

The R package “corrplot” was utilized to calculate the correlation coefficients among prognostic DEPs, DELs and DEMs. Then, node pairs with a correlation coefficient > 0.1 were selected for further analysis. A coexpression network containing protein-coding genes, lncRNAs and miRNAs was constructed, and further analysis was conducted to identify a subnetwork by using Cytoscape software (https://cytoscape.org/).

Construction of ceRNA Network

Interactions of lncRNAs with microRNAs and mRNAs were predicted by using the miRcode database. Finally, the lncRNA-miRNA-mRNA ceRNA network analysis was conducted by using Cytoscape software.

Cell Lines and Cell Culture

Two grade III glioma cell line HS683 and SW1088 were purchased from American Type Tissue Culture Collection. The cell lines were maintained in DMEM supplemented with 10% fetal bovine serum at 37°C under a humidified atmosphere of 5% CO₂.

Cell Transfection

miR-184 mimics, miR-184 inhibitor, or negative control were transfected by Lipofectamine3000 (Invitrogen Inc., Carlsbad, CA, USA), which were synthesized by GenePharma (Shanghai, China). The sequences were listed as follows:²³ miR-184 mimic, 5′-UGGACGGAGAACUGAUAAGGGUCCUUAUCAGUUCUCCGUCCAUU-3′; the negative control for miR-184 mimic, 5′-UUCUCCGAACGUGUCACGUTTACGUGACACGUUCGGAGAATT-3′; miR-184 inhibitor, 5′-ACCCUUAUCAGUUCUCCGUCCA-3, negative control for miR-184 inhibitor, 5′-CAGUACUUUUGUGUAGUACAA-3′. siRNAs targeting COL6A2 and HOXA-AS2 were synthesized by GenePharma, and eukaryotic vectors expressing HOXA-AS2 were obtained from Genechem (Shanghai, China). The siRNAs and vectors were used to transfected to glioma cells by Lipofectamine3000. The siRNA sequences were listed as follows:^24,25 siHOXA-AS2, si1: 5′-AAACCUUGUAGAUAGCUUGAGCUGG-3′, si2: 5′-CAAGCUUGACAAGUUCAGCUCAA-3′; siCOL6A2: si1: 5′-GGGCCUCCUUCAUCAAGAATT-3′; si2: 5′-GCAGGCCUGGAUUCAGCUATT-3′.

Western Blotting

Total protein was gathered by RIPA Lysis Buffer. The protein samples were separated by 8% SDS-polyacrylamide gels and transferred to PVDF membranes. After blocking with 5% BSA for 2 h, the membranes were incubated overnight with primary antibodies including mouse anti-GAPDH monoclonal antibody (Catalog number: 60,004-1-Ig, proteintech, Wuhan, China) and rabbit anti-COL6A2 monoclonal antibody (Catalog number: ab180855, Abcam, Cambridge, UK) at 4°C. After being washed in TBST for three times, membranes were incubated with HRP-conjugated anti-rabbit (Catalog number: SA00001-2, proteintech, Wuhan, China) or HRP-conjugated anti-mouse (Catalog number: SA00001-1, proteintech, Wuhan, China) secondary antibodies at a 1:5000 dilution at room temperature for 1h. Then, the binding secondary antibody was detected and recorded using the enhanced ECL Western blotting detection system.

Quantitative Real-Time PCR (qPCR)

The total RNA of glioma cells were extracted with Trizol reagent (Invitrogen, USA) according to the manufacturer’s instruction. Total RNA (1 μg) was used for cDNA synthesis using Reverse Transcription kit (TOYOBO, Osaka, Japan). The primers were: GAPDH, forward: 5′-AACGGATTTGGTCGTATTG-3′, reverse: 5′-GGAAGATGGTGATGGGATT-3′; HOXA-AS2, forward: 5′-GTTCAGCTCAAGTTGAACATA-3′, reverse: 5′-AAACCTTGTAGATAGCTTGAG-3′; COL6A2, forward: 5′- TGCTCCGTGCTCCTGCTCTG-3′, reverse: 5′- ATGGTGACGCTCTCCGAGGTG −3′. The qRT-PCR assays were carried out using SYBR Select Master Mix (Applied Biosystems, Foster City, CA, USA) on ABI 7500 Real-Time PCR system (Applied Biosystems). The PCR conditions were as follows: step 1: 95°C for 5 min; step 2: 40 cycles consisting of 95°C for 10 s, 60°C for 30 s, and 72°C for30 s; step 3: 72°C for 5 min. Relative RNA level of each gene was calculated as 2^−ΔΔCt. Each sample was performed in triplicate.

Cell Counting Kit-8 (CCK-8) Proliferation Assay

The transfected HS683 and SW1088 cells were seeded into 96-well plates at a density of 5000 cells/well. After incubated at 37 °C for 2 h, the absorbance of each well was measured at a wavelength of 450 nm.

Luciferase Assay

Glioma cells were seeded in 24-well plate, and allowed to grow for 24 h. To test the binding specificity of miR-184 in COL6A2, we constructed the wild type reporter vector (COL6A2-3′UTR_WT: 5′-CCCGUCCA-3′) and the corresponding mutant (COL6A2-3′UTR_MUT: 5′-CUUAGUUA-3′). To test the binding specificity of miR-184 in HOXA-AS2, we constructed the wild type reporter vector (HOXA-AS2_WT: 5′-UCCGUCC-3′) and the corresponding mutant (HOXA-AS2_MUT: 5′- GUUAGUU −3′). Then, 100ng pGL3-COL6A2-3′UTR-luciferase plasmid (or MUT) and pGL3-HOXA-AS2-luciferase plasmid (or MUT) was transfected into glioma cells using the Lipofectamine3000 reagent. Subsequently, the cells were transfected with miR-184 mimics, miR-184 inhibitor, or negative control. The luciferase activities in cell lysates were measured at 48 h after transfection by the Dual Luciferase Reporter Assay Kit (Promega, Madison, WI, USA) according the manufacturer’s protocol. Three independent experiments were performed.

Immunohistochemistry (IHC) and Evaluation

Four primary LGG tissues and four recurrent LGG tissues were obtained from Department of Neurosurgery of Shengjing Hospital affiliated to China Medical University and histologically diagnosed with grade II or III glioma in 2017. These experiments were approved by Human Ethics Committee of China Medical University and all patients signed the informed consent. Tissues were fixed by formalin, embedded by paraffin, and finally cut into 4μm sections. The sections were deparaffinized by three 10-min washes in xylene, and then rehydrated by 100%, 95%, 85% and 75% ethanol solutions (5min/time). Antigen retrieval was performed in 0.01 M citrate buffer (pH 6.0) at 120°C for 5 min. Endogenous peroxidase was blocked by 3% H₂O₂ in PBS for 15 min and washed for 2×5 min in PBS. Then, non-specific binding were blocked by 10% bovine serum albumin (BSA) for 45 min. Subsequently, sections were incubated with the primary antibody against COL6A2 (1:100) at 4 °C for overnight. After washed using PBS, the sections were incubated with HRP-conjugated goat anti-rabbit secondary antibody (1:1000) for 45 min. Finally, the sections were stained using 3,3-diaminobenzidine (DAB) for 5 min, counterstained with hematoxylin, mounted and captured representative images using an optical microscope (Olympus, Tokyo, Japan). As for the IHC scoring, protein expression was calculated by multiply the staining intensity scores and the positive cell rate scores. Staining intensity scores: 0, negative; 1, weak; 2, intense; and 3, strongly intense. Positive cell rate scores: 0, 0%; 1, 1–10%; 2, 11–25%; 3, 26–50%; 4, 51–75%; and 5, 76–100%.

Statistical Analysis

The main statistical analysis in this study is stated above. Using a 3-gene expression value, hierarchical cluster analyses were conducted with the “pheatmap” package in R. To explore the prognostic significance of DEMs, ROC curves were plotted, and optimal cutoff values were set to divide all LGG samples as low and high expression groups. The survival curves were plotted using GraphPad Prism 6.0 (La Jolla, CA, United States). When we performed differential expression analysis, only protein-coding genes, lncRNAs or miRNAs with |log2FC| > 1 and P<0.05 were considered statistically significant. Univariate Cox regression analyses were conducted using the “survival” package in R, The Student’s t-test, or the chi-squared test, was used to evaluate significant differences between groups of data. All data are represented as the means ± SD. Student’s t-test using GraphPad Prism 6.0 was used to determine significant differences between two groups. All data are represented as the means and 95% CI. A P value < 0.05 was considered statistically significant.

Results

Recurrence-Related RNA Expression Patterns Vary in GBM and LGG

To explore whether GBM and LGG recurrences contribute to distinct mechanisms, we compared the differentially expressed RNAs between recurrent and primary samples of GBM or LGG, including protein-coding genes, pseudogenes, and long noncoding and short noncoding RNAs. In the CGGA database, there were 183 upregulated and 38 downregulated RNAs in recurrent GBM samples (|log2 (fold change)| > 1, P<0.05, Figure 1A and B, Table S1). A total of 45 upregulated and 7 downregulated RNAs were observed in recurrent LGG samples (|log2 (fold change)| > 1, P<0.05, Figure 1C and D, Table S2). Venn diagrams showed that there were only 6 overlapping RNAs among the upregulated RNAs and 1 overlapping RNA among the downregulated RNAs in recurrent LGG and GBM (Figure 1E and F). These results revealed distinct recurrence mechanisms between GBM and LGG at the transcriptional level.

Figure 1 Recurrent GBM and LGG showed distinct transcriptomic characteristics in the CGGA database. Notes: Volcano plot (A) and heatmap (B) showing differentially expressed RNAs (|log2 (fold change)| > 1, P<0.05) between recurrent and primary GBM samples in the CGGA database (dataset 2). Volcano plot (C) and heatmap (D) showing differentially expressed RNAs (|log2 (fold change)| > 1, P<0.05) between recurrent and primary LGG samples in the CGGA database. (E) Six overlapping genes were upregulated in recurrent GBM and LGG. (F) One overlapping gene was downregulated in recurrent GBM and LGG Abbreviations: GBM_P, primary GBM; GBM_R, recurrent GBM; LGG_P, primary LGG; LGG_R, recurrent LGG.

To validate the above results, the same analyses were conducted in the TCGA database. There were 1332 upregulated and 466 downregulated RNAs in recurrent GBM samples (|log2 (fold change)| > 1, P<0.05, Figure 2A and B, Table S3). A total of 401 upregulated and 616 downregulated RNAs were observed in recurrent LGG samples (|log2 (fold change)| > 1, P<0.05, Figure 2C and D, Table S4). Venn diagrams showed that there were only 6 overlapping RNAs among the upregulated RNAs and 14 overlapping RNAs among the downregulated RNAs in recurrent LGG and GBM (Figure 2E and F).

Figure 2 Transcriptomic characteristics of recurrent GBMs and LGGs in the TCGA database. Notes: Differentially expressed RNAs between primary and recurrent GBM samples are displayed by volcano plot (A, |log2 (fold change)| > 1, P<0.05) and heatmap (B, |log2 (fold change)| > 1, P<0.05). Differentially expressed RNAs between primary and recurrent LGG samples are displayed by volcano plot (C, |log2 (fold change)| > 1, P<0.05) and heatmap (D, |log2 (fold change)| > 1, P<0.05). (E) Six overlapping genes were upregulated in recurrent GBM and LGG. (F) Fourteen overlapping genes were downregulated in recurrent GBM and LGG. Abbreviations: GBM_P, primary GBM; GBM_R, recurrent GBM; LGG_P, primary LGG; LGG_R, recurrent LGG.

Recurrent Mechanisms of GBM and LGG are Different

GO analysis was performed using DEPs. The dotplot of the top 10 enriched GO terms is displayed. In recurrent GBM, there were no identical enriched GO terms between the TCGA and CGGA databases (Figure 3A and B). However, we observed the same enriched GO terms involved in LGG recurrence (Figure 3C and D). The results showed that LGG recurrence was positively associated with GO terms involved in the extracellular matrix (such as collagens, P<0.001).

Figure 3 GO analysis of DEPs in the CGGA and TCGA databases. Notes: GO analysis of DEPs (|log2 (fold change)| > 1, P<0.05) in recurrent GBM in the CGGA (A) and TCGA (B) databases (P<0.001). GO analysis of DEPs (|log2 (fold change)| > 1, P<0.05) in recurrent LGG in the CGGA (C) and TCGA (D) databases (P<0.001). Abbreviations: GBM_P, primary GBM; GBM_R, recurrent GBM; LGG_P, primary LGG; LGG_R, recurrent LGG; BP, biological process; MF, molecular function; CC, cellular component.

GO and PPI Network Analysis Revealed Three Core Genes Involved in LGG Recurrence

To screen critical genes that mediate LGG recurrence, we generated a circular diagram to show the top 5 related GO terms with the enriched genes. Referencing the analysis results from the CGGA and TCGA databases, a total of 20 genes were enriched in GO terms associated with the extracellular matrix (Figure 4A and B). Using these 20 genes to construct a PPI network, the top 10 hub genes were identified (Figure 4C). Among them, three genes were also intersected in both TCGA and CGGA analysis results: TIMP1, COL1A1 and COL6A2. Expression of the three genes was increased more than double in recurrent LGGs in the CGGA (Figure 4D, all P<0.001) and TCGA databases (Figure 4E, all P<0.05).

Figure 4 Identification of TIMP1, COL1A1 and COL6A2 as the core genes involved LGG recurrence. Notes: Circular diagram showed 5 most related GO terms involved in LGG recurrence in CGGA (A) and TCGA (B) database (P<0.001). (C) PPI network analysis showed top10 hub protein-coding genes among genes in GO term associated with extracellular matrix. The colors nodes of the network from light yellow, dark yellow to red indicated gradually crucial members in the network. Of the 10 hub protein-coding genes, TIMP1, COL1A1 and COL6A2 were upregulated in LGG recurrence comparing primary tumor in CGGA (D, all P<0.001) and TCGA (E, all P<0.05) database. *P<0.05, ***P<0.001. Abbreviations: GBM_P, primary GBM; GBM_R, recurrent GBM; LGG_P, primary LGG; LGG_R, recurrent LGG.

To explore the prognostic significance of the 3 genes in primary LGGs, we performed hierarchical cluster analysis. In the TCGA project and CGGA dataset (dataset 2), all primary LGG samples were clustered into two groups: High_risk and Low_risk. The High_risk group showed poorer prognosis than the Low_risk group in the CGGA (dataset 2, HR=2.97 (2.11–5.17), P<0.001, Figure 5A) and TCGA (HR=2.69 (1.86–3.75), P<0.001, Figure 5B) databases. Another two CGGA datasets (dataset 3 and dataset 4) were used as external validation, and the same results were observed (dataset 3, HR=2.20 (1.42–4.09), P=0.001, Figure 5C; dataset 4, HR=1.87 (1.08–3.01), P=0.025, Figure 5D).

Figure 5 A 3-gene signature for overall survival prediction in LGG. Notes: Hierarchical cluster analysis of 3 genes by heatmap showed two clusters: Low_risk and High_risk group. High_risk group with high 3-gene signature showed poorer prognosis than Low_risk group in primary LGGs in CGGA (A, dataset 2) and TCGA (B) database. Using another two datasets in CGGA database, two clusters (Low_risk and High_risk group) were also identified. High_risk group showed poorer prognosis than Low_risk group of primary LGGs (C, dataset 3; D, dataset 4).

Recurrence-Related miRNAs are Distinct Between LGG and GBM

To further explore the roles of miRNAs in LGG and GBM recurrence, we used CGGA microRNA array data for the next analysis. In the CGGA dataset (dataset 6), 29 increased and 59 decreased miRNAs were observed in recurrent GBM (|log2 (fold change)| > 1, P<0.05, Figure 6A and B, Table S5). In LGG, we found 9 upregulated and 8 downregulated miRNAs in recurrent samples (|log2 (fold change)| > 1, P<0.05, Figure 6C and D, Table S6). Venn diagrams showed only 2 overlapping miRNAs among the upregulated miRNAs and 1 overlapping miRNA among the downregulated miRNAs in recurrent LGG and GBM (Figure 6E and F).

Figure 6 Recurrent GBMs and LGGs showed different changed microRNAs in CGGA database. Notes: Volcano plot (A, |log2 (fold change)| > 1, P<0.05) and heatmap (B, |log2 (fold change)| > 1, P<0.05) showing DEMs between recurrent and primary GBM samples in CGGA database (dataset 6). Volcano plot (C, |log2 (fold change)| > 1, P<0.05) and heatmap (D, |log2 (fold change)| > 1, P<0.05) showing DEMs between recurrent and primary LGG samples in CGGA database (dataset 6). (E) Two overlapped microRNAs were up-regulated in recurrent GBM and LGG. (F) One overlapped microRNAs was down-regulated in recurrent GBM and LGG. Abbreviations: GBM_P, primary GBM; GBM_R, recurrent GBM; LGG_P, primary LGG; LGG_R, recurrent LGG.

Prognostic miRNAs and lncRNAs in LGG

We further screened prognostic miRNAs in LGG patients by using TCGA miRNA-Seq data. First, we compared differentially expressed miRNAs between recurrent and primary LGG samples from 14 LGG patients. There were 23 upregulated and 20 downregulated miRNAs in recurrent LGGs (Figure 7A and B, Table S7). Among them, 8 of 23 upregulated miRNAs were significantly associated with inferior prognosis, and 11 of 20 downregulated miRNAs were positively associated with better prognosis (Figure 7C and Figure S1, P<0.01).

Figure 7 Identification of recurrence-related microRNAs in LGGs in TCGA database. Notes: Volcano plot (A, |log2 (fold change)| > 1, P<0.05) and heatmap (B, |log2 (fold change)| > 1, P<0.05) showing DEMs between recurrent and primary GBM samples in TGGA database. (C) Among DEMs, prognosis analysis showed that 6 downregulated miRNAs were positively associated with better prognosis (HR <0.5, P<0.01) and 5 upregulated miRNAs were significantly associated with inferior prognosis (HR >2, P<0.01). Abbreviations: GBM_P, primary GBM; GBM_R, recurrent GBM; LGG_P, primary LGG; LGG_R, recurrent LGG.

According to the above results on differentially expressed RNAs, we identified 35 upregulated lncRNAs and 71 downregulated lncRNAs in recurrent LGGs. By Cox univariate regression analysis, 31 of 35 upregulated lncRNAs were significantly related to worse prognosis, and 60 of 71 downregulated lncRNAs were positively associated with better prognosis (Table S8, P<0.01).

Identification and Validation of LncRNA HOXA-AS2/miR-184/COL6A2 ceRNA Regulation Associated with LGG Recurrence

The TCGA database supplies multidimensional data, which includes somatic mutations, copy number variants (CNVs), gene expression, lncRNA expression, and miRNA expression. Therefore, we could establish a protein-coding gene-lncRNA-miRNA coexpression network using these data. There were 12 protein-coding gens enriched in the GO terms involved in the extracellular matrix, 19 prognostic miRNAs and 91 prognostic lncRNAs in the TCGA database. Using the protein-coding genes and prognostic lncRNAs and miRNAs, the node pairs with correlation coefficients > 0.1 and P value < 0.05 were selected for further analysis. Accordingly, we constituted a complex co-expression network containing 12 protein coding genes, 18 miRNAs and 29 lncRNAs (Figure 8A). To better understand whether the above 29 hub lncRNAs regulated 12 protein coding genes via ceRNA regulatory mechanism, we performed ceRNA network analysis by using miRcode database. Finally, we constructed a ceRNA network including 9 protein coding genes, 6 miRNAs and 12 lncRNAs. As shown in Figure 8B, HOXA-AS2 may act as a competing endogenous RNA (ceRNA) for miR-184 to regulate the hub gene COL6A2 expression. COL6A2 expression was positively related to LncRNA HOXA-AS2 expression, and negatively related to miR-184 expression in glioma in TCGA database (Figure 8C). In CCLE database, the four grade III glioma cell lines showed different levels of expression of LncRNA HOXA-AS2, COL6A2 and miR-184 (Figure 8D). Because HS683 showed high expression of them and SW-1088 showed low levels of them, we further selected these two cell lines for the next experiments. We also detected the COL6A2 protein expression by IHC, and COL6A2 expression was significantly higher in recurrent LGG tissues (n=4) than in primary LGG tissues (n=4) (P<0.01, Figure 8E). By referencing to miRcode database, we obtained two potential binding site of miR-184 in LncRNA HOXA-AS2 mRNA (Figure 8F). In order to further validate the interaction, luciferase reporter assays were performed. The miR-184 inhibitors enhanced the luciferase activity of HOXA-AS2_WT reporter systems in SW1088, but not HOXA-AS2_MUT (Figure 8G). The miR-184 mimics inhibited the luciferase activity of HOXA-AS2_WT reporter systems in HS683, but not HOXA-AS2_MUT (Figure 8G). Furthermore, we also obtained the binding site of miR-184 in 3′UTR of COL6A2 mRNA (Figure 8H). Similarly, the miR-184 inhibitors enhanced the luciferase activity of COL6A2-3′UTR_WT reporter systems in SW1088, but not COL6A2-3′UTR_MUT (Figure 8I). The miR-184 mimics inhibited the luciferase activity of COL6A2-3′UTR_WT reporter systems in HS683, but not COL6A2-3′UTR_MUT (Figure 8I).

Figure 8 Identification and validation of HOXA-AS2/miR-184/COL6A2 ceRNA regulation in LGGs. Notes: (A) Construction of co-expression networks analysis using 12 protein-coding genes, 18 miRNAs and 29 lncRNAs. The node pairs with correlation coefficients > 0.1 and P value < 0.05. Blue circle: miRNA; Red sexangle: genes coding protein; Green diamond: LncRNA. (B) ceRNA network analysis among protein-coding genes, lncRNAs and miRNAs in the above co-expression network. Green circle: miRNA related to good prognosis; Red circle: miRNA related to poor prognosis; Blue sexangle: genes coding protein; Green diamond: LncRNA involved in good prognosis; Red diamond: LncRNA involved in poor prognosis. (C) Correlation analysis between COL6A2 and HOXA-AS2 or miR-184 using TCGA data (n=525). (D) The mRNA levels of COL6A2, HOXA-AS2 and miR-184 in grade III glioma cell lines in CCLE database (SW1783, SW1088, HS683 and GOS3). (E) COL6A2 protein level was significantly higher in recurrent LGG tissues (n=4) than in primary LGG tissues (n=4) by IHC (200×, P<0.01). (F) The putative binding sites by miRcode and corresponding mutant site for HOXA-AS2 within miR-184. (G) The miR-184 inhibitors enhanced the luciferase activity of HOXA-AS2_WT in SW-1088 cells by luciferase reporter assay, but not HOXA-AS2_MUT. The miR-184 mimics reduced the luciferase activity of HOXA-AS2_WT in HS683 cells by luciferase reporter assay, but not HOXA-AS2_MUT. (H) Schematic representation of the mature miR-184 targeting site in the 3′-UTR predicted by miRcode and mutant site of COL6A2 mRNA. (I) The miR-184 inhibitors enhanced the luciferase activity of COL6A2-3′UTR_WT in SW1088 cells by luciferase reporter assay, but not COL6A2-3′UTR_MUT. The miR-184 mimics reduced the luciferase activity of COL6A2-3′UTR_WT in HS683 cells by luciferase reporter assay, but not COL6A2-3′UTR_MUT. **P < 0.01, ***P < 0.001.

LncRNA HOXA-AS2/miR-184/COL6A2 ceRNA Regulation Promoted Glioma Cell Proliferation

In order to explore whether LncRNA HOXA-AS2/miR-184/COL6A2 ceRNA regulation involved in glioma cell proliferation, rescue experiments were performed. We first knocked down (Figure S2A) or overexpressed HOXA-AS2 (Figure S2B), and knocked down COL6A2 (Figure S2C) in glioma cells. The more effective siRNA sequence was selected for the further experiments. HOXA-AS2 upregulation increased COL6A2 expression in SW1088 cells, which was partly reversed by miR-184 mimics (Figure 9A). HOXA-AS2 downregulation suppressed COL6A2 expression in HS683, which was restored by miR-184 inhibitors to some extent (Figure 9A). Forcing expression of HOXA-AS2 enhanced cell proliferation ability of SW1088 cells, and the effect was attenuated by miR-184 mimics (Figure 9B). In HS683 cells, HOXA-AS2 knockdown mediated cell proliferation inhibition was reversed by miR-184 inhibitors (Figure 9B). HOXA-AS2-mediated COL6A2 upregulation could be partly blocked by COL6A2 knockdown in SW1088 cells (Figure 9C). Double knockdown of HOXA-AS2 and COL12A1 more apparently suppressed expression of COL6A2 in HS683 cells than single knockdown of HOXA-AS2 and COL6A2 (Figure 9C). Cell proliferation ability of SW1088 was enhanced by HOXA-AS2, which was attenuated by COL6A2 downregulation (Figure 9D). Double knockdown of HOXA-AS2 and COL6A2 more apparently inhibited cell proliferation in HS683 than single knockdown of them (Figure 9D).

Figure 9 HOXA-AS2/miR-184/COL6A2 ceRNA regulation promotes glioma cell proliferation. Notes: (A) Overexpression of HOXA-AS2 increased COL6A2 protein level in SW1088 cells, which was blocked by miR-184 mimics. Knockdown of HOXA-AS2 decreased COL6A2 protein level in HS683 cells, which was partly reversed by miR-184 inhibitors. (B) HOXA-AS2-mediated cell proliferation was attenuated by miR-184 mimics in SW1088 cells. HOXA-AS2 knockdown inhibited cell proliferation in HS683 cells, which was partly reversed by miR-184 inhibitors. (C) Overexpression of HOXA-AS2 increased COL6A2 RNA and protein level in SW1088 cells, which was blocked by COL6A2 knockdown. Double knockdown of HOXA-AS2 and COL12A1 more apparently decreased COL6A2 RNA and protein level in HS683 cells than single knockdown of HOXA-AS2 and COL6A2. (D) HOXA-AS2-mediated cell proliferation was attenuated by COL6A2 downregulation in SW1088 cells. Double knockdown of HOXA-AS2 and COL12A1 more apparently inhibited cell proliferation of HS683 cells than single knockdown of HOXA-AS2 and COL6A2. *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion

Although glioma patients receive aggressive treatment including surgery, chemotherapy and radiotherapy, some patients suffer tumor recurrence.^26,27 Therefore, exploring the recurrence mechanisms of glioma will be conducive to improving the diagnosis and treatment of glioma. An increasing number of studies have inferred the distinctions between GBM and LGG.^6–9 Hence, we explored the recurrence mechanisms of GBM and LGG separately. By comparing them to primary samples, we found that recurrent samples of GBM and LGG showed distinct changes in transcriptomic characteristics due to less overlapping DEPs, DELs and DEMs. Further GO analysis using DEPs also indicated that the recurrent mechanisms of GBM were different from those of LGG. However, in GBM, we did not obtain similar results of GO analysis from the CGGA and TCGA databases, which may be due to the large heterogeneity of GBM.²⁸ For example, the markedly enriched GO terms (MF) in recurrent GBMs were protein heterodimerization activity and structural constituent of ribosome in CGGA database and channel activity and passive transmembrane transporter activity in the TCGA database, although these GO terms have been reported to contribute to GBM progression.^28–31 Glioma stem cells play critical roles in glioblastoma multiforme chemotherapy resistance and recurrence, such as CD133 cells.³² Glioma stem cells can self-renew and differentiate into other specific GBM subpopulations, which is involved in activation of Wnt/β-catenin pathway.³³ Maintenance of stem cell characteristics may be involved in activation of the above pathway related GBM recurrence.^33,34 Therefore, they have been held responsible for malignant relapse after primary standard therapy and the poor prognosis of recurrent GBM. Nevertheless, the enriched GO terms in recurrent LGG in the CGGA database were similar to those in the TCGA database, which were involved in the extracellular matrix (mainly including collagens). The extracellular matrix plays pivotal roles in multiple biological processes in cancers, including metastasis, immune evasion and drug resistance.^35–37

Of the genes in the top 5 most relevant GO terms involved in LGG recurrence, TIMP1, COL1A1 and COL6A2 were upregulated in recurrent samples in both the CGGA and TCGA databases. In pancreatic cancer, high expression of TIMP1 can attenuate the gemcitabine therapy response and accelerate tumor growth and liver metastasis.³⁸ COL1A1 is a member of the collagen I family and has been reported to be associated with tumor cell proliferation and invasion in many cancers, such as breast, lung and renal cancers.³⁹ As a member of the Collagen VI family, COL6A2 is widely expressed in a variety of cancers and favors cancer progression.⁴⁰ These three core genes were found to interact with each other and with other genes involved in the extracellular matrix as demonstrated using the STRING database. Furthermore, hierarchical cluster analysis using the three genes indicated that there were two clusters of LGGs, and the High_risk group with high expression of these genes showed poor prognosis in primary LGGs. These results indicated the critical roles of these genes in LGG recurrence. These three genes have been reported to be significantly associated with poor prognosis in other cancers, such as bladder cancer and colorectal cancer.^41,42 Many previous studies have been designed to acquire prognostic signatures of microRNAs, pseudogenes and lncRNAs, and the studies screened the genes by prognostic analyses.^43–45 However, in this study, we screened the hub genes by analyzing DEPs between recurrent and primary samples. Deng et al also screened DEPs between recurrent and primary samples using TCGA data and identified hub genes (APOA2, COL3A1, COL1A1, TYR, COL1A2, COL5A1, PAPOLB, IGF2BP1 and ANHX) by PPI analysis, but did not intersect these data with those from other databases.¹⁸ Our screened hub genes were validated in another dataset. Because some cases had both primary and recurrent lesions in the TCGA database, in order to obtain more reliable results, we analyzed DEPs, DELs and DEMs using mRNA-Seq and miRNA-Seq data from these cases. LGG patients were Chinese in CGGA and Westerners in TCGA. In this study, TIMP1, COL1A1 and COL6A2 were considered critical genes in the process of LGG recurrence and were identical in different geographic locations.

It is well known that cancer progression is orchestrated by a complex regulatory network.^46,47 Recently, many studies have reported that miRNAs and lncRNAs are involved in gene expression regulation. Song et al revealed a regulatory network of CLDN4 in gastric cancer, which was mediated by miRNAs and lncRNAs, resulting in cancer cell proliferation and invasion.⁴⁸ In GBM, AGAP2-AS1 epigenetically decreases TFPI2 expression by binding to EZH2 and LSD1, and this biological process promotes cancer cell proliferation and invasion.⁴⁹ Therefore, exploring complex regulatory networks among protein-coding genes, lncRNAs and miRNAs will contribute to clarifying the recurrence mechanisms of cancer. In this study, we constructed a coexpression networks among 12 protein coding genes, 12 miRNAs and 29 lncRNAs. Of the lncRNAs, HOXA10-AS was upregulated in LUAD tissues and correlated with poor prognosis in lung cancer patients, and HOXA-AS2 could promote cancer progression via complicated regulatory mechanisms, including induction of EMT, direct inhibition of gene expression and methylation.^50,51 For the screened miRNAs, miR-139-3p was identified as a key miRNA related to hepatocellular carcinoma, and miR-326 could inhibit cell proliferation, migration, and invasion by targeting oncogenic CCND1.^52,53 In skin fibroblasts, lncRNA AC067945.2 overexpression inhibited the expression of COL1A1, COL1A2, COL3A1.⁵⁴ Cao et al demonstrated that lncRNA PVT1 promotes atrial fibrillation by increasing collagen expression via the miR-128-3p-SP1-TGF-β1-Smad axis in atrial fibroblasts.⁵⁵ In this study, we identified a complex regulatory network of extracellular matrix-related genes, lncRNAs and miRNAs in the process of LGG recurrence.

The ceRNA network is also a critical regulatory mechanism between lncRNAs and Mrna^56,57 We identified a complex ceRNA network among 9 protein coding genes, 6 miRNAs and 12 lncRNAs. The results indicated that HOXA-AS2 could act as a competing endogenous RNA (ceRNA) for miR-184 to regulate the hub gene COL6A2 expression. In osteosarcoma cells, LncRNA HOXA-AS2 could promote cancer cell migration and invasion by acting as a ceRNA of miR-520c-3p.⁵⁸ Although the research is not innovative enough, we confirmed that LncRNA HOXA-AS2 sponge miR-184 with complementary binding sites to increase expression of COL6A2, which resulted in LGG proliferation.

Conclusion

We demonstrated that the recurrence mechanisms of GBM and LGG are different at the transcriptional level, including protein coding genes, lncRNAs and miRNAs. Further analysis revealed three critical genes, TIMP1, COL1A1 and COL6A2, and the 3-gene signature was identified to improve prognosis prediction in LGG. Finally, LncRNA HOXA-AS2 acted as a ceRNA for miR-184 to regulate COL6A2 expression, and this ceRNA regulation induced cell proliferation in LGG.