Identification of miRNA signature associated with BMP2 and chemosensitivity of TMZ in glioblastoma stem-like cells

Abstract

Glioblastoma multiform (GBM) is the most lethal intracranial tumor in adults. Glioblastoma stem-like cells (GSCs) are responsible for tumorigenesis and chemotherapy resistance. BMPs are known to increase temozolomide (TMZ) response in GSCs, however, the intracellular molecular mechanism remains largely unknown. In this study, we built a GSC cell model called U87S, and performed RNA sequencing to identify differentially expressed (DE) miRNA profiles in U87S cells treated with BMP2, TMZ or combined BMP2 and TMZ respectively. Bioinformatics analysis revealed that most DE miRNAs were involved in the cancer pathways, suggesting their crucial roles in gliomagenesis. Eight miRNAs from RNA-seq were validated. Four out of these miRNAs (has-miR-199a-3p, hsa-miR-374b-5p, hsa-miR-320d, and hsa-miR-339-5p) were found significantly up-regulated in GBM tumor tissues. One of them, hsa-miR-199a-3p, was significantly correlated with the survival of GBM patients, and differentially expressed in U87S cells. Expression of hsa-miR-199a-3p was up-regulated by BMP. Overexpression of hsa-miR-199a-3p in U87S cells inhibited cell viability and enhanced the cytotoxicity of TMZ. And activation of BMP boosted the effect of hsa-miR-199a-3p on cell viability and TMZ-mediated cytotoxicity. Besides, expressions of five predicted targets of hsa-miR-199a-3p were evaluated. Four of them were differentially expressed in GBM tumors. And one of them, SLC22A18, was associated with the survival of GBM patients. In the end, a hsa-miR-199a-3p-mediated ceRNA network was constructed for the convenience of future study. Together, our data provided DE miRNA expression profiles associated with BMP2 and TMZ in GSCs, which might lead to finding out miRNA-based target therapies that specially target GSCs.

Keywords: BMP2; Cancer stem cell; CeRNA; Gliobalstoma; RNA-seq; TMZ.

Figure 1

Characterization of U87S GSCs derived from U87MG. (A) Morphologies of U87MG vs. U87S cells. Most of the adherent growth U87MG cells are showing astrocyte-like morphology, flat cell body with short processes. The U87S cells are cultured in the serum-free neurosphere culture medium. Spheres are at the size of approximately 50–100 μM, showing positive immunofluorescence staining of stemness marker CD133 (Prominin-1) (green). (B) Flow cytometry analysis shows the percentage of CD133+ positive cells in U87S spheres was up to 89.34%, compared to 1.28% in U87MG. (C) U87S cells possess a multipotential differentiation capacity. U87S cells were dissociated and seeded into the differentiation medium containing 2% FBS for 10 days. Representative immunofluorescence staining for the astroglial marker, GFAP, and the neuronal marker, β tubulin-Ⅲ (TUBB3) respectively. Scale bar in A is 100 μM, scale bare in C is 50 μM.

Figure 2

BMP2 sensitizes U87S cells to temozolomide. (A) Cytotoxic effect of TMZ on U87S cells. U87S cells were treated with TMZ at assigned concentrations (0–4000 μM) for 24 h. Cell viability was measured by CCK-8 assay. IC50 of TMZ is 725 μM, 479 μM, and 275 μM in U87S, U87MG, and HA. The data are given as the means ± SD of three replicates. (B) BMP2 sensitizes U87S cells to TMZ. U87S cells were treated with a series of TMZ (0–1000 μM) alone (black curve), together with 20 ng/ml recombinant BMP2 (red curve) or 200 nM LDN 193189. Cell viability was measured by CCK-8 assay. The data are given as the means ± SD of three replicates. (C) BMP signaling pathway is active in U87S cells. mRNA expression levels of intracellular BMP signaling effectors, SMAD1 and SMAD5, and transcriptional targets of BMP signaling, SMAD6, and ID1 by qRT-PCR. Data is normalized vs. endogenous control: GAPDH. Results are shown in relative expression and data are means ± SD from three independent experiments. *p < 0.05, and **p < 0.01.

Figure 3

Differentially expressed miRNAs by RNA-seq. (A) Heatmap of hierarchical clustering analysis of differentially expressed miRNAs in three treated groups (U87S-BMP, U87S-TMZ, and U87S-BMP-TMZ), (fold change > 2, compared to untreated U87S respectively, FDR (false-discovery rate) < 0.001). Red color represents up-regulated miRNAs and blue color represents down-regulated miRNAs. Expression data are represented as log2 fold change versus U87S for each group. (B) Venn diagram reports the differentially expressed miRNAs in three treated groups, U87S-BMP, U87S-TMZ, and U87S-BMP-TMZ. (C-E) Functional enrichment analysis of putative target genes of DE miRNAs by GO. The ar chart represents the classification of GO Biological Processes, Cellular Component or Molecular Function. Bars represent the number of genes in the specified category, organized by p-value. (F–H) The top 20 enriched signaling pathways for putative target genes of DE miRNAs by KEGG. Enriched pathways are selected based on the enrichment factor value with the number of genes in a pathway ≥4. The blue color from dark to light represents the Q value. The size of the dots represents the gene number in each pathway.

Figure 4

Validation of 8 differentially expressed miRNAs and GO and KEGG analysis. (A) qRT-PCR was performed to assess relative miRNA expression of hsa-miR-199a-3p, hsa-miR-33a-5p, hsa-miR-409-3p, hsa-miR-374b-5p, hsa-miR-27a-5p, hsa-miR-376b-5p, and has-miR-320d. Data are presented as mean ± SD from three independent experiments. *p < 0.05 compared with untreated group, and **p < 0.01 compared with the untreated group. (B) Functional enrichment analysis of putative 7379 targets genes of 8 validated miRNAs by GO. The bar chart represents the classification of GO Biological Processes, Cellular Component or Molecular Function. Bars represent the number of genes in the specified category, organized by p-value. (C) The top 20 enriched signaling pathways for putative target genes of 8 validated miRNAs by KEGG. Enriched pathways are selected based on the enrichment factor value with the number of genes in a pathway ≥4. The blue color from dark to light represents the Q value. The size of the dots represents the gene number in each pathway.

Figure 5

Clinic relevance of 8 validated miRNAs. (A) Heatmap of 8 validated miRNAs expression patterns in 7 normal brain tissues and 16 GBM patient tissues from GEO data. The color from blue to red represents the gene count number. has-miR-320d, has-miR-339-5p, has-374b-5p and has-miR-199a-3p were up-regulated in GBM tumors (fold change ≥ 2, p < 0.05). (B) Boxplot of the hsa-199a-3pand has-miR-409-3p expression distribution in four subtypes of GBM of 196 TCGA patients. hsa-miR-199a-3p shows significantly elevated expression in M subtype of GBM patients (p = 0.001). has-miR-409-3p does not show a significant difference among subtypes of GBM. (C) Heatmap of hsa-miR-199a-3p and hsa-miR-409-3p expression in four subtypes of GBM patients. hsa-miR-199a-3p and hsa-miR-409-3p have distinct expression patterns in GBM subtypes. (D) The relevance of has-miR199a-3p expression to the survival time of GBM patients in four subtypes by Kaplan-Meier analysis. Higher expression of hsa-miR-199a-3p in C subtype (p = 0.02) and M subtype (p = 0.029) is associated with better survival outcome. (E) The relevance of hsa-miR-409-3p expression to the survival time of GBM patients in four subtypes by Kaplan-Meier analysis. The expression level of hsa-miR-409-3p in each GBM subtype is not significantly associated with the survival time of patients. (F) Survival analysis with expressions of hsa-199a-3p and hsa-miR-409-3p both considered. P subtype patients with lower than the median level of hsa-miR-199a-3p and hsa-miR-409-3p have a longer survival time than those with the upper level than the median (p = 0.005).

Figure 6

hsa-miR-199a-3p inhibits the proliferation and enhances TMZ chemosensitivity in U87S cells. (A) Detection of endogenous hsa-miR-199a-3p levels in HA, U87 MG, and U87S cells by RT-PCR. Data are presented as mean ± SD from three independent experiments. *p < 0.05 compared with HA. (B) The effect of hsa-miR-199a-3p on the proliferation and TMZ chemosensitivity in U87S. Cell viability was measured by CCK-8 assay in the U87S cells transfected with the hsa-miR-199a-3p mimics, hsa-miR-199a-3p inhibitor, or their corresponding negative control (miR-NC and anti-NC) with or without TMZ treatment respectively. Data are presented as mean ± SD from three independent experiments. *p < 0.05 compared to the NC group. (C) BMP activation boosts the inhibitory effect of hsa-miR-199a-3p on cell proliferation. Cell viability was measured by CCK-8 assay in the U87S cells transfected with the hsa-miR-199a-3p mimics, hsa-miR-199a-3p inhibitor or their corresponding negative control (miR-NC and anti-NC) with the endogenous BMP signaling (WT, black bars), activation of BMP signaling (BMP2, red bars), and inhibition of BMP signaling (LDN 193189, blue bars) respectively. Data are presented as mean ± SD from three independent experiments. **p < 0.001, *p < 0.05 compared to NC group. (D) hsa-miR-199a-3p and BMP sensitize TMZ chemosensitivity in an additive way. The U87S cells were transfected with the hsa-miR-199a-3p mimics, hsa-miR-199a-3p inhibitor or their corresponding negative control (miR-NC and anti-NC) with the endogenous BMP signaling (WT, black bars), activation of BMP signaling (BMP2, red bars), and inhibition of BMP signaling (LDN 193189, blue bars) respectively. 500 μM TMZ was used to treat cells for 24 h. Cell viability was measured with CCK-8 assay. Data are presented as mean ± SD from three independent experiments. *p < 0.05 compared to NC group.

Figure 7

Putative targets of hsa-miR-199a-3p. (A) Functional enrichment analysis of putative 1000 target genes of hsa-199a-3p by GO. The bar chart represents the classification of GO Biological Processes, Cellular Component or Molecular Function. Bars represent the number of genes in the specified category, organized by p-value. (B) The top 20 enriched signaling pathways for putative target genes of hsa-199a-3p by KEGG. Enriched pathways are selected based on the enrichment factor value with the number of genes in a pathway ≥4. The blue color from dark to light represents the Q value. The size of the dots represents the gene number in each pathway. (C) Heatmap shows expression patterns of has-miR-199a-3p and 5 of its predictive targets in RNA-seq and qRT-PCR (fold change > 2 in the treated group compared to untreated U87S, p < 0.005). A negative correlation of hsa-miR-199a-3p and these five putative targets are shown both in RNA-seq and qRT-PCR data. Data are presented as mean ± SD (n = 3). Red color represents up-regulated and blue color represents down-regulated. (D) Heatmap of expressions of AGAP3, TXNIP, SLC22A18, ANGPTL2 and CITED2 in 156 TCGA GBM patient samples compared to normal brain tissue. (E) Boxplot of the expression pattern of AGAP3. (F) Boxplot of the expression pattern of TXNIP. (G) Boxplot of the expression pattern of ANGPTL2. (H) Boxplot of the expression pattern of CITED2. (I) Boxplot of the expression pattern of SLC22A18. (J) Kaplan-Meier plot shows the association of SLC22A18 expression with patient survival. Lower expression of SLC22A18 than medium leads to longer survival time for patients.

Figure 8

The hsa-199a-3p-mediated ceRNA network. The hsa-miR-199a-3p mediated ceRNA network contained 189 mRNAs and 269 lnRNAs, which share common hsa-miR-199a-3p binding elements, can work as ceRNAs competing binding to hsa-miR-199a-3p. 5 putative targets of hsa-miR-199a-3p (AGAP3, TXNIP, SLC22A18, ANGPTL2 and CITED2) validated by qRT-PCR were marked by blue circles.