Mesenchymal glioblastoma constitutes a major ceRNA signature in the TGF-β pathway.

Rationale: Competitive endogenous RNA (ceRNA) networks play important roles in posttranscriptional regulation. Their dysregulation is common in cancer. However, ceRNA signatures have been poorly examined in the invasive and aggressive phenotypes of mesenchymal glioblastoma (GBM). This study aims to characterize mesenchymal glioblastoma at the mRNA-miRNA level and identify the mRNAs in ceRNA networks (micNET) markers and their mechanisms in tumorigenesis.
Methods: The mRNAs in ceRNA networks (micNETs) of glioblastoma were investigated by constructing a GBM ceRNA network followed by integration with a STRING protein interaction network. The prognostic micNET markers of mesenchymal GBM were identified and validated across multiple datasets. ceRNA interactions were identified between micNETs and miR181 family members. LY2109761, an inhibitor of TGFBR2, demonstrated tumor-suppressive effects on both primary cultured cells and a patient-derived xenograft intracranial model.
Results: We characterized mesenchymal glioblastoma at the mRNA-miRNA level and reported a ceRNA network that could separate the mesenchymal subtype from other subtypes. Six genes (TGFBR2, RUNX1, PPARG, ACSL1, GIT2 and RAP1B) that interacted with each other in both a ceRNA-related manner and in terms of their protein functions were identified as markers of the mesenchymal subtype. The coding sequence (CDS) and 3'-untranslated region (UTR) of TGFBR2 upregulated the expression of these genes, whereas TGFBR2 inhibition by siRNA or miR-181a/d suppressed their expression levels. Furthermore, mesenchymal subtype-related genes and the invasion phenotype could be reversed by suppressing the six mesenchymal marker genes. Conclusions: This study suggests that the micNETs may have translational significance in the diagnosis of mesenchymal GBM and may be novel therapeutic targets.

Introduction

Competitive endogenous RNAs (ceRNAs) were recently reported to play a role in posttranscriptional regulation. Recent findings indicate that RNA transcripts can communicate with each other via common microRNA (miRNA) response elements (MREs) 1. This hypothesis provides a powerful explanation for the inconsistencies in RNA levels compared to protein levels. This new biological language has given rise to novel insights on human development and disease. The muscle-specific transcription factors (TFs) MAML1 and MEF2C are upregulated by the long noncoding RNA linc-MD1 in a ceRNA-related manner, leading to muscle differentiation 2. Dysregulation of the ceRNA network has been reported in several cancers 3-5. PTENP1, a pseudogene of PTEN, functions as a tumor suppressor by upregulating PTEN in a ceRNA-related manner 6, 7. Further studies are needed beyond profiling the imbalance of ceRNA networks in cancers to investigate the role of ceRNAs in tumor progression. Regulation of established oncogenic pathways in glioblastoma by miRNA-RNA interactions was first described by Sumazin et al., whom reported that PTEN is a core-hub regulated by 13 miRNAs in glioblastoma multiform (GBM) 5.

GBM is the most common malignant primary brain tumor in adults and has a median survival of 12.2 to 18.2 months despite tumor resection followed by radiation and chemotherapy 8, 9. Temozolomide (TMZ) is a widely used and virtually unique chemotherapeutic drug for GBM 10. Molecular alteration profiles of GBM could shed light on its pathogenesis and potential therapeutic targets. Chen et al. systematically identified somatic genomic alterations in glioblastoma. GBM cells commonly exhibit amplification of chromosome 7 and deletion of chromosome 10, along with sustained activation of the PI3K/AKT/RAS, p53 and cell cycle pathways 11. The most widely accepted molecular subtypes were presented in 2010 as proneural, neural, classical and mesenchymal based on integrated genomic analysis 12. Although all subtypes share the same essential characteristics, there are some differences in features among subtypes. The proneural subtype has alterations in PDGFRA and a point mutation in IDH1. The proneural subtype might be the only subtype that benefits from first-line bevacizumab treatment 13. The mesenchymal subtype correlates with poor outcome and resistance to irradiation 14, 15. In vivo, hypoxic stress stimulates 3444glioblastoma to transform into the mesenchymal subtype to adapt to the microenvironment 16, 17. The transcriptional network for mesenchymal transformation was first reported by Maria et al. in 2011, whom identified C/EBPβ and STAT3 as initiators and master regulators 18. However, systematic descriptions of the mesenchymal subtype are still needed.

Transforming growth factor-β (TGF-β) is a pluripotency cytokine that plays a dominant role in epithelial-to-mesenchymal transition (EMT) 19. The TGF-β pathway is the putative oncogenic driver of group 3 in medulloblastoma 20. TGF-β treatment induced a mesenchymal subtype in serous ovarian cancer 21. In glioma, TGF-β drives proliferation, promotes migration and invasion, and stimulates angiogenesis 22. TGF-β receptor type 2 (TGFBR2) is the sole ligand-binding receptor of TGF-β1-3. This function gives TGFBR2 a key position in the TGF-β pathway and indicates that TGFBR2 is a potential therapeutic target that cannot be ignored.

In this study, we characterized mesenchymal glioblastoma at the mRNA-miRNA level and reported a ceRNA network that could separate the mesenchymal subtype from other subtypes. Six genes (TGFBR2, RUNX1, PPARG, ACSL1, GIT2 and RAP1B) that interacted with each other in both a ceRNA-related manner and in terms of their protein functions were identified as markers of the mesenchymal subtype. The coding sequence (CDS) and 3'-untranslated region (UTR) of TGFBR2 upregulated the expression of these genes, whereas TGFBR2 inhibition by siRNA or miR-181a/d suppressed their expression levels. Furthermore, mesenchymal subtype-related genes and the invasion phenotype could be reversed by suppressing the six mesenchymal marker genes, which indicated that these marker genes are key mRNAs in ceRNA networks (micNETs) of the mesenchymal subtype and identified the TGF-β pathway as a therapeutic target for patients with mesenchymal GBM.

