Analysis of microRNA expression in CD133 positive cancer stem‑like cells of human osteosarcoma cell line MG-63

Cell culture

MG-63 and Saos-2 cells were provided by the China Center for Type Culture Collection as in our previous study (Shu et al., 2020). MG-63 and Saos-2 cells were cultured in RPMI‑1640 and McCoy’s 5A medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA), respectively, supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA), 1% L‑glutamine, and 1% penicillin‑streptomycin sulfate (Thermo Fisher Scientific, Inc., Waltham, MA, USA) at 37 °C in a 5% CO₂ humidity-controlled incubator (Chen et al., 2019; Zhou et al., 2016).

Fluorescence‑activated cell sorting

MG-63 and Saos-2 cell suspensions were prepared by trypsinization, and 1 × 10⁶ cells in 500 μl were stained with phycoerythrin-labeled anti-CD133 (1:50, Miltenyi Biotec GmbH, Bergisch Gladbach, Germany) at 4 °C for 60 min. After washing, the CD133(+) and CD133(−) MG-63 cells were sorted by fluorescence-activated cell sorting (FACS) using a BD FACS AriaIII system (BD Biosciences).

RNA extraction and miRNA array

Total RNA was extracted from MG-63 cells using TRIzol reagent (Takara, Japan) and the miRNeasy mini kit (Qiagen, Wset Sussex, United Kingdom) following the manufacturer’s instructions. The extracted RNAs were labeled based on the miRCURY™ Hy3™/Hy™ Power labeling kit (Exiqon, Vedbaek, Denmark) and hybridized on a miRCURY™ LNA Array (version 18.0, Exiqon, Vedbaek, Denmark). After washing, the slides were scanned using an Axon GenePix 4000B microarray scanner (Axon Instruments, Foster City, CA, USA). The scanned images were imported into the GenePix Pro 6.0 platform (Axon Instruments) for grid alignments and analysis. Replicated miRNAs were averaged, and miRNAs with expression intensities of ≥50 in all of the samples were used to calculate the normalized expression using the median normalization method. The DEmiRNAs were identified by volcano plot filtering. Finally, hierarchical clustering was performed to identify DEmiRNAs using MEV software (Version 4.6; TIGR, Microarray Software Suite 4, Boston, MA, USA).

Prediction of target genes of DEmiRNAs and collection of validated target genes

DEmiRNA–mRNA interactions were predicted using TargetScan (http://www.targetscan.org/), miRwalk (http://www.ma.uni-heidelberg.de/apps/zmf/mirwalk/), and the miRDB (http://www.mirdb.org/) database. The validated DEmiRNA–mRNA interactions were collected from miRTarBase (https://mirtarbase.cuhk.edu.cn/).

GO and KEGG pathway enrichment analyses of target genes of DEmiRNAs

Target genes of DEmiRNAs were obtained, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and Gene Ontology (GO) analysis were performed using DAVID (Guo et al., 2019). The cut-off criterion for both analyses was P < 0.05.

Construction of miRNA–mRNA pathway network

The association among the differently expressed miRNAs, mRNAs and the mRNA miRNAs, mRNAs and pathways in the network were represented by nodes of different shapes and colors.

qRT-PCR

Total RNA was extracted using TRIzol reagent (Takara, Japan). Complementary DNA (cDNA) was then synthesized through reverse transcription of the RNA using the Prime-Script RT regent Kit and gDNA Eraser (TaKaRa) (Wang et al., 2020a). Subsequently, quantitative real-time PCR (qRT-PCR) was performed using the SYBR Premix Ex Taq II Kit (Takara, Kusatsu, Japan) and an ABI 7500 qRT-PCR system (Applied Biosystems, Waltham, MA, USA), using the DNA as template. U6 was used as the internal control. The relative miRNA expression levels were calculated based on the 2^–ΔΔCt equation. The primers (Sangon Biotech, Shanghai, China) are shown in Table 1.

Sphere formation assay

MG‑63 and Saos‑2 cells (500 cells/well) were plated in Ultra‑Low Attachment 24‑well plates (Corning, Inc., Corning, NY, USA) with 0.8% methyl cellulose (Sigma‑Aldrich, St. Louis, MO, USA; Merck KGaA, Darmstadt, Germany) supplemented with 20 μl/ml B27, 20 ng/ml bFGF, 10 ng/ml epidermal growth factor, 1% L‑glutamine and 1% penicillin‑streptomycin sulfate (all were obtained from Thermo Fisher Scientific, Inc., Waltham, MA, USA). Every 3 days, each well was examined under a light microscope (IX71; Olympus Corporation, Tokyo, Japan).

Invasion assay

To assess the invasive ability, 2 × 10⁵ serum-starved cells were seeded in 200 μl medium without serum and then plated in the top of a Transwell™ chamber (24-well insert; pore size, 8 μm; Corning, Corning, NY, USA) that was coated with diluted Matrigel (BD Biosciences). After 8 h, the infiltrating cells were stained with 4,6-diamidino-2-phenylindole (DAPI) and counted using a microscope. The invasive ability of cells was quantified.

MiR-4284 knockdown

FACS‑sorted CD133(+) and CD133(−) cells were seeded into 6-well plates and incubated overnight. They were transfected with siRNA against miR-4284 or negative control (Ribobio Co) at 60 nM using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instructions. After 72 h, the cells were harvested and RNA and protein were isolated.

Western blot

Proteins were isolated using a Total Protein Extraction Kit (Thermo Fisher Scientific) according to the manufacturer’s protocol. The following primary antibodies were used to detect protein expression: anti-β-catenin (#8480; Cell Signaling), anti-phospho-β-catenin (#9561; Cell Signaling), and anti-β-actin (#4970; Cell Signaling). Antibodies were diluted as specified in the specifications.

Statistical analysis

Data were analyzed using GraphPad Prism 8 (GraphPad Software, Inc., San Diego, CA, USA). Continuous data are expressed as mean ± SD. Differences between groups were analyzed using the two-tailed Student unpaired t-test. Results were considered statistically significant if P < 0.05.

DEmiRNAs between CD133(+) and CD133(−) MG-63 cells

FACS-sorted CD133(+) and CD133(−) MG-63 cells (Fig. 1A) were grown under serum-free conditions for the miRNA array. DEmiRNAs between CD133(+) and CD133(−) MG-63 cells were identified. In total, 20 DEmiRNAs (fold change ≥ 1.2, P < 0.05) were found (16 upregulated in CD133(+) cells and 4 downregulated). Hierarchical clustering graphs and volcano plots of the 20 DEmiRNAs between CD133(+) and CD133(−) MG-63 cells are shown in Figs. 1B and 1C. The DEmiRNAs are listed in Table 2.

Functional enrichment analysis of DEmiRNA target genes

In total, 9,687 predicted target genes of the 20 DEmiRNAs were obtained. Of these predicted target genes, experimental evidence reported in the literature was obtained for 2,099 validated target genes.

To obtain a deeper understanding of the roles of the 20 DEmiRNAs, GO enrichment analysis and KEGG pathway analysis were performed on predicted and validated target genes. GO analysis indicated that among the predicted target genes, several biological processes, including “regulation of cell morphogenesis” and “cell junction organization” (P < 0.05), were significantly enriched (Fig. 2A). KEGG pathway enrichment analysis indicated that among the predicted target genes, several KEGG pathways were significantly enriched (P < 0.05), including “Wnt signaling pathway,” “MAPK signaling pathway,” “PI3K-Akt signaling pathway,” and “Ras signaling pathway” (Fig. 2B).

GO analysis indicated that among the validated target genes, several biological processes, including “chromatin assembly” and “regulation of protein stability” (P < 0.05), were significantly enriched (Fig. 3A). KEGG pathway enrichment analysis indicated that among the validated target genes, several pathways were statistically enriched (P < 0.05), including “Cell cycle,” “Fc epsilon RI signaling pathway,” “ErbB signaling pathway,” and “Rap1 signaling pathway” (Fig. 3B).

Validation and functional exploration of the top three upregulated miRNAs and top 3 downregulated miRNAs

The expression of the top three upregulated miRNAs and top three downregulated miRNAs was further validated via RT-PCR. The results confirmed that the expression levels of hsa-miR-4284, hsa-miR-4485-3p, and hsa-miR-3656 were significantly increased in CD133(+) cells compared with CD133(−) cells (Figs. 4A–4C), while the expression levels of hsa-miR-487b-3p, hsa-miR-493-5p, and hsa-miR-431-5p were significantly decreased in CD133(+) cells compared with CD133(−) cells (Figs. 4B–4F).

After validation by RT-PCR, GO analysis and KEGG pathway enrichment analysis of the predicted and validated target genes of the top three upregulated miRNAs and top three downregulated miRNAs were also performed (Figs. 5 and 6). GO analysis indicated that among the predicted target genes, several biological, including biological processes such as “Wnt signaling pathway,” “calcium modulating pathway,” “transcription from RNA polymerase II promoter,” and “positive regulation of GTPase activity” (P < 0.05), were significantly enriched (Fig. 5A). KEGG pathway enrichment analysis indicated that among the predicted target genes, several KEGG pathways were significantly enriched (P < 0.05), including “Wnt signaling pathway,” “MAPK signaling pathway,” “Rap1 signaling pathway,” “PI3K-Akt signaling pathway,” “cGMP-PKG signaling pathway,” and “Ras signaling pathway” (Fig. 5B). GO analysis indicated that among the validated target genes, several biological processes, including “cell cycle,” “stem cell population maintenance,” and “mesenchymal cell differentiation” (P < 0.05), were significantly enriched (Fig. 6A). KEGG pathway enrichment analysis indicated that among the validated target genes, several pathways were statistically enriched (P < 0.05), including “p53 signaling pathway” and “PPAR signaling pathway” (Fig. 6B). The results of the functional enrichment analysis of these six confirmed that DEmiRNA target genes are in agreement with the functional enrichment analysis of the total number of DEmiRNA target genes.

We observed a significant interconnection between the molecular roles and signaling pathways of predicted and validated target mRNAs (Tables S1 and S2). For instance, in the interconnection among predicted target mRNAs and multiple signaling pathways, WNT7A and FZD5 regulates “Wnt signaling pathway,” “Proteoglycans in cancer,” and “Pathways in cancer”. IGF1R regulates “Proteoglycans in cancer,” “PI3K-Akt signaling pathway,” “Pathways in cancer,” “Ras signaling pathway,” “Rap1 signaling pathway,” and “Transcriptional misregulation in cancer”. AKT3 regulates “PI3K-Akt signaling pathway,” “Ras signaling pathway,” “Pathways in cancer,” “Rap1 signaling pathway,” “cGMP-PKG signaling pathway,” “Proteoglycans in cancer,” and “MAPK signaling pathway”. IGF1R is the predicted target gene for miR-493-5p and miR-431-5p. AKT3 is the predicted target gene of hsa-miR-493-5p. In the interconnection among validated target mRNAs and multiple signaling pathways, CDKN1A and CDK4 regulate the p53 signaling pathway. WNT7A and FZD5 are predicted target genes of hsa-miR-4284. CDKN1A is a validated target gene for hsa-miR-493-5p. CDK4 is a validated target gene for hsa-miR-431-5p and hsa-miR-4284.

MiR-4284 regulates the self-renewal capacity and invasion ability in CD133(+) cells through Wnt signaling

To evaluate the effect of miR-4284 on the cells’ self-renewal ability, we initially performed a sphere formation assay of CD133(+) and CD133(−) cells transfected with siRNA targeting miR-4284. MiR-4284 knockdown significantly reduced the sphere forming capacity of CD133(+) from MG-63 and Saos-2 cells (Fig. 7A). Next, we investigated the role of miR-4284 in the regulation of the invasion ability of OSCs. After 48 h, the number of invading cells in each group was calculated. The results showed that the number of invasive CD133(+) MG-63 and Saos-2 cells was markedly decreased after miR-4284 knockdown (Fig. 7B).

Next, we investigated the relationship between miR-4284 and β-catenin phosphorylation. As shown in Fig. 7C, the level of phospho-β-catenin was upregulated, while the expression level of total β-catenin remained unaltered. Collectively, these data demonstrate that miR-4284 may promotes the self-renewal ability and invasion of OCSCs through Wnt signaling.