This article is mentioned in:

Introduction

FK506 (also known as tacrolimus) is a potent immunosuppressive agent that is used to prevent graft rejection and to treat autoimmune diseases, both experimentally and clinically (1). Over the past decade, several studies have reported that FK506 treatment may increase bone mass and may be a potential treatment for bone defects (2,3). FK506 exposure has been reported to induce osteogenic differentiation by promoting the expression of osteoblastic transcription factors, including core-binding factor α1 (also termed runt related transcription factor 2, Runx2), osteopontin (OPN) and osteocalcin (OCN) (2,3). A previous in vivo study using bone tissue/hydroxyapatite composite demonstrated increased bone formation and higher osteogenic parameters when cultured with FK506 (4). However, other studies reported that FK506 may result in osteoporosis by increasing the expression levels of osteoclast differentiation factor or by inhibiting osteoblast differentiation by suppressing the activity of calcineurin and the expression of Runx2 (5,6). These results suggested that the function of FK506 during osteogenic differentiation was contradictory and required for further study.

MicroRNAs (miRNAs) are small, 22-nucleotides-long, endogenous noncoding RNAs that bind to target mRNAs at the 3′-untranslated region to mediate translation and reduce protein expression levels (7). They are key negative regulators of diverse biological and pathological processes, such as proliferation, migration, apoptosis and differentiation (7,8). A number of miRNAs have been demonstrated to be involved in the differentiation of osteoblasts from bone marrow stromal cells (BMSCs) by regulating specific gene expression (9). One previous study has indicated that miRNA (miR)-96 may promote osteogenic differentiation by suppressing the expression of heparin-binding epidermal growth factor-like growth factor in mouse BMSCs (10). During the process of osteogenic differentiation of BMSCs, miR-205 was downregulated and inhibition of miR-205 enhanced osteogenic differentiation by regulating the expression of special AT-rich sequence-binding 2 (SATB2) and Runx2 (11). However, identification of specific miRNAs and their regulatory roles in the FK506-induced osteodifferentiation have been poorly investigated.

The aim of the present study was to examine the effects of various concentrations of FK506 (5–5,000 nM) on the osteogenic differentiation of rat BMSCs. Differentially expressed miRNAs were profiled by miRNA array, verified by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) and subjected to gene ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

Materials and methods

Isolation and culture of rat BMSCs