Methods

Glioblastoma dataset

The glioblastoma dataset was collected from TCGA data portal (https://tcga-data.nci.nih.gov/tcga/), CGGA data portal (http://www.cgga.org.cn/), and Rembrandt dataset (https://caintegrator.nci.nih.gov/rembrandt/). The sequencing platform of the database and the number and subtype of the patients are recorded in Table .

ceRNA network construction

First, 583 GBM gene expression datasets (Agilent microarray) were downloaded from TCGA legacy database. Then, 560 out of 583 GBMs with matched whole-genome gene expression (Agilent microarray platform) and miRNA expression data were downloaded from TCGA database. A Spearman correlation test was applied to calculate the correlation coefficient and P values of the genes and miRNAs, and the P values were corrected by the Benjamini and Hochberg method. Gene-miRNA pairs with corrected P values less than or equal to 0.05 and correlation coefficients less than 0 were retained. Gene-miRNA pairs that were not included in public miRNA target databases, such as miRanda, miRTarBase, or TargetScan, were removed. Gene pairs (ceRNA pairs) with less than 10 shared miRNAs were filtered out. For each ceRNA pair, a hypergeometric test was applied using the following formula:

, where N is the total number of miRNAs, K is the number of miRNAs interacting with the specified gene, n is the number of miRNAs interacting with the candidate ceRNA of the chosen gene, and c is the number of miRNAs shared between these two genes. A Benjamini and Hochberg multiple comparison was applied to calculate the corrected P value. ceRNA pairs with corrected P values greater than 0.05 were removed.

Cell lines, clinical samples and patient-derived xenografts

The human GBM cell lines U87 and LN229 were purchased from ATCC. U87 and LN229 cells were cultured in complete DMEM containing 10% FBS, 100 units/mL penicillin, and 50 μg/mL streptomycin and incubated at 37 °C and 5% CO2. Cells at ~70% confluence were transfected with CDS, the 3'UTR or siRNA of TGFBR2, or a miR-181a/d overexpression vector, as indicated. The TGFBR2-CDS and TGFBR2-3'-UTR plasmids were synthesized by Shanghai Integrated Biotech Solutions. siRNA of TGFBR2 and mimics of miR-181a and miR-181d were ordered from GenePharma. The transfection reagent Lipofectamine® 3000 was purchased from Thermo Fisher.

The tumors were surgically resected from patients diagnosed with GBM. After resection, tumor tissues were cut into 1-mm3 blocks. A portion of the blocks was treated with DMSO or LY2109761 (Selleck) for 48 h in complete DMEM. The remaining portion of the blocks was surgically implanted into mice to establish patient-derived xenografts (PDX). All patients signed a consent form.

Real-time PCR

After the indicated treatments, GBM cells or tissue blocks were lysed in TRIzol reagent (Sigma-Aldrich). The lysate was mixed well with chloroform and centrifuged for 15 min at 13000 ×g at 4 °C. RNA in the top aqueous phase was collected and mixed with propan-2-ol to precipitate total RNA. The cDNA was synthesized with 1 µg of total RNA following the protocol of a Promega reverse transcriptional kit. The results of real-time PCR were normalized to the corresponding GAPDH mRNA levels. Relative gene expression was analyzed by the 2∆∆-Cq method. The primers were as follows:

TGFBR2-CDS: forward: 5'-GTCTGTGTGGCTGTATGGAGA-3', reverse: 5'- ACAAGTCAGGATTGCTGGTGT-3'

TGFBR2-3'UTR: forward: 5'-CTCTAGGCACCCTCCTCAGT-3', reverse: 5'- CCCTACGGTGCAAGTGGAAT-3'

RUNX1 forward: 5'-TTTCCAGTCGACTCTCAACGG-3', reverse: 5'- GCCGTAGTACAGGTGGTAGG-3'

PPARG: forward: 5'-AACTTTGGGATCAGCTCCGT-3', reverse: 5'- TGAGGGAGTTGGAAGGCTCT-3'

GIT2: forward: 5'-TTCAACGCCCATGAGTTTGC-3', reverse: 5'- GCTGTCATAGTCGGGCTGA-3'

ACSL1: forward: 5'-ACCTCAAGGTGCTTCAACCC-3', reverse: 5'- CACCATCAGCCGGACTCTTC-3'

RAP1B: forward: 5'-TGCACAACAGTGTATGCTTGAAA-3', reverse: 5'- ATGTGGACTGTGCTGTGATGG-3'

RUNX2 forward: 5'-GAACCCAGAAGGCACAGACA-3', reverse: 5'- CCTAGGCACATCGGTGATGG-3'

LIF: forward: 5'-TCTTGGCGGCAGTACACAG-3', reverse: 5'- GAGGTGCCAAGGTACACGAC-3'

SHC1: forward: 5'-CCGGGACCTGTTTGACATGA-3', reverse: 5'- TGCAAGCCAGTGAGCACATA-3'

NRP2: forward: 5'-CACCAGAACTGCGAGTGGAT-3', reverse: 5'- ATGGAGCCCGAGGAGATGAT-3'

SOCS3: forward: 5'-CAGCTCCAAGAGCGAGTACC-3', reverse: 5'- TGTCGCGGATCAGAAAGGTG-3'OSMR: forward: 5'-TACGCGTCAGAGTTTGCACT-3', reverse: 5'- CTCAGGGAACTTGGCATCGT-3'

BACE2: forward: 5'-CCTTCTTCGACTCCCTGGTG-3', reverse: 5'- CGCATCAAACACCTTCTGGG-3'PLAU: forward: 5'-AGGCTTAACTCCAACACGCA-3', reverse: TTTCCACTGTGGGTCAGCAG-3'

SERPINH1: forward: 5'-CCTGCTAGTCAACGCCATGT-3', reverse: 5'- GGATGATGAGGCTGGAGAGC-3'

COL4A1: forward: 5'-CCACCAGGACAAAAGGGTGA-3', reverse: 5'- GCATCCTGGAATACCTGGGG-3'

ACTN1: forward: 5'-GCTGCACCCTTCAACAACTG-3', reverse: 5'- GATCTCCTGAGGCGTGATGG-3'

PDPN1: forward: 5'-GAACCAGCGAAGACCGCTAT-3', reverse: 5'- CCTTCCCGACATTTTTCGCA-3'

SERPINE1: forward: 5'-GAGAACCTGGGAATGACCGA-3', reverse: 5'- CCACAAAGAGGAAGGGTCTG-3'

ITGA5: forward: 5'-CACCACCAGCAAAAACGGGA-3', reverse: 5'- TGCCTTGGTCCATTGCACA-3'

RNA-binding protein immunoprecipitation and RIP-qPCR assays

RNA immunoprecipitation (RIP) experiments were performed following the instructions of the RIP™ RNA-Binding Protein Immunoprecipitation system (Merck Millipore). Rabbit anti-AGO2 antibody (1:50, Abcam) was used to precipitate RNAs in U87 cells. Lysates were incubated with the antibody for 4 h at 4 °C. The coprecipitated RNAs were detected by real-time PCR 23.

Preparation of plasmids and determination of 3'-UTR activity

The 3'-UTR reporters of the 6 core genes were obtained from Ibsbio (Shanghai, China). The sequences of the reporters are shown as follows:

TGFBR2: GCTAGCAGTTTCTTTTGCTTATATGTTAATAGTTTTACCCTCTGCATTGGAGAGAGGAGTGCTTTACTCCAAGAAGCTTTCCTCATGGTTACCGTTCTCTCCATCATGCCAGCCTTCTCAACCTTTGCAGAAATTACTAGAGAGGATTTGAATGTGGGACACAAAGGTCCCATTTGCAGTTAGAAAATTTGTGTCCACAAGGACAAGAACAAAGTATGAGCTTTAAAACTCCATAGGAAACTTGTTAATCAACAAAGTCGAC

PPARG: GCTAGCCAGAGAGTCCTGAGCCACTGCCAACATTTCCCTTCTTCCAGTTGCACTATTCTGAGGGAAAATCTGACACCTAAGAAATTTACTGTGAAAAAGCATTTTAAAAAGAAAAGGTTTTAGAATATGATCTATTTTATGCATATTGTTTATAAAGTCGAC

RUNX1: GCTAGCTTGTTTATTTTAGACTTTAAAATGAGCTACTTCTTATTCACTTTTGTAAACAGCTAATAGCATGGTTCCAATTTTTTTTAAGTTCACTTTTTTTGTTCTAGGGGAAATGAATGTGCAAAAAAAGAAAAAGAACTGTTGGTTATTTGTGTTATTCTGGATGTATAAAAATCAATGGAAAAAAATAAACTTTCAAATTGAAATGACGGTATAACACATCTACTGAAAAAGCAACGGGAAATGTGGTCCTATTGTCGAC

GIT2: GCTAGCAAGAGTTCCTTAGAGAGTGTGCTCTAGCTGGCTGCCTGCTACCCCAAAGCTCTGAAAGAGTAAGGTTGAAGTGCTGGAACATGCTCAGATTGTTTTTTAGTTACTCAAAAGCCCCTTCTTGTTCTGAATGTGAAAGCTCTTCCCGTTCCTGCTTACTGCAGCCATTCTGTGTGTGTGGTTCAGACCCTGGCCTAGGGCCGGTGATGTGGATGAGTGAAGGGAGCCAGGTACACGGCAGCAGGGCAGAGGAGTCGAC

ACSL1: GCTAGCTATATACAAGCACAACAGGGCTTGCACTAAAGAATTGTCATTGTAATAACACTACTTGGTAGCCTAACTTCATATATGTATTCTTAATTGCACAAAAAGTCAATAATTTGTCACCTTGGGGTTTTGAATGTTTGCTTTAAGTGTTGGCTATTTCTATGTTTTATAAACCAAAACAAAATTTCCAAAAACAATGAAGGAAAGTCGAC

RAP1B: GCTAGCCTTGCTCATGCATAAGTGTATTTGCAATACCAAATATACAGGTTTAGTATTTTTGCCTGTTAGTGATTGTTTCACATGTGTAACGTTTTGGTTGAGATGTTAAATGGTGGACGAGTACTGTGGATGTGAATGTGGGAAGTAATTTTAATCATATGTAATTGGTCACAAGGCCTAATTTGCAGTAACTATTGCTGTTTTATTTAACAATGCCTTGTTGCTTTGTATGCATTAATGTTTGGATGTAAAGATTGTCGAC

The vector of the reporter plasmids was pmirGLO. To analyze 3'-UTR activity, ~95% confluent U87 cells in 96-well plates were cotransfected with the 3'-UTR reporters of the six genes and either TGFBR2-3'-UTR or si-TGFBR2. Twenty-four hours after transfection, the cells were lysed and used to determine firefly and Renilla luciferases activities using the Dual-GLO® Luciferase Assay System from Promega (E2920).

Immunohistochemical staining

Immunohistochemistry was performed on PDX tissues by subjecting 4-µm paraffin sections to a three-step process and a DAB staining kit (ZSGB-BIO). In brief, formalin-fixed, paraffin-embedded tissue sections were dried at 80 °C for 15 min and dewaxed in xylene, rinsed in graded ethanol, and rehydrated in double-distilled water. For antigen retrieval, slides were pretreated by steaming them in sodium citrate buffer (10 mM sodium citrate, pH 6.0) for 15 min at 100 °C. After the sections were washed with phosphate-buffered saline for 3 min, they were immunostained with primary antibodies against Ki-67, TGF-β and CD34 from ZSGB-BIO and incubated at 4 °C overnight. After 3 washes in PBS buffer, the tissues were covered by HRP-conjugated anti-mouse/rabbit antibodies for 30 min at room temperature (RT). Staining reactions were performed by covering the tissue samples with the prepared DAB chromogen solution and incubating them for approximately 10 min to allow proper brown color development 24.

Immunofluorescence staining

U87 cells were seeded into a confocal Petri dish until ~90% confluence was reached. Then, the cells were fixed in 4% paraformaldehyde for 30 min at RT. Immunofluorescence staining was conducted with antibodies against TGFBR2 (R&D, 1:200), RUNX1 (Abcam, 1:200), PPARG (Proteintech, 1:200), GIT2 (Abcam, 1:200), ACSL1 (Proteintech, 1:200) and RAP1B (R&D, 1:200). The cells were washed with PBS and incubated with Alexa Fluor 594 (Life Technologies, 1:200). Then, phalloidin (Life Technologies, 1:200) was employed to stain the F-actin in the cells. Nuclei were stained with DAPI, and the cells were visualized by an FV-1200 laser scanning confocal microscope.

Statistical analysis

We developed a predictive model for the mesenchymal subtype and a risk score model for patient survival based on the expression level (expr) of six micNETs and the corresponding β values from logistic regression analysis or the regression coefficients derived from univariate Cox regression analysis. The predictive model score = 2.268 × exprTGFBR2 + 0.827 × exprRUNX1 - 0.019 × exprPPARG + 0.469 × exprGIT2 + 0.279 × exprRAP1B + 0.454 × exprACSL1 - 0.835. The survival risk score = 0.001 × exprTGFBR2 + 0.057 × exprRUNX1 + 0.048 × exprPPARG + 0.238 × exprGIT2 + 0.024 × exprRAP1B + 0.041 ×exprACSL1 25.

Receiver operating characteristic (ROC) curves were constructed to determine the discriminatory capacity of micNETs or TFs for diagnosis. Kaplan-Meier analysis was performed to estimate the survival time of different subgroups, and a log-rank test was used to determine prognostic differences. Cox proportional hazards regression analyses and Pearson correlation tests were used to determine significant differences. Comparisons of binary and categorical patient characteristics between subgroups were performed using Fisher's exact test. STRING, a database of protein-protein interactions, was also used. STRING takes a list of gene products and returns a diagram of known and potential interactions 26. Gene Ontology (GO) and KEGG Pathway analyses were performed when gene sets (Pearson r > 0.4) were submitted to the DAVID website (http://david.abcc.ncifcrf.gov/home.jsp). ClueGO 27, a Cytoscape APP, was used to perform functional analysis. Biological processes were further analyzed through gene set enrichment analysis (GSEA) 28. Single-sample GSEA (ssGSEA) was used to calculate the enrichment score of every gene set for every sample 29. Heatmaps were constructed and produced by Gene Cluster 3.0 and Gene Tree View software. All statistical computations were performed with the statistical software environment R version 3.4.1 (http://www.r-project.org/), GraphPad Prism Version 6.01 or Microsoft Excel 2013. P < 0.05 was considered statistically significant.

Results

Identification of a ceRNA network enriched in the mesenchymal subtype

We collected 560 GBMs with matched whole-genome gene expression and miRNA expression data from TCGA database to establish the ceRNA network for GBM. Correlations between genes and miRNAs were calculated using Spearman correlations, and then, gene-miRNA pairs that were not included in public databases, miRanda, miRTarBase, or TargetScan were removed (Figure ). We used a stringent cutoff of 10 miRNAs shared between any two ceRNA pairs and relaxed the number of miRNAs shared between ceRNA pairs to eight and five. For the cutoff of five miRNAs shared between any two ceRNA pairs, we obtained a ceRNA network of 600397 connections and 5355 genes (Table ). For the cutoff of eight miRNAs shared between any two ceRNA pairs, we obtained a ceRNA network of 101545 connections and 2587 genes (Table ). The ceRNA networks under different numbers of shared miRNAs between two ceRNA pairs also retain key features: the ceRNA networks are characterized by two distinct, highly connected subnetworks and sparse connections between them. Interestingly, 58.5% of the genes (blue nodes) in the subnetwork were significantly enriched in the mesenchymal subtype, whereas 41.5% of the genes (yellow nodes) were mutually exclusive in the mesenchymal subtype (Figure ). Additionally, genes in the ceRNA network were not enriched in specific chromosomes but were instead uniformly distributed throughout the genome, suggesting their universal and powerful roles in tumorigenesis (Figure ). Overall, the ceRNA networks we constructed are robust with regard to the number of miRNAs shared between two ceRNA pairs.

After gene pairs with fewer than 10 shared miRNAs were filtered out, gene pairs in the ceRNA network were finally determined based on hypergeometric test results. This process led to a GBM ceRNA network of 1525 mRNAs and 188 miRNAs (Figure ), which were characterized by two distinct, highly connected subnetworks and sparse connections between them (Figure ). Interestingly, 831 genes (88% of the blue nodes) in the subnetwork were significantly enriched in the mesenchymal subtype, while 446 genes (77% of the yellow nodes) were mutually exclusive in the mesenchymal subtype (P = 1.32×10-7, Table ). GO and KEGG pathway enrichment analyses revealed that the subnetwork with genes highly expressed in the mesenchymal subtype was enriched in EMT-related pathways (Figure and Table ), whereas the subnetwork of lowly expressed genes was associated with proneural functions (Figure and Table ). Additionally, genes in the ceRNA network were not enriched in specific chromosomes but were instead uniformly distributed throughout the genome, suggesting their universal and powerful roles in tumorigenesis (Figure ).

We collected three ceRNA networks, including the Hermes GBM ceRNA network (Figure ) 5, hallmark-associated ceRNA networks in 20 major cancers (including GBM) (Figure ) 30, and cancer-associated ceRNA networks from starBase v2 (Figure ) 31. The Hermes GBM ceRNA network was established by the information-theoretic Hermes algorithm 5, which predicts ceRNA coupling between two mRNAs based on the relative size of their shared miRNA regulatory program, as predicted by the Cupid algorithm 5, 32, 33, and the conditional mutual information between each mRNA and each of their shared miRNAs 34. In the hallmark-associated cancer ceRNA network, AGO-CLIP data were used to support miRNA-target interactions during ceRNA network construction 30. The ceRNA network from starBase v2 31 was manually extracted with the default parameter settings for genes covered in our ceRNA network. We first compared the three ceRNA networks and discovered that a small percentage of connections was shared between any two ceRNA networks (6.25% of connections shared between the Hermes GBM ceRNA network and starBase v2, 2.2% between the Hermes GBM ceRNA network and the hallmark-associated cancer ceRNA network, 4.55% between the hallmark-associated cancer ceRNA network and starBase v2). We also compared our ceRNA network with the three networks and found that 6.83%, 0.5% and 6.66% of the connections were shared by our ceRNA network with the Hermes GBM ceRNA network, the hallmark-associated cancer ceRNA network, and the cancer ceRNA network in starBase v2, respectively. Considering the different methods used to establish the three cancer ceRNA networks, we established a high-confidence ceRNA network based on our ceRNA network by extracting connections reported in at least one of the three other ceRNA networks (Table ). The high-confidence ceRNA network consists of two distinct subnetworks (Figure ): the subnetwork of genes highly expressed in the mesenchymal subtype that are enriched in EMT-related pathways, and the subnetwork of lowly expressed genes that are associated with proneural functions (Figure ). Additionally, genes in the ceRNA network were not enriched in specific chromosomes but were instead uniformly distributed throughout the genome (Figure ), suggesting their universal and powerful roles in tumorigenesis.

Genes with high degrees of connectivity tended to be at the center of each ceRNA subnetwork (Figure ). A total of 26 out of 27 genes with connection degrees higher than 300 were associated with the mesenchymal subtype (Figure ). A protein-protein interaction subnetwork including 26 genes was extracted from the STRING database. Six genes (TGFBR2, RUNX1, PPARG, GIT2, ACSL1 and RAP1B) were mutually connected with each other in both the ceRNA network (Figure ) and the protein interaction network (Figure ). Pearson correlation analysis indicated that TGFBR2 might play a key role in the six micNETs, whereas RAP1B might weakly interact with the other micNETs (Figure ).

As shown in Figure , the 26 highly expressed genes in the ceRNA network were predicted to bind 51 miRNAs. Their six core genes were targeted by a total of 44 miRNAs, including all members of the miR-181 family (Figure ). According to a report by Kim et al. 35, none of these 44 miRNAs were enriched in the neuro-mesenchymal precursor cluster, which was consistent with our results on all 44 miRNAs highly expressed in the nonmesenchymal subtypes. miR-181 family genes were characteristically downregulated in the mesenchymal subtype and were negatively associated with the six core genes in GBMs (Figure ). TGFBR2 was found to have a high degree of connectivity (> 300) in the ceRNA network, and its expression levels were significantly correlated with those of the other five micNETs. TGFBR2 is also an indispensable receptor in the TGF-β pathway 36. We overexpressed the CDS and 3'-UTR of TGFBR2 in U87 and LN229 cells to investigate whether TGFBR2 could regulate the other micNETs. Real-time PCR analysis indicated that the expression of micNETs was enhanced by TGFBR2 in both a CDS and 3'-UTR transfection-dependent manner (Figure ). To further confirm the ceRNA interactions of the six micNETs, we referred to a previous work 5 and constructed 3'-UTR luciferase reporters of TGFBR2, RUNX1, PPARG, GIT2, ACSL1 and RAP1B. TGFBR2 siRNA inhibited the activity of the 3'-UTRs of these genes, while overexpression of TGFBR2 3'-UTR significantly upregulated the luciferase activity of these constructs (Figure ), suggesting the presence of ceRNA interactions among the six micNETs.

As the 3'-UTR of TGFBR2 could not enhance its own expression, upregulation of the other micNETs by TGFBR2-3'-UTR tended to affect the ceRNA network. miRNAs function as bridges in the ceRNA network. Then, we knocked down TGFBR2 by siRNA or overexpressed miR-181a/d, which were predicted to target key micNETs in the ceRNA network. The mRNA levels of the six micNETs were inhibited by downregulation of TGFBR2 or overexpression of miR-181a/d (Figure ). It is known that miRNAs and their target mRNAs can form a complex with AGO2 protein. miR-181a/d was predicted to target the six micNETs and exhibited low expression levels in the mesenchymal subtype, suggesting that the miR-181 family may participate in EMT by mediating the expression levels of TGFBR2 and the other five micNETs. We overexpressed miR-181a and miR-181d in U87 cells and examined micNET-AGO2-miR-181a/d complex binding to further investigate ceRNA regulation between TGFBR2 and other micNETs. RIP analysis indicated that miR-181a/d promoted the binding of micNETs to the AGO2 complex, which partly explained the opposite expression patterns of micNETs and miRNAs in mesenchymal subtypes and confirmed the ceRNA interactions among micNETs (Figure ).

Six micNETs, TGFBR2, RUNX1, PPARG, GIT2, ACSL1 and RAP1B, predict the GBM mesenchymal subtype

We applied the six micNETs to TCGA (1119 samples), CGGA (272 samples) and Rembrandt (227 samples) databases to investigate whether the six micNETs could predict the mesenchymal subtype. Previous studies have reported that six TFs (STAT3, C/EBPβ, RUNX1, FOSL2, BHLHE40 and ZNF238) 18 are marker genes of the mesenchymal subtype. ssGSEA enrichment analysis revealed that the six micNETs and six TFs were significantly enriched in the mesenchymal samples (Figure and Figure ), suggesting that the micNETs could act as surrogate markers of mesenchymal subtypes.

We used the risk score method to construct a model for the prediction of survival or subtype based on these six genes and investigated their prognostic and predictive value 25, 37. The signature risk score was calculated for each of the 481 patients with glioblastoma in TCGA Agilent dataset and was then applied to validate samples in other datasets.

ROC curves were generated based on the risk score model to evaluate the performance of the micNETs and TFs. The micNETs performed slightly better than the TFs in classifying samples as mesenchymal subtypes in TCGA Agilent GBM database (micNETs: AUC = 0.8883, sensitivity = 80.13%, specificity = 84%; six TFs: AUC = 0.8647, sensitivity = 83.33%, specificity = 76.92%, Figure and Figure ). The optimal cutoff for the risk score was determined according to the ROC curve, revealing that 94.3% and 87.2% of the mesenchymal samples were successfully classified as the mesenchymal subtype in TCGA HiSeq and Rembrandt databases, respectively; meanwhile, the frequencies of the six TFs were 79.2% and 80.7%, respectively, in TCGA HiSeq and Rembrandt databases (Figure and Figure ).

Survival analysis revealed that patients with a high-risk score of micNETs had poor overall survival in all databases (Figure , TCGA Agilent: P = 0.0332, TCGA u133a: P = 0.0234, TCGA RNA-seq: P = 0.0171 and CGGA RNA-seq: P = 0.0105). The distribution of risk score and the corresponding cutoff are shown in Figure . After adjustments for chemotherapy, age at diagnosis and Karnofsky performance score (KPS), Cox regression analysis revealed that the risk score of micNETs was still correlated with the overall survival of patients with GBM (Figure , TCGA Agilent: P = 4.30×10-2, HR = 1.626, TCGA u133a: P = 3.66×10-2, HR = 1.301, TCGA RNA-seq: P = 7.03×10-3, HR = 1.290 and CGGA RNA-seq: P = 3.72×10-2, HR = 2.276).

Patients with a low-risk score significantly benefited from TMZ chemotherapy in TCGA Agilent GBM dataset (Figure , P < 0.0001). In addition, we confirmed the results in TCGA RNA-seq (P = 0.0461), TCGA u133a (P = 0.0005) and CGGA RNA-seq (P < 0.0001) datasets, further demonstrating the good prognosis of groups with low-risk scores after TMZ treatment (Figure ). Although patients in the high-risk score group also benefited from TMZ treatment, their treatment effect was much worse than that of the low-risk score group in TCGA Agilent (P = 0.0222), TCGA RNA-seq (P = 0.3139), TCGA u133a (P = 0.0663) and CGGA RNA-seq (P = 0.0025) datasets (Figure ).

TCGA subtypes are associated with particular signaling pathways and distinct immune response statuses 12, 38. Thus, exploring the roles of micNETs involved in the immune response could provide new potential agents to target mediators of tumor-induced immunosuppression in GBM. We parsed TCGA and CGGA datasets to show that patients with high micNET scores had a significantly higher correlation with immunosuppressors and immune effectors than patients with low micNET scores 12, 39 (Figure ). For further comprehensive study of micNET-related immunological processes in malignant brain tumors, we downloaded the classifier gene sets for specific immune cell lineages 40, 41, which were subsequently defined as metagenes. Based on the single-sample GSEA, we observed that the micNET score was positively associated with immune cell lineages, such as B cells (r = 0.669), helper T cells (r = 0.587), cytotoxic T cells (r = 0.488), myeloid cells (r = 0.807), monocytes (r = 0.683), NK cells (r = 0.647), dendritic cells (r = 0.691) and T-cell lineages (r = 0.557) in TCGA GBM Agilent dataset (Figure ). A similar pattern of GBM in TCGA GBM RNA-seq, TCGA GBM u133a and CGGA GBM RNA-seq datasets was observed (Figure ). As revealed above, the six micNETs play a key role in immunosuppression in GBM. We overlapped micNET-related genes (r > 0.4) of the CGGA and TCGA datasets with immunosuppressors (Figure ). The six micNETs were significantly correlated with ARRB2, BTK, CASP4, CD4, CD59, ELF4, LAIR1, LYN, MR1, MSR1, NFKB1, NLRP3, PML, RAB27A, REL, SWAP70, TGFB1, TGFBR2, and TRPM2. Our results indicated that the six micNETs provide new evidence regarding immunotherapy and the infiltration of immune cells into glioma tumors.

The TGF-β pathway is a potential therapeutic target for the GBM mesenchymal subtype

In 2010, Carro et al. reported that C/EBPβ, STAT3 and four other TFs participate in the activation of mesenchymal gene expression in malignant glioma and showed that activation of this small regulatory module is necessary to initiate and maintain an aberrant phenotypic state in cancer cells 18. MicNETs with connection degrees of more than 200 were subjected to ClueGO analysis and shown to be enriched in positive regulation of EMT (P = 8.27×10-3) and the transforming growth factor beta receptor complex (P = 4.41×10-4) (Figure and Table ). Based on the COMPARTMENTS subcellular localization database 42, we predicted the subcellular locations of the six micNET proteins (Figure ). The TGFBR2 protein was primarily located in the cell membrane (evidence: 100%); most of the RAP1B, ASCL1 and GIT2 proteins were enriched in the cytosol (evidence: 100%, 99.7%, and 99.9%, respectively); and most of the RUNX1 and PPARG proteins were in the nucleus (evidence: 94.0% and 99.8%, respectively). We also used immunofluorescence to validate the subcellular localization of proteins in glioblastoma cells (Figure ).

The six TFs reported by Carro et al. can bind the promoter regions of the 41 “mesenchymal” genes. ClueGO analysis of these target genes (“mesenchymal” genes) revealed that they are involved in cell adhesion mediated by integrin (P = 7.11×10-6), positive regulation of the JAK-STAT cascade (P = 6.63×10-4), and the integrin-mediated signaling pathway (P = 8.06×10-5) (Figure and Table ). After overlapping these two gene sets, we selected 14 mesenchymal subtype-related genes that potentially drive the phenotypic transition. Overexpression of the CDS and 3'-UTR of TGFBR2 upregulated the expression levels of these mesenchymal-subtype related genes. LIF, SERPINE1 and RUNX2 are micNETs with high connectivity in the ceRNA network and are sensitive to the expression levels of the 3'-UTR of TGFBR2. Furthermore, inhibition of TGFBR2 expression levels by siRNA led to downregulation of these micNETs (Figure ), confirming that the TGF-β pathway participates in the regulation of mesenchymal genes.

We generated xenografts from tumor tissues resected from two patients with GBM to further validate the functional roles of the TGF-β pathway in mesenchymal glioblastoma (Figure ). Thirteen 1-mm3 tumor tissues were cultured in complete DMEM and treated with either DMSO or LY2109761, an inhibitor of TGFBR2. Real-time PCR analysis indicated that the mesenchymal-related gene expression pattern was downregulated by LY2109761 (Figure ). Tissue blocks of TBD0207 and TBD0220 were subjected to RNA-seq and miRNA-seq analysis. The sequencing results revealed that LY2109761 could upregulate miR-181a/b expression and downregulate the expression of mesenchymal genes, such as MET, LIF, IL-6 and OSM (Figure ). GO and KEGG pathway analyses of the differentially expressed genes (FDR < 0.001, log2FC > 1) revealed that they are involved in EMT (P = 5.77×10-6, Figure and Table ), positive regulation of mesenchymal cell proliferation (P = 5.49×10-7), the TGF-β signaling pathway (P = 4.23×10-2), cytokine-cytokine receptor interactions (P = 2.80×10-6), and the PI3K-Akt signaling pathway (P = 3.15×10-3).

Then, we employed the PDX model to test the efficiency of LY2109761 in vivo. The drug was intraperitoneally administered (Figure ), and the mice were sacrificed on the 31st day. Tumor volume was clearly reduced in the LY2109761 treatment group and was accompanied by decreased Ki-67, TGF-β and CD34 staining (Figure ). Then, the tumors were excised and subjected to RNA-seq. All groups had similar expression density distributions. Hierarchical clustering showed that the same processing groups had the closest distances (Figure ). The expression profiles of clusters determined by Mfuzz showed that the genes of mesenchymal subtype-related micNETs were downregulated by LY2109761 (Figure and Table ).

Discussion

In the present study, a ceRNA network was constructed for GBM based on TCGA GBM Agilent microarray data. mRNAs in ceRNA networks in glioblastoma were designated micNETs. In the ceRNA network, these micNETs were assembled into two subnetworks that were each characterized by hubs linked with genes with a low degree of connectivity. Recently, Yang's team performed paired-end RNA sequencing and microarray analyses to obtain the expression profiles of mRNAs, lncRNAs, circRNAs and miRNAs of three glioblastoma and matched normal samples. They found that 501 lncRNAs, 1999 mRNAs, 2038 circRNAs and 143 miRNAs were differentially expressed between glioblastoma and matched normal brain tissues 43. Cao et al. analyzed the lncRNAs regulating mRNAs as ceRNAs in glioblastoma and constructed a long noncoding RNA-mediated ceRNA network by integrating genome-wide lncRNA and mRNA expression profiles. Survival analyses demonstrated that the lncRNAs and network module are potential prognostic biomarkers for GBM.

In our work, the two subnetworks of micNETs had different gene expression levels in the mesenchymal subtype, suggesting that they might have different biological functions in shaping GBM phenotypes. In addition, the expression patterns of micNETs in the two ceRNA subnetworks are completely opposite between the proneural and mesenchymal subtypes, suggesting different biological functions in the two GBM subtypes. Most non-GCIMP mesenchymal GBMs were reported to arise from a proneural-like precursor. Loss of NF-1 could convert the proneural subtype to a mesenchymal subtype 44. It is not surprising that NF-1 is present in the proneural subtype-enriched subceRNA network. Therefore, low expression of NF-1 appears to not only promote the mesenchymal subtype by eliminating the suppression of the RAS pathway but also participate in posttranscriptional regulation.

There were more highly expressed micNETs in the mesenchymal subtype than in the proneural subtype. A total of 26 out of the 27 genes with degrees greater than 300 were in the hubs of the ceRNA subnetwork associated with the mesenchymal subtype. However, a set of 26 genes is still too large for clinical applications. We employed STRING protein-protein interaction networks to narrow down the micNET and establish core micNET signatures. TGFBR2, RUNX1, PPARG, GIT2, ACSL1 and RAP1B were identified as the core micNET signature as they interacted with each other in both the ceRNA and protein-protein networks. TGF-β was reported to stimulate the expression and function of RUNX1 45. Thus, it makes sense that the CDS of TGFBR2 upregulated the expression of RUNX1 (Figure ). The six micNETs are directly sponged by 44 miRNAs in the network, accounting for a quarter of the total miRNAs in the ceRNA network. Notably, the miR-181 family and several other miRNAs had decreased expression in the mesenchymal subtype. Many studies have been conducted on the miR-181 family in glioblastoma 46. Specifically, in our previous study, miR-181d was shown to have great prognostic value and was recommended for use as a molecular biomarker for gliomas by Chinese Glioma Cooperative Group (CGCG) clinical practice guidelines 47. As the promoter of the miR-181 family has many methylation sites, epigenetic regulation of miR-181ds might be present in the mesenchymal subtype. Tri-methylation of H3K27 was reported to be facilitated by the RTK/EGFR/AKT pathway in GBM 48, whereas EGFR overexpression plus NF-1 mutation were common in the mesenchymal subtype 12. These genomic phenotypes might trigger the silencing of miR-181s and other miRNAs, leading to a deregulated ceRNA network (Figure ).

The six micNETs were enriched in the mesenchymal subtype and could also predict the mesenchymal subtype, as verified in 1618 clinical patient samples (Figure and Figure ), suggesting their great potential for clinical diagnosis. Microarray chips or PCR kits with these micNETs may rapidly characterize glioblastoma samples. The proteins of these six micNETs have different cellular locations. As a cytokine receptor, TGFBR2 is indispensable for TGF-β1, TGF-β2 and TGF-β3 signal transduction 49. RUNX1 is a TF that regulates multiple functions and is a potential therapeutic target for leukemia 50-52. RUNX1 functions as a miRNA hub in hematopoiesis in both a transcriptional and posttranscriptional manner 53. In glioblastoma, RUNX1 is also involved in a miRNA-mRNA interaction network that regulates oncogenic pathways 5 and is considered a key TF in the mesenchymal subtype 18. PPARγ is a nuclear hormone receptor transcription factor that has a well-known role in lipid metabolism 54. Recently, PPARγ was found to regulate EMT after internal radiation (IR) of glioma, which leads to radioresistance 55. ACSL1, which initiates cellular long-chain fatty acid metabolism, was reported to enhance the proliferation of hepatoma cells 56. GIT2 is a multifunctional protein that was recently discovered to promote DNA repair induced by IR 57 and indicates a potentially poor prognosis in radiotherapy. RAP1B mediates TMZ resistance in GBM cells by cytoskeletal remodeling, and this effect could be reduced by miR-181 58. Thus, the marker micNETs selected by bioinformatics analysis reveal complex tumor-promoting functions.

As the TGF-β pathway impacts the micNET signature and major ceRNA networks, targeting the TGF-β pathway could be valuable in the mesenchymal subtype. Although TMZ is still a first-line chemotherapy for patients with GBM, our results indicated that patients with high expression levels of micNET signatures responded poorly to TMZ. TMZ induces cytotoxicity by methylating the O6 and N7 positions of guanine and causing O6-methylguanine lesions. Meanwhile, DNA damage repair confers resistance to TMZ 59, 60. Indeed, TGF-β was shown to promote an effective DNA damage response 61, 62. Blockade of the TGF-β pathway enhanced the radiation response and prolonged the survival of glioblastoma in a U87 subcutaneous model and GBM-derived CSLC 63. Inhibiting TGF-β can abolish the radiation resistance mediated by the tumor microenvironment in GBM 64. Tumors with active TGF-β signaling exhibit resistance to DNA-damaging agents by regulating DNA repair genes 65, 66. Our research confirmed the effect of LY2109761 in glioblastoma tissue blocks and GBM PDX. PDX models have the advantage of retaining key characteristics of clinical tumors 67. Indeed, the PDX models in this study possessed the same characteristics of invasiveness and rapid progression as the source tumor in the patient. LY2109761 suppressed the expression of mesenchymal signatures in tissue blocks and PDX in vivo. Genes (such as LIF, RUNX2 and SERPINE1) were inhibited in both tissue blocks, which was consistent with our measurements in GBM cell lines.

Although TMZ remains a first-line chemotherapy for patients with GBM, our results indicated that patients with high expression levels of micNET signatures responded poorly to TMZ. Therefore, a treatment combining TMZ with TGF-β pathway inhibitors is a better therapeutic strategy for patients with the mesenchymal subtype than TMZ treatment alone. Recently, immunotherapy with PD-1/PD-L1 inhibitors showed great promise for cancer treatment. Nivolumab, pembrolizumab and atezolizumab are the listed drugs for tumor therapy 68-70. PD-1/PD-L1 are highly expressed and have predictive value for the mesenchymal glioblastoma subtype 71, 72. The mesenchymal subtype of GBM tends to exhibit low purity and exhibits an intensive immune phenotype 73. Together, this evidence indicates that the mesenchymal glioblastoma subtype might benefit from anti-PD-1/PD-L1 immunotherapy. However, the TGF-β pathway tends to drive immune evasion when tumors are treated with PD-1/PD-L1 blockade 74, 75. Therefore, combined treatment targeting the TGF-β pathway and PD-1/PD-L1 immune check point will provide new possibilities for mesenchymal glioblastoma therapy.

In conclusion, our study revealed that the ceRNA signature of mesenchymal glioblastoma is enriched in the TGF-β pathway and characterized this subtype at the mRNA-miRNA level. We found that key micNET signatures could be used to diagnose mesenchymal GBM and discovered that the TGF-β pathway is a potential therapeutic target for mesenchymal glioblastoma.

Supplementary figures and table legends.

Click here for additional data file.

Supplementary table 1.

Click here for additional data file.

Supplementary table 2.

Click here for additional data file.

Supplementary table 3.

Click here for additional data file.

Supplementary table 4.

Click here for additional data file.

Supplementary table 5.

Click here for additional data file.

Supplementary table 6.

Click here for additional data file.

Supplementary table 7.

Click here for additional data file.

Supplementary table 8.

Click here for additional data file.

Supplementary table 9.

Click here for additional data file.

Supplementary table 10.

Click here for additional data file.

Supplementary table 11.

Click here for additional data file.