Primary BMSCs were harvested from 5 male Sprague Dawley rats (age, 4–5 weeks; weight, 80–100 g) using previously described methods (2). The male Sprague-Dawley rats were purchased from Shanghai SLAC Laboratory Animal Co., Ltd. (Shanghai, China). The rats were maintained in a temperature controlled room (22±3°C) with a 12-h light/dark cycle and were allowed free access to drinking water and food. All procedures were approved by The Animal Research Committee of Zhongshan Hospital, Fudan University (Shanghai, China). Briefly, rats were sacrificed by cervical dislocation followed by soaking in 75% ethanol for 10 min. The tibias and femurs from both legs of each rat were dissected and the two ends of each bone were cut under aseptic conditions. The marrow was flushed out with rat Mesenchymal Stem Cell Growth medium (Cyagen Biosciences, Inc., Santa Clara, CA, USA) containing 10% fetal bovine serum, 1% penicillin/streptomycin and 1% glutamine. The cells were centrifuged at 500 × g for 5 min at 4°C, resuspended with Mesenchymal Stem Cell Growth medium, plated in 25 cm2 culture flasks and incubated at 37°C in 5% CO2. The medium was replaced every 3 days and non-adherent cells were removed. The cells reached ~80% confluence at 10–14 days of culture, and were subsequently dissociated using TrypLE Express (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and subcultured to passage 3 for further analysis. The cells were analyzed using an Olympus SZX12 (Olympus Corporation, Tokyo, Japan) inverted microscope equipped with a digital camera and connected to a PC using MagnaFire 2.0 camera software (Optronics, Goleta, CA, USA).

Phenotype analysis

BMSC cell-surface markers were analyzed by flow cytometry. Briefly, cells at passage 3 with 80% confluence were trypsinized using TrypLE Express (Gibco; Thermo Fisher Scientific, Inc.), washed twice with PBS and centrifuged at 500 × g for 5 min at 4°C. Cells were diluted with stain buffer (BD Bioscience, Franklin Lakes, NJ, USA) and the cell number was determined by cell counting assay following the manufacturer's instructions (Precision Instrument Co. Ltd., Shanghai, China). Following counting, a 100 µl cell suspension containing 5×105 cells was added to flow tubes, mixed with rat antibodies, including 1 µl FITC-labeled anti-CD29 (1:100; cat. no. 102205; BioLegend, Inc., San Diego, CA, USA), anti-CD45 (1:100; cat. no. 202205; BioLegend, Inc.) and anti-CD90 (1:100; cat. no. 206105; BioLegend, Inc.) and 2.5 µl phycoerythrin (PE) -labeled anti-CD34 (1:40; cat. no. ab187284; Abcam, Cambridge, UK), and incubated in dark for 40 min at 4°C. Following two washes with PBS, the stained cells were resuspended in 300 µl PBS and immediately analyzed using a BD FACSCalibur flow cytometer (BD Biosciences, Franklin Lakes, NJ, USA). Identification of rat BMSCs was performed in triplicate. The mouse immunoglobulin G1, κ was used as an isotype control (BioLegend, Inc.) and stain buffer was used as blank control. Data were analyzed using FlowJo software version 10 (Tree Star, Inc., Ashland, OR, USA).

Osteogenic differentiation

Cells at passage 3 with 80% confluence were trypsinized using TrypLE Express and seeded into six plates at 1×104 cells/cm2 and incubated with Mesenchymal Stem Cell Growth medium at 37°C for 24 h. Following incubation, the growth medium was replaced with 2 ml osteo-induction medium was prepared as previously described (12) containing 10 mM Na β-glycerophosphate (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany), 0.25 mM l-ascorbic acid (MP Biomedicals, LLC, Santa Ana, CA, USA) and 0, 5, 50, 500 or 5,000 nM FK506 (Sigma-Aldrich; Merck KGaA). The osteo-induction medium was replaced every other day; cells were collected at 3, 7 and 14 days for RT-qPCR analysis, and phenotype staining was done at 7 or 14 days. All biochemical assays were carried out at least three times.

Alkaline phosphatase (ALP) and alizarin red S (ARS) staining

ALP staining was performed following 7 days osteogenic stimulation. Briefly, cells were rinsed twice with PBS and fixed with 4% formaldehyde for 30 min at room temperature. The fixed cells were stained with 5-bromo-4-chloro-3-indolyl phosphate/nitroblue tetrazolium solution (Beyotime Institute of Biotechnology, Shanghai, China) for 30 min at room temperature and washed 3 times with PBS. ARS staining was used to assess mineral deposition produced at late-stage bone formation. Briefly, following 14 days osteogenic stimulation, cells were washed twice with PBS, fixed with 4% formaldehyde for 30 min at room temperature, washed twice with PBS and then stained with the Alizarin Red S Solution (Cyagen Biosciences, Inc., Chicago, USA) for 30 min at room temperature. Following ARS staining, the cells were washed 3 times with PBS. Images were captured using an Olympus SZX12 inverted microscope equipped with a digital camera and connected to a PC using MagnaFire 2.0 camera software.

RNA extraction and RT-qPCR analysis

According to the manufacturer's protocol, total RNA was extracted from 1×105 control and FK506-induced BMSCs with TRIzol Reagent (Invitrogen, Carlsbad, USA), purified using the RNeasy Mini Kit (Qiagen, CA, USA) and quantitated by using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). RNA (500 ng) was reverse transcribed into cDNA using the PrimeScript™ RT-PCR kit (Takara Bio, Inc., Otsu, Japan). The cDNA was amplified using the SYBR Premix Ex Taq II kit (Takara Bio, Inc.) with the GeneAmp-PCR system 7500 (Thermo Fisher Scientific, Inc.). Total RNA for miRNA analysis was isolated using a miRCute miRNA Isolation kit (Tiangen Biotech Co., Ltd., Beijing, China) and the expression of mature miRNA was determined by miRCute miRNA PCR Detection kit (Tiangen Biotech Co., Ltd.) with the GeneAmp-PCR system 7500. All primers used for RT-qPCR are listed in the Table I, and each sample was measured in triplicate. Gene expression results of mRNA or miRNA were evaluated by the 2−ΔΔCq method (13) and normalized to GAPDH or U6, respectively.

Table I. Primer sequences for reverse transcription-quantitative polymerase chain reaction analysis.

Table I.

Primer sequences for reverse transcription-quantitative polymerase chain reaction analysis.

Gene	Primer sequence (5′→3′)
Sp7	F: GAGGCACAAAGAAGCCATACA
	R: GGGAAAGGGTGGGTAGTCAT
OCN	F: ACAAGTCCCACACAGCAACTC
	R: CCAGGTCAGAGAGGCAGAAT
Runx2	F: CACCTCTGACTTCTGCCTCTG
	R: GATGAAATGCCTGGGAACTG
OPN	F: CTTGGCTTACGGACTGAGG
	R: GCAACTGGGATGACCTTGAT
Smad7	F: TCGGAAGTCAAGAGGCTGTG
	R: CTGGACAGTCTGCAGTTGGTT
MAPK9	F: AACTCGCTACTATCGGGCTC
	R: TGGGAACAGGACTTTATGGAGG
Smad5	F: CCAGTGTTAGTGCCTCGTCA
	R: GTGGAAGGAATCAGGAAACG
Bambi	F: CATTGCTGGCGGACTGAT
	R: TCCCTTCTTGGAGTGGTGTC
Dusp2	F: GGTGGTCCTGTGGAAATCTT
	R: GAATGCTCTTGTAGCGGAAAA
Jag1	F: TGCTTGGTGACAGCCTTCTA
	R: TGGGGTTTTTGATTTGGTTC
GAPDH	F: CAGTGCCAGCCTCGTCTCAT
	R: AGGGGCCATCCACAGTCTTC

[i] Bambi, BMP and activin membrane-bound inhibitor; Dusp2, dual-specificity phosphatase 2; F, forward; Jag1, jagged 1; MAPK9, mitogen-activated protein kinase 9; OCN, osteocalcin; OPN, osteopontin; R, reverse; Runx2, runt-related 2.

miRNA microarray analysis

miRNA microarray analysis was performed by KangChen Bio-tech Inc. (Shanghai, China). Briefly, total RNA was extracted from 1×107 control and FK506-induced BMSCs using TRIzol reagent and purified with RNeasy Mini kit (Qiagen, Inc., Valencia, CA, USA), according to manufacturer's protocol. The quality and quantity of RNA were measured by using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and RNA integrity was determined using 1.5% gel electrophoresis. miRNAs were labeled with Hy3/Hy5 fluorescence dyes using the miRCURY LNA microRNA Array Power Labeling kit (Exiqon A/S, Vedbaek, Denmark) and hybridized in a Hybridization Chamber II (Ambion; Thermo Fisher Scientific, Inc.) using the miRCURY LNA microarray kit (Exiqon A/S) according to the manufacturer's protocol. Following hybridization, the microarray slides were processed with an Axon GenePix 4000B Microarray Scanner (Axon Instruments; Molecular Devices, LLC, Sunnyvale, CA, USA), and data analysis was performed using GenePix Pro software version 6.0 (Molecular Devices, LLC). The experiment was performed with 3 samples both in control and treatment group, and miRNAs were defined as differentially expressed if changes in expression levels were ≥2-fold than those of the controls.

GO analysis and KEGG pathway annotation based on miRNA expression profile

Target genes of the differentially expressed miRNAs were predicted using miRBase (http://www.mirbase.org), miRanda (http://www.microrna.org), TargetScan (http://www.targetscan.org) and miRDB (http://mirdb.org/miRDB) databases. GO analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov) to determine the functions of the predicted target genes and to uncover the miRNA-target gene regulatory network based on the predicted biological processes and molecular functions. KEGG pathway analysis (http://www.genome.jp/kegg) was used to determine the potential pathways that the predicted target genes may be a part of. Fisher's exact test and the χ2 test were used to classify the GO category and KEGG pathway, and the false discovery rate (FDR) (14) was calculated to correct the P-value. Thresholds of P<0.05 and FDR<0.05 were used to select significant GO categories. Fisher P-value stands for the enrichment P-value used Fisher's exact test and Fisher-P value<0.05 was used as a threshold to select significant KEGG pathways.

Statistical analysis

Statistical analysis of two groups was determined by unpaired Student's t-test. All data were represented as the mean ± standard error of the mean. P<0.05 was considered to indicate a statistically significant difference.

Results

Growth characteristics and phenotype of BMSCs using flow cytometry

Following 3 days of primary culture, the isolated BMSCs were mostly spindle shaped, single or small colony size and grew at a low density. (Fig. 1A). Primary BMSCs cultured for 12 days replicated rapidly and reached ~80% confluence (Fig. 1B). Cells were expanded by successive subculture and exhibited colony-like growth at passage 3. BMSCs phenotypes were verified by flow cytometric analysis. The subcultured BMSCs at passage 3 were positive for MSCs markers CD29 (99.9%) and CD90 (99.3%), and negative for hematopoietic lineage markers CD34 (0.998%) and CD45 (0.064%) (Fig. 1C).

Figure 1. Growth characteristics and phenotypic analysis of rat BMSCs. Inverted microscope images (40x) reveal that BMSCs in primary culture were (A) mostly spindle-shaped and single cells at day 3 and (B) reached almost 80% confluence at day 12. Flow cytometric analysis indicated that BMSCs at passage 3 were positive for CD29 and CD90, and negative for CD34 and CD45 (C). BMSCs, bone marrow stromal cells; FITC, fluorescein isothiocyanate; PE, phycoerythrin.

FK506-induced osteogenic differentiation of BMSCs

To detect the effects of FK506 treatment on the osteogenic differentiation of BMSCs, cells were collected at passage 3 and cultured in osteo-induction medium with a concentration of FK506 ranging between 5 and 5,000 nM. ALP and ARS staining were strongest in cells treated with 50 nM FK506 (Fig. 2A). RT-qPCR analysis of BMSCs treated with 50 nM FK506 revealed that the expression of Sp7, Runx2, and OPN were significantly increased following 7 days stimulation and OCN expressed significantly at 7 and 14 days. (P<0.05 vs. negative control; Fig. 2B). These results indicated that FK506 treatment has the potential to stimulate the BMSCs differentiation into osteoblasts at a concentration of 50 nM for 7 days of stimulation, and thus this concentration was used for further miRNA array experiments.

Figure 2. ALP and ARS staining in BMSCs after FK506-induced osteogenic differentiation. (A) ALP and ARS staining were strongest in rat BMSCs when stimulated with FK506 at the concentration of 50 nM. (B) Expression levels of osteoblastic specific genes were significantly increased at 50 nM FK506 treatment for 7 days. *P<0.05 and **P<0.01 compared with NC. Data are expressed as the mean ± standard error of the mean; n=3 per group. ALP, alkaline phosphatase; ARS, alizarin red S; d, days; NC, negative control; OCN, osteocalcin; OPN, osteopontin; Runx2, runt-related 2.

miRNA expression analysis following FK506-induced osteogenic differentiation

To examine the differential expression of miRNAs that are involved in FK506-induced osteodifferentiation of rat BMSCs, miRCURY LNA microarrays were used to profile miRNA expression 7 days following FK506 (50 nM) induction. Comparison with the control group identified a total of 56 miRNAs, including 17 that were upregulated and 39 that were downregulated (Tables II and III, respectively). RT-qPCR analysis was used to verify the microarray results and demonstrated that five miRNAs were upregulated (miR-106b-5p, miR-101b-3p, miR-193a-3p, miR-485-3p and miR-142-3p) and four were downregulated (miR-27a-3p, miR-207, miR-218a-2-3p and let-7a-5p) compared with untreated control cells (Fig. 3).

Figure 3. Reverse transcription-quantitative polymerase chain reaction validation of microRNA expression between untreated and FK506-induced rat bone marrow stromal cells. The five upregulated and four downregulated microRNAs were confirmed to be consistent with the results identified by microarray analysis. *P<0.05 compared with NC. Data are expressed as the mean ± standard error of the mean; n=3 per group. miR, microRNA; NC, negative control.

Table II. Upregulated microRNAs in FK506-treated rat bone marrow stromal cells vs. untreated cells.

Table II.

Upregulated microRNAs in FK506-treated rat bone marrow stromal cells vs. untreated cells.

microRNA	Fold change
rno-miR-101b-3p	2.14
rno-miR-193a-3p	2.15
rno-miR-485-3p	3.03
rno-miR-33-5p	3.03
rno-miR-324-5p	3.17
rno-miR-142-3p	3.18
rno-miR-106b-5p	2.45
rno-miR-185-3p	2.67
rno-miR-339-5p	2.34
rno-miR-1843a-5p	2.32
rno-miR-1839-3p	2.22
rno-let-7i-3p	2.19
rno-miR-3068-5p	2.19
rno-miR-455-5p	2.06
rno-miR-296-3p	2.08
rno-miR-345-3p	2.42
rno-miR-3084b-5p	3.18

[i] miR, microRNA; rno, Rattus norvegicus.

Table III. Downregulated microRNAs in FK506-treated rat bone marrow stromal cells vs. untreated cells.

Table III.

Downregulated microRNAs in FK506-treated rat bone marrow stromal cells vs. untreated cells.

microRNA	Fold change
rno-miR-27a-3p	2.70
rno-miR-218a-2-3p	3.70
rno-miR-665	5.00
rno-miR-542-5p	10.0
rno-miR-207	2.63
rno-miR-329-3p	3.45
rno-miR-10a-5p	3.03
rno-miR-450a-5p	4.12
no-miR-672-3p	3.45
rno-miR-878	3.03
rno-miR-2985	2.56
rno-miR-370-5p	2.38
rno-miR-488-3p	3.23
rno-miR-193a-5p	6.67
rno-miR-29b-3p	2.17
rno-miR-465-5p	2.33
rno-miR-3560	2.13
rno-miR-3596d	3.70
rno-miR-1956-5p	50.0
rno-miR-653-3p	4.76
rno-miR-203a-3p	7.14
rno-miR-652-5p	6.25
rno-miR-200c-5p	4.17
rno-miR-181b-5p	2.13
rno-miR-1-3p	2.04
rno-miR-143-5p	2.70
rno-miR-664-2-5p	2.38
rno-let-7a-5p	2.22
rno-miR-3557-5p	2.08
rno-miR-331-5p	3.85
rno-let-7f-5p	3.57
rno-miR-18a-3p	12.5
rno-miR-30c-1-3p	2.33
rno-miR-10a-3p	25.0
rno-miR-487b-3p	2.38
rno-miR-497-5p	11.1
rno-miR-582-5p	2.56
rno-miR-493-5p	3.45
rno-miR-194-3p	5.00

[i] miR, microRNA; rno, Rattus norvegicus.

GO category and KEGG pathway analysis

To investigate the specific roles of the nine identified differentially expressed miRNAs in the regulation of osteodifferentiation in BMSCs, target genes were predicted and subjected to GO term and KEGG pathway analysis. GO analysis demonstrated that the highly enriched GO terms targeted by these dysregulated miRNAs were involved in multiple biological processes such as cell differentiation, cell growth and cell migration (Tables IV and V). KEGG pathway analysis revealed the target genes to be highly enriched in 22 upregulated and 21 downregulated signaling pathways, including mitogen-activated protein kinase (MAPK) signaling pathway, calcium signaling pathway and Toll-like receptor signaling pathway (Fig. 4).

Figure 4. KEGG pathway analysis for predicted target genes. (A) Pathways that may be upregulated in FK506-induced rat BMSCs. (B) Pathways that may be downregulated in FK506-induced rat BMSCs. The horizontal axis is the enrichment of pathways and the vertical axis is the pathway category. Fisher-P-value<0.05 was used as a threshold to select significant KEGG pathways. LgP is the negative logarithm of the Fisher-P-value. ABC, ATP-binding cassette; BMSCs, bone marrow stromal cells; GABA, γ-aminobutyric acid; MAPK, mitogen-activated protein kinase; TNF, tumor necrosis factor.

Table IV. Predicted GO term enrichment of the upregulated microRNAs in FK506-induced osteodifferentiation in rat bone marrow stromal cells.

Table IV.

Predicted GO term enrichment of the upregulated microRNAs in FK506-induced osteodifferentiation in rat bone marrow stromal cells.

GO.ID	Category name	Fold enrichment	P-value	FDR
GO:0051234	Establishment of localization	1.60	1.93×10−10	5.02×10−7
GO:0048522	Positive regulation of cellular process	1.48	3.49×10−7	3.13×10−4
GO:0048518	Positive regulation of biological process	1.42	1.01×10−6	7.80×10−4
GO:0009605	Response to external stimulus	1.71	1.34×10−6	8.68×10−4
GO:0006464	Cellular protein modification process	1.55	1.61×10−6	8.68×10−4
GO:0048731	System development	1.42	2.65×10−6	1.30×10−3
GO:2000278	Regulation of DNA biosynthetic process	8.45	3.52×10−6	1.46×10−3
GO:0007275	Multicellular organismal development	1.38	6.08×10−6	2.04×10−3
GO:0006873	Cellular ion homeostasis	2.40	2.69×10−5	5.16×10−3
GO:0048514	Blood vessel morphogenesis	2.38	3.13×10−5	5.81×10−3
GO:0016043	Cellular component organization	1.33	6.33×10−5	9.08×10−3
GO:0051247	Positive regulation of protein metabolic process	1.85	8.64×10−5	1.08×10−2
GO:0006875	Cellular metal ion homeostasis	2.36	9.88×10−5	1.19×10−2
GO:0008283	Cell proliferation	1.60	1.02×10−4	1.19×10−2
GO:0007267	Cell-cell signaling	1.86	1.26×10−4	1.33×10−2
GO:0034220	Ion transmembrane transport	1.88	1.50×10−4	1.39×10−2
GO:0051352	Negative regulation of ligase activity	1.34	1.60 10−4	1.39×10−2
GO:0012501	Programmed cell death	1.60	1.76×10−4	1.48×10−2
GO:0007167	Enzyme linked receptor protein signaling pathway	1.95	2.04×10−4	1.66×10−2
GO:0009653	Anatomical structure morphogenesis	1.46	2.43×10−4	1.74×10−2
GO:1902531	Regulation of intracellular signal transduction	1.61	3.03×10−4	1.98×10−2
GO:0051338	Regulation of transferase activity	1.94	3.54×10−4	2.08×10−2
GO:0065009	Regulation of molecular function	1.45	3.69×10−4	2.09×10−2
GO:0051270	Regulation of cellular component movement	1.93	4.61×10−4	2.32×10−2
GO:0006915	Apoptotic process	1.55	5.82×10−4	2.57×10−2
GO:0030334	Regulation of cell migration	1.95	9.55×10−4	3.71×10−2
GO:0010647	Positive regulation of cell communication	1.62	1.02×10−3	3.86×10−2
GO:0051046	Regulation of secretion	1.90	1.12×10−3	4.21×10−2
GO:0006874	Cellular calcium ion homeostasis	2.24	1.29×10−3	4.66×10−2
GO:0051928	Positive regulation of calcium ion transport	3.73	1.35×10−3	4.76×10−2

[i] FDR, false discovery rate; GO, gene ontology.

Table V. Predicted GO term enrichment of the downregulated microRNAs in FK506-induce osteodifferentiation in rat bone marrow stromal cells.

Table V.

Predicted GO term enrichment of the downregulated microRNAs in FK506-induce osteodifferentiation in rat bone marrow stromal cells.

GO.ID	Category name	Fold enrichment	P-value	FDR
GO:0044710	Single-organism metabolic process	1.55	4.56×10−10	8.19×10−7
GO:0071702	Organic substance transport	1.98	5.32×10−9	5.73×10−6
GO:0010648	Negative regulation of cell communication	2.32	4.75×10−7	1.47×10−4
GO:0042981	Regulation of apoptotic process	2.00	3.06×10−6	5.15×10−4
GO:0044237	Cellular metabolic process	1.21	2.25×10−5	2.02×10−3
GO:0048583	Regulation of response to stimulus	1.52	3.12×10−5	2.58×10−3
GO:0006464	Cellular protein modification process	1.53	3.37×10−5	2.63×10−3
GO:0030182	Neuron differentiation	1.92	3.82×10−5	2.94×10−3
GO:0050767	Regulation of neurogenesis	2.31	5.02×10−5	3.57×10−3
GO:0007399	Nervous system development	1.65	5.04×10−5	3.57×10−3
GO:0050804	Regulation of synaptic transmission	3.04	5.12×10−5	3.58×10−3
GO:0065009	Regulation of molecular function	1.59	5.81×10−5	3.89×10−3
GO:0016049	Cell growth	2.62	9.00×10−5	5.25×10−3
GO:2000026	Regulation of multicellular organismal development	1.73	9.04×10−5	5.25×10−3
GO:0044249	Cellular biosynthetic process	1.34	9.74×10−5	5.35×10−3
GO:0030799	Regulation of cyclic nucleotide metabolic process	4.18	1.44×10−4	7.08×10−3
GO:0006954	Inflammatory response	2.38	1.69×10−4	7.71×10−3
GO:0031100	Organ regeneration	4.52	1.72×10−4	7.71×10−3
GO:0043269	Regulation of ion transport	2.25	1.97×10−4	8.51×10−3
GO:0042127	Regulation of cell proliferation	1.72	2.86×10−4	9.77×10−3
GO:0045595	Regulation of cell differentiation	1.68	4.15×10−4	1.31×10−2
GO:0023014	Signal transduction by phosphorylation	2.07	4.87×10−4	1.47×10−2
GO:0055012	Ventricular cardiac muscle cell differentiation	10.00	5.82×10−4	1.67×10−2
GO:0051592	Response to calcium ion	3.77	6.68×10−4	1.84×10−2
GO:0000165	MAPK cascade	1.99	1.43×10−3	3.16×10−2
GO:0017156	Calcium ion-dependent exocytosis	4.78	1.56×10−3	3.31×10−2
GO:0030099	Myeloid cell differentiation	2.33	2.16×10−3	4.00×10−2
GO:0035051	Cardiocyte differentiation	3.35	2.78×10−3	4.78×10−2
GO:0002520	Immune system development	1.76	2.81×10−3	4.78×10−2
GO:0045622	Regulation of T-helper cell differentiation	10.00	2.98×10−3	4.85×10−2

[i] FDR, false discovery rate; GO, gene ontology; MAPK, mitogen-activated protein kinase.

Verification of the target genes of differentially expressed miRNAs

Osteodifferentiation is driven by a series of complex events that result in the activation of specific transcription factors. Of the predicted target genes and differentially regulated function analysis, only those closely related to osteoblast differentiation, ossification and signal transduction are listed in Table VI. RT-qPCR demonstrated that MAPK9, Smad5 and jagged 1 (Jag1), which are positive regulators of osteogenesis, were significantly upregulated in FK506-induced rat BMSCs undergoing osteogenic differentiation, compared with the untreated control (Fig. 5); whereas dual-specificity phosphatase 2 (Dusp2), Smad7 and BMP and activin membrane-bound inhibitor (Bambi), which are negative regulators of osteodifferentiation, were demonstrated to be significantly downregulated. The results suggested that FK506 treatment may induce osteogenic differentiation by regulating the expression of specific target genes of these differentially expressed miRNAs.

Figure 5. Six osteogenic genes that were predicted to be target genes of the differentially expressed microRNAs were confirmed by reverse transcription-quantitative polymerase chain reaction to be significantly regulated in FK506-induced rat bone marrow stromal cells. *P<0.05 and **P<0.01 compared with NC. Data are expressed as the mean ± standard error of the mean; n=3 per group. Bambi, BMP and activin membrane-bound inhibitor; Dusp2, dual-specificity phosphatase 2; Jag1, jagged 1; MAPK9, mitogen-activated protein kinase 9; NC, negative control.

Table VI. Predicted target genes and functions of microRNAs associated with FK506-induced osteodifferentiation in rat bone marrow stromal cells.

Table VI.

Predicted target genes and functions of microRNAs associated with FK506-induced osteodifferentiation in rat bone marrow stromal cells.

microRNA	Target genes	Function
rno-miR-101b-3p	Smad5	Osteoblast differentiation, ossification
	TGFβR1	Receptor for growth factor, protein metabolism
rno-miR-106b-5p	Bambi	Stem cell differentiation, TGF-β receptor signaling
	Dusp2	Tissue development, signal transduction
	MAP3K2, MAP3K8 and MAPK9	MAPK cascade, proliferation and differentiation
	Smad7	Regulator of BMP/TGF-β signaling, osteogenesis
rno-miR-142-3p	TGFβR1	Receptor for growth factor, protein metabolism
rno-miR-27a-3p	MAPK9	MAPK cascade, proliferation and differentiation
	Smad5	Osteoblast differentiation, ossification
	CBFβ	Osteoblast differentiation, osteogenesis
rno-miR-218a-2-3p	Wnt5a	Stem cell differentiation and proliferation
	Jag1	Mesenchymal cell differentiation, osteogenesis
rno-let-7a-5p	TGFβR1	Receptor for growth factor, protein metabolism
	AcvR1C	Receptor for growth and differentiation factor
	MAPK6 and MAPK9	MAPK cascade, proliferation and differentiation

[i] AcvR1c, activin A receptor type 1C; Bambi, BMP and activin membrane-bound inhibitor; BMP, bone morphogenetic protein; CBFβ, core-binding factor β; Dusp2, dual-specificity phosphatase 2; Jag1, jagged 1; MAPK, mitogen-activated protein kinase; MAP3K, mitogen-activated protein kinase kinase kinase; miR, microRNA; rno, Rattus norvegicus; TGF, transforming growth factor; TGFβR1, TGFβ receptor 1.

Discussion

FK506 is a powerful and clinically useful immunosuppressant, which has been widely used as a treatment in organ transplantation, atopic dermatitis and rheumatoid arthritis (15). During the treatment of rheumatoid arthritis, a number of studies demonstrated that FK506 treatment could increase bone formation, inhibit osteoclastic bone resorption and improve clinical symptoms (16,17). In addition, it has been reported that FK506 exposure may promote osteogenic and chondrogenic differentiation (18,19), and a change of immunosuppressive monotherapy from cyclosporin A to FK506 may be able to reverse bone loss (20), indicating that FK506 may be a potential treatment for bone defects.

However, other studies demonstrated that FK506 treatment may result in osteoporosis, prevent osteoblast differentiation and/or cause significant bone loss (5,6,21). These conflicting results may be due to different drug concentrations administered or induction time, which is crucial in the process of bone metabolism. Low-dose and short-term application of FK506 may be able to promote the early stage of osteogenesis. By contrast, long-term administration of FK506 led to accelerated bone formation and bone resorption (22). The present study demonstrated that FK506 treatment in rat BMSCs had the most significant osteogenic potential when administered at a dose of 50 nM for 7 days; whereas high concentrations of FK506 or longer periods of stimulation reduced the osteogenic potential, suggesting that the possible side effects of drug cytotoxicity may inhibit the osteoinduction of BMSCs.

miRNAs serve a key role in regulating osteodifferentiation (9). The present study identified and verified five upregulated and four downregulated miRNAs that were differentially expressed in FK506-induced rat BMSCs undergoing osteogenic differentiation, some of which have been reported as positive and/or negative regulators of osteogenesis. For example, miR-27a, an early negative regulator of osteogenesis, was revealed to delay osteoblast differentiation by suppressing SATB2 expression (23). By contrast, a different study reported that miR-27a may be able to promote osteoblast differentiation by repressing the expression of secreted frizzled-related protein at the transcriptional level (24). A previous study revealed that miR-193a was downregulated during osteogenic differentiation of human adipose-derived stem cells and the osteoporotic fractures; however, the mechanism has not been elucidated (25). In addition, miR-218 was reported to be upregulated during osteodifferentiation and its overexpression could significantly increase human adipose tissue-derived stem cells osteogenic differentiation (26). However, results from the present study are not fully in line with the results of previous studies, which may be explained by the variability of mesenchymal stem cells from different tissues and culture microenvironments.

To further investigate the molecular mechanisms of FK506-induced osteogenic differentiation in rat BMSCs, functional analysis were performed. Predicted miRNA target genes, such as MAPK9 and Dusp2, were verified by RT-qPCR. Highly enriched KEGG pathway analysis indicated that MAPK9 and Dusp2 may serve a role through the MAPK signaling pathway, which is involved in a number of cellular functions, such as differentiation, proliferation and cell death (27). Of these, the physiological functions of MAPK9 [(also known as c-Jun N-terminal kinase 2 (JNK2)] in osteodifferentiation and bone formation have been investigated previously (28). MAPK9 activity is essential for the expression of activating transcription factor 4, and its overexpression may enhance mineral deposition that is essential in the late stages osteogenic differentiation (29). Activation of JNKs is also essential for bone morphogenetic protein (BMP) 9-induced osteodifferentiation of mesenchymal stem cells (30). The expression of Dusp2, an inhibitor of osteoblast differentiation, was revealed to be downregulated during osteogenic differentiation (31,32). The DUSP subfamily (DUSP1, DUSP2, DUSP4 and DUSP5) of nucleus-inducible MAPK phosphatases are rapidly induced by growth factors or stress signals, such as mitogens, oxidative stress, heat shock and hypoxia and localize exclusively to the nucleus (32,33). In the present study, MAPK9 expression was revealed to be significantly upregulated in rat BMSCs stimulated with 50 nM FK506 for 7 days, whereas the expression of Dusp2 was downregulated. These data suggested that MAPK9 and Dusp2 are important regulatory components in the process of FK506-induced osteogenesis in rat BMSCs.

Verification of the predicted target genes also demonstrated that Smad5 expression was upregulated, and the expression levels of Smad7 and Bambi were downregulated. Smad5, Smad7 and Bambi belong to the transforming growth factor (TGF) -β/BMP superfamily, which are involved in numerous cellular processes and bone formation (34). Previous studies have indicated that FK506 may be able to induce osteogenic and chondrogenic differentiation by regulating Smads (18,19), the downstream effectors of TGF-β/Smad signaling pathway. Smad5 is a positive regulator in BMP/Smad signaling, and was previously demonstrated to be highly expressed in genistein-induced osteogenic differentiation of human BMSCs (35). Smad5 protein degradation in the Smurf1-overexpression C2C12 mouse myoblast cells was demonstrated to cause the inability to undergo BMP-induced osteoblast differentiation and Smad1/5-knockout mice exhibited severe chondrodysplasia (36,37). Smad7 and Bambi act as negative regulators of osteogenic differentiation: Smad7 was reported to inhibit proliferation, differentiation and mineralization of mouse osteoblastic cells (38); and Bambi, a target gene of miR-20a, was demonstrated to inhibit the osteogenesis of human BMSCs (39). In the present study, although the target genes of differentially expressed miRNAs were not highly enriched in the TGF-β/Smad signaling pathway, it may be concluded that these differentially expressed miRNAs may promote FK506-induced osteogenic differentiation by regulating the expression of their target genes, such as Smad5, Smad7 and Bambi.

In conclusion, the results of the present study demonstrated that FK506 treatment was able to significantly promote osteogenic differentiation of rat BMSCs at the concentration of 50 nM. Following 7 day stimulation, several miRNAs were identified as differentially expressed and may regulate osteodifferentiation by targeting genes that are involved in bone formation. These preliminary results may provide an experimental basis for further research on the functions of miRNAs in FK506-induced osteogenic differentiation.

Acknowledgements

The present study was supported by a grant from the Major Scientific and Technological Research Projects of The Science and Technology Commission of Shanghai Municipality (grant no. 14JC1490600) and a grant from The Natural Science Foundation of Shanghai (grant no. 15ZR1406200).

References