Bioinformatics Analysis Identified miR-584-5p and Key miRNA-mRNA Networks Involved in the Osteogenic Differentiation of Human Periodontal Ligament Stem Cells

Wang, Chengze; Dong, Lingling; Wang, Ying; Jiang, Zhiwei; Zhang, Jing; Yang, Guoli

doi:10.3389/fgene.2021.750827

ORIGINAL RESEARCH article

Front. Genet., 27 September 2021

Sec. Stem Cell Research

Volume 12 - 2021 | https://doi.org/10.3389/fgene.2021.750827

Bioinformatics Analysis Identified miR-584-5p and Key miRNA-mRNA Networks Involved in the Osteogenic Differentiation of Human Periodontal Ligament Stem Cells

$\r\nChengze Wang&#x;$ Chengze Wang^†

Lingling Dong^†

Ying Wang

Zhiwei Jiang

Jing Zhang

Guoli Yang^*

Key Laboratory of Oral Biomedical Research of Zhejiang Province, The Affiliated Hospital of Stomatology, School of Stomatology, Zhejiang University School of Medicine, Hangzhou, China

Human periodontal ligament cells (PDLCs) play an important role in periodontal tissue stabilization and function. In the process of osteogenic differentiation of PDLSCs, the regulation of molecular signal pathways are complicated. In this study, the sequencing results of three datasets on GEO were used to comprehensively analyze the miRNA-mRNA network during the osteogenic differentiation of PDLSCs. Using the GSE99958 and GSE159507, a total of 114 common differentially expressed genes (DEGs) were identified, including 62 up-regulated genes and 52 down-regulated genes. GO enrichment analysis was performed. The up-regulated 10 hub genes and down-regulated 10 hub genes were screened out by protein-protein interaction network (PPI) analysis and STRING in Cytoscape. Similarly, differentially expressed miRNAs (DEMs) were selected by limma package from GSE159508. Then, using the miRwalk website, we further selected 11 miRNAs from 16 DEMs that may have a negative regulatory relationship with hub genes. In vitro RT-PCR verification revealed that nine DEMs and 18 hub genes showed the same trend as the RNA-seq results during the osteogenic differentiation of PDLSCs. Finally, using miR-584-5p inhibitor and mimics, it was found that miR-584-5p negatively regulates the osteogenic differentiation of PDLSCs in vitro. In summary, the present results found several potential osteogenic-related genes and identified candidate miRNA-mRNA networks for the further study of osteogenic differentiation of PDLSCs.

Introduction

The high incidence of periodontal disease is an important cause of tooth loss. Periodontal disease can cause varying degrees of periodontal bone defects. To date, various conventional therapies for periodontal regeneration have shown limited clinical outcomes. In clinical work, the main goal of periodontal surgery is to remove infections and restore defective tissues. Periodontal regenerative therapy with membranes and bone grafting materials called guided bone regeneration (GBR) and guided tissue regeneration (GTR) has been employed with distinct levels of clinical success. Surgery can help remove the source of infections but it is often difficult to recover periodontal ligament and alveolar bone. Literature studies (Kämmerer et al., 2017; Hu et al., 2018; Liu et al., 2019b) have reported that injection therapy based on PDLSCs promoted the regeneration of alveolar bone and periodontal ligament. Repairing alveolar bone defects and periodontal ligament requires the participation of PDLSCs. Understanding the signal pathway regulation in the process of osteogenic differentiation of PDLSCs can provide a theoretical basis for periodontal tissue regeneration.

Human PDLSCs were first isolated and reported by Seo et al. (2004). Recent studies have shown that PDLSCs may offer a more reliable strategy for the treatment of periodontal defects by a cell-based tissue engineering approach (Ma et al., 2017; Zhang et al., 2019; Aprilianti et al., 2020). This treatment method relies on the multiple functions of PDLSCs such as anti-inflammatory properties (Nagata et al., 2017; Liu et al., 2019a; Qiu et al., 2020), osteogenic differentiation ability, and chemokines secretion capacity (Lee J. S. et al., 2017; Wang et al., 2018). The osteogenic differentiation of PDLSCs is the most critical among them. Many studies have shown that PDLSCs also have an impressive effect on the regeneration of bone tissue. For example, PDLSCs are used for alveolar bone regeneration (Lee J. et al., 2017), which can restore bone defect by mixing hydroxyapatite (Aprilianti et al., 2020), tricalcium phosphate (Xu X. Y. et al., 2019) and hydrogels (Ma et al., 2017). In addition, the PDLSCs membrane sheet technology can promote the osseointegration of alveolar bone and implant (Iwasaki et al., 2019). In particular, one study (Zhang et al., 2018) has shown that human PDLSCs have similar osteogenic differentiation ability compared with bone marrow-derived MSCs (BMSCs) and adipose-derived stem cells (ADSCs).

In the past, it was reported that the osteogenic differentiation of PDLSCs was mainly related to the Wnt (Gu et al., 2017), integrin (Tang et al., 2017), and PI3K-AKT (Xu X. et al., 2019) signaling pathways. However, the specific molecular regulation mechanism of osteogenic differentiation of PDLSCs is still not clear. Recently, it has been reported that the miRNA-mRNA network also plays an important regulatory role in the osteogenic differentiation process of PDLSCs (Whitfield et al., 2004; Ahmad et al., 2021). Therefore, the analysis of multiple RNA-seq and miRNA-seq datasets is conducive to detailed screening of the changes in signal pathways during the osteogenic differentiation of PDLSCs.

Microarray techniques and bioinformatics analysis have been widely used to screen for the DEGs, functional pathways and PPI associated with osteogenic differentiation of stem cells (Yang et al., 2019; Fan et al., 2020; Bini et al., 2021; Shin et al., 2021). In this study, three microarray datasets (total 16 samples) were downloaded from Gene Expression Omnibus (GEO) for analysis to identify DEGs and DEMs between the control group and the induction group. In summary, 18 hub genes and nine miRNAs were verified consistently by qPCR in vitro. Among the nine miRNAs, we predict that miR584-5p can bind to ALPL and may regulate the expression of ALPL. ALPL is a very important osteogenic regulatory gene. Liu reported (Liu et al., 2020) that miR-584-5p was also down-regulated during osteogenic differentiation of dental pulp stem cells, which was in line with our sequencing analysis results. However, the article did not do further functional verification experiments of miR-584-5p. So we further verified the role of miR-584-5p in PDLSCs Osteogenic differentiation.

Materials and Methods

Microarray Data

Gene Expression Omnibus¹ is a public functional genomics data repository. Before the study, we searched for related datasets on GEO’s website. On GEO’s website, we conducted a search using keywords “periodontal ligament stem cell” or “PDLSC” or “dental stem cell” or “periodontal ligament tissue” and “osteogenic differentiation.” We selected datasets including human PDLSCs samples with osteogenic differentiation and samples under normal medium as control. We selected GSE99958, GSE159507, and GSE159508 for further analysis. This study used three datasets, two mRNA-seq (GSE99958, GSE159507) and one miRNA-seq (GSE159508). GSE99958, GSE159507, GSE159508 were downloaded from GEO. GSE99958 contains four mRNA-seq samples, which are cultured in control medium GSM2666465, in osteogenic medium for 3 days GSM2666463, 7 days GSM2666464, and 14 days GSM2666462. GSE159507 contains 6 samples, Group Control (GSM4831419, GSM4831420, and GSM4831421) are cultured in control medium, Group Induced (GSM4831416, GSM4831417, and GSM4831418) are cultured in osteogenic medium for 14 days. GSE159508 contains six samples, three samples are under osteogenic induction for 0 days (GSM4831425, GSM4831426, and GSM4831427), and three samples are under osteogenic induction for 14 days (GSM4831422, GSM4831423, and GSM4831424). Database GSE99958 was performed on GPL17303, database GSE159507 was performed on GPL16956 Agilent-045997 Arraystar human lncRNA microarray V3 (Probe Name Version) and database GSE159508 was performed on miRCURY LNA microRNA Array, 7th generation; lot 35,106–35,106 (miRBase 21.0). Based on the platform annotation information, probes were transformed into corresponding gene symbols in the R software environment.

Differentially Expressed Analysis

For GSE99958, GSE159507, GSE159508, the probe sets without corresponding gene symbols or the genes with multiple probe sets were removed or averaged, respectively. For the dataset GSE99958, without biological replicates, CORNAS (default 1.5-fold count change) was used to obtain DEGs (Hu et al., 2018). The DEGs of 3, 7, and 14 days compared to 0 days were obtained and the common DEGs were determined. For the dataset GSE159507, the limma package was used to identify DEGs. LogFC > 1 or LogFC < −1 and P-value < 0.05 were considered to indicate statistical significance. At last, the common up-regulated genes and down-regulated genes between GSE99958 and GSE159507 were selected. For the dataset GSE159508, the limma package was also used to identify the differentially expressed miRNAs (DEMs). Up-regulated miRNA (LogFC > 1 and P-value < 0.05) and down-regulated miRNA (LogFC < −1 and P-value < 0.05) were selected.

Gene Ontology and KEGG Enrichment Analyses of Differentially Expressed Genes

We performed gene ontology (GO) analysis of common DEGs (62 up-regulated genes and 52 down-regulated genes) between GSE99958 and GSE159507. GO and KEGG enrichment analyses were performed via the package clusterProfiler on R platform. GO enrichment analysis predicted the functional roles of DEGs based on three aspects, including biological processes (BP), cellular components (CC), and molecular functions (MF). The threshold was P < 0.05.

Protein-Protein Interaction Network Network Analysis and Analysis of Hub Genes

The protein-protein pairs of DEGs (62 up-regulated genes and 52 down-regulated genes) were identified via STRING (Szklarczyk et al., 2015) v11.0². Then, the PPI networks were visualized in Cytoscape version 3.8.2 software³. CytoHubba, a plugin in Cytoscape, is used to screen hub genes, and the Degree algorithm is to select genes with the top 10 nodes ranked by degree.

Targets Prediction for Differentially Expressed miRNAs and Candidate miRNA Selection

Potential targets for DEMs were predicted by bioinformatics algorithms in the miRWalk database (Sticht et al., 2018): miRWalk⁴. Because there are many target genes from miRNA prediction and the results of the algorithm prediction have large false positives. A negative regulatory relationship between hub genes and miRNAs is used to narrow the size of candidate miRNAs from DEMs.

Isolation and Characteristics of Human Periodontal Ligament Stem Cells

The study protocol was approved by the Medical Ethical Committee of School of Stomatology, Zhejiang University and written informed consent was obtained from each individual patient before teeth extraction. Three healthy human premolars, which were extracted for orthodontic reasons, were used for periodontal ligament (PDL) cell isolation. Human PDL tissue was scraped from the middle third of the root surface and isolated as described by our previous article (Wang et al., 2017, 2020). Briefly, Each PDL tissue sample was cut into 1–2 mm fragments, and each fragment was placed in a T25 culture flask containing a minimum amount of α-MEM medium and 10% fetal bovine serum. The medium was changed every three days until the growth of PDLSC was observed. The cells (Supplementary Figure 1) were detached with trypsin-EDTA (0.25%) and cultured into fresh plates.

For adipogenic differentiation assays (Cyagen, China), cells were exposed to an adipogenic medium; the medium was a-MEM containing 10% FBS, 2 mM dexamethasone, 0.2 mM indomethacin, 0.01 g/L insulin, and 0.5 mM isobutyl-methylxanthine. The adipogenic medium was refreshed every 3 days. After 4 weeks, cells were stained with oil red O (Supplementary Figure 1E).

For chondrocyte differentiation (Cyagen, China), the medium contained a-MEM supplemented with 10% FBS, 2 ng/mL transforming growth factor-b1, ITS + Premix, 50 mg/mL L-ascorbic acid, 100 mg/mL sodium pyruvate, 100 nM dexamethasone, and 100 U/mL penicillin/streptomycin. After 4 weeks, chondrogenic differentiation was assessed via Alcian Blue (Solarbio) staining (Supplementary Figure 1C).

The inhibitor NC, miR584-5p inhibitor, mimics NC and miR584-5p mimics were synthesized by Sangon (Shanghai, China). PDLSCs were seeded at a density of 5 × 10⁴ cells/ml. The inhibitor NC, miR584-5p inhibitor, mimics NC and miR584-5p mimics were transfected into cells using jetPRIME according to manufacturer’s instructions (Invitrogen, United States). Total 50 nM artificial miRNA, 50 μL buffer and 0.5 μL jetPRIME per well were used in a 12-well plate. The medium was replaced after 24 h transfection.

Alizarin Red and ALP Staining Assay

Periodontal ligament stem cells were plated into 12-well plates at 2 × 10⁵ cells per well for ALP staining assays and at 1.5 × 10⁵ cells per well for Alizarin Red staining assays. Cells were cultured in the osteogenic medium, which contains Alpha Modification Minimum Essential Medium Eagle (α-MEM) containing 10% Fetal Bovine Serum (FBS, GEMINI Bio, Liverpool, United Kingdom), 50 μM ascorbic acid, 10 μM Dexamethasone and 10 mM β-glycerophosphate. For ALP staining assays, PDLSCs were stained at 3 and 7 days with an alkaline phosphatase kit (Beyotime, China). For Alizarin Red staining assays, PDLSCs were stained with 2% alizarin red S (ScienCell, United States) at 14 and 21 days.

RNA Isolation and Real-Time qPCR Analysis

After the cells were rinsed with PBS at 4°C, total mRNA was extracted using Trizol^® reagent and reverse transcribed into cDNA by PrimeScript^TM RT reagent Kit (RR037A). For miRNAs analysis, miRNA were reverse transcribed to cDNAs using the miRNA First-Strand cDNA Synthesis (Sangon, China) with the Stem-loop Method by primers (Supplementary Table 2). RT- qPCR was performed using SYBR^® Primix Ex Taq^TM by primers (Supplementary Table 3). U6 small nuclear RNA and GAPDH was used as an endogenous control for miRNA and mRNA. Relative differences in the PCR product amounts were evaluated by using the 2^–ΔΔCT method.

Protein Extraction and Western Blot Analysis

Periodontal ligament stem cells were washed with cold PBS and lysed with RIPA containing 1% PMSF at 4°C for 30 min. Then cells were centrifuged at 4°C and 12,000 rpm for 30 min. For western blot analysis, 20 μg of each protein was separated on 10% SDS-polyacrylamide gels and transferred onto polyvinylidene difluoride membranes (Millipore, Bedford, MA, United States). The membranes were blocked with Tris-buffered saline containing 0.1% Tween-20 (TBST) and 5% fat-free milk, then incubated overnight at 4°C with primary antibodies against ALPL (1:1,000) (11187-1-AP), RUNX2 (1:1,000) (ab92336), SP7(1:1,000) (ab209484), COL1A1(1:500) (ab6308) and GAPDH (1:5,000). The membranes were washed three times with TBST and incubated with HRP-conjugated secondary antibodies (1:5,000) for 1 h at room temperature. The bands were visualized using enhanced chemiluminescence reagents (Millipore). The bands were measured using ImageJ software (National Institutes of Health, Bethesda, MD, United States).

Statistical Analysis

Statistical analysis was performed using GraphPad prism 8. Student’s t-test or two-way ANOVA was used to analyze significant differences between groups. A P-value less than 0.05 was considered statistically significant. For biological replicates, three sets of cell-based assays were performed.

Results

Identification of Differentially Expressed Genes During Periodontal Ligament Stem Cells Osteogenic Differentiation

For GSE99958, we used the package CORNAS (Low et al., 2017), which is specially used to analyze non-duplicate samples, and got 276 up-regulated genes and 285 down-regulated genes (Figure 1A). For GSE159507, the limma package was conducted to select 1,908 up-regulated genes and 2,174 down-regulated genes (Figure 1C). LogFC > 1 or LogFC < −1 and P-value < 0.05 were considered to indicate statistical significance. The representative DEGs are shown in a heat map (Figures 1B,D). Then 62 up-regulated genes and 52 down-regulated genes were identified by intersecting DEGs between GSE99958 and GSE159507 (Figures 2A,C).

FIGURE 1

Figure 1. Analysis of differentially expressed genes (DEGs) of GSE99958 and GSE159507. For GSE99958, the CONORS package was used to analyze the differentially expressed genes (DEGs) between the Group induced 3, 7, 14 days and the Group control. A venn graph (A) shows the total 276 up-regulated and 285 down-regulated DEGs and a heat map (B) shows representative DEGs. The number in the heat map grid represents multiples of gene counts. For GSE159507, the limma package was used to analyze the difference genes between the Group induced and the Group control. A volcano plot (C) and heat map (D) reveal the 4,082 DEGs (total 1,908 up-regulated and 2,174 down-regulated DEGs). Red color indicates up-regulated genes, and blue indicates down-regulated genes. Group Control: PDLSCs treated with control medium. Group induced: PDLSCs treated with osteogenic medium.

FIGURE 2

Figure 2. Identification and GO enrichment analysis of common differentially expressed genes. After the common DEGs of GSE99958 and GSE159507 were determined, there were 62 up-regulated genes and 52 down-regulated genes (A,C). Functional enrichment analysis shows biological processes (BP), cellular components (CC), and molecular functions (MF) enrichment from analyses of the 62 up-regulated genes (B) and 52 down-regulated genes (D). P value is 0.05. The GO enrichment analysis of up-regulated genes is mainly concentrated in “extracellular matrix organization,” “extracellular structure organization,” and “ossification.” The GO enrichment of down-regulated expressed genes mainly includes “morphogenesis of a branching structure,” “morphogenesis of a branching epithelium,” and “activation of protein kinase activity.”

Gene Ontology and KEGG Enrichment Analyses of Different Expression miRNAs

Differentially expressed genes were separately enriched in GO terms, and the top most significant terms of up-regulated and down-regulated DEGs are shown (Figures 2B,D). Moreover, up-regulated DEGs were enriched in “extracellular matrix organization,” “extracellular structure organization,” and “ossification” Down-regulated DEGs were involved in “morphogenesis of a branching structure,” “morphogenesis of a branching epithelium,” and “activation of protein kinase activity.” The results of KEGG and detailed GO enrichment analyses of DEGs are provided in Supplementary Figure 2 and Supplementary Table 1.

Protein-Protein Interaction Network Construction and Module Analysis

The PPI network was analyzed using STRING (see text footnote 2). Analysis of functional interactions between 114 DEGs was performed in order to elucidate the mechanisms of osteogenic differentiation of PDLSCs. An interaction with a combined score was selected and used to construct a PPI network with Cytoscape software (Figures 3A,B). The top 10 up-regulated hub genes (APOB, ASPN, ALPL, COL3A1, A2M, COL4A5, OMD, COL11A1, ZBTB16, and CX3CL1 and down-regulated hub genes (FGF2, NGF, NRG1, FST, EDN1, SEMA3A, TNC, NEFM, GRIN2A, and NPTX1) were screened out separately through CytoHubba, a plug-in of Cytoscape (Figures 3C,D).

FIGURE 3

Figure 3. PPI network construction and hub genes screening. PPI network construction and module analysis (A) and (B) shows the PPI network of DEGs. The 62 up-regulated genes (35 nodes, 46 edges) are marked in red, while the 52 down-regulated genes (19 nodes, 21 edges) are marked in blue. The 10 up-regulated hub genes (C) and 10 down-regulated hub genes (D) were identified in the densest connected regions with the Degree algorithm, using cytoHubba. The score is indicated in red color. Darker color indicates a higher score.

Identification of Different Expression miRNAs During Periodontal Ligament Stem Cells Osteogenic Differentiation and miRNA-mRNA Network Prediction

For GSE159508, we used the limma package to obtain six up-regulated miRNAs (LogFC > 1 and P-value < 0.05) and 10 down-regulated miRNAs (LogFC < −1 and P-value < 0.05) (Figures 4A,B). Using the miRwalk website (see text footnote 4), the predicted target genes of DEMs were determined. After the predicted target genes of DEMs and the obtained hub genes were intersected, a negatively regulated PPI network of miRNA-mRNA was obtained (Figures 4C,D). We obtained five up-regulated (miR-337-3p, miR-376c-3p, miR-4288, miR-483-5p, and miR-654-3p) and 6 down-regulated (miR-25-5p, miR-3940-5p, miR-584-5p, miR-642b-3p, miR-663a, and miR-874-3p) DEMs that potentially interacted to hub genes.

FIGURE 4

Figure 4. Differentially expressed miRNA analysis and miRNA-mRNA potential interaction network construction. For the miRNA-seq dataset GSE159508, limma package was used to analyze the DEMs between the Group induced and the Group Control. A volcano plot (A) and heat map (B) shows the 16 DEMs. Red color indicates six up-regulated miRNA, and blue indicates 10 down-regulated miRNA. The predicted target genes of DEMs is calculated by miRWalk website. After the predicted target genes of DEMs and the obtained hub genes were intersected, a negatively regulated PPI of miRNA-mRNA network was obtained (C,D). We obtained five up-regulated (miR-337-3p, miR-376c-3p, miR-4288, miR-483-5p, and miR-654-3p) and 6 down-regulated (miR-25-5p, miR-3940-5p, miR-584-5p, miR-642b-3p, miR-663a, and miR-874-3p) DEMs that potentially interacted to hub genes. Group Control: PDLSCs treated with control medium. Group induced: PDLSCs treated with osteogenic medium.

Hub Genes and Key Different Expression miRNAs Change During the Osteogenic Differentiation of Periodontal Ligament Stem Cells by qPCR

Periodontal ligament stem cells were isolated, passaged, and identified (Supplementary Figure 1). At 3, 7, 14, and 21 days after osteogenic differentiation, ALP and ARS staining confirmed the osteogenic differentiation of PDLSCs (Figure 5A). The result of qPCR (Figures 5B,C) shows that the expression of 10 up-regulated DEGs (ALPL, COL4A5, COL3A1, COL11A, OMD, APOB, ZBTB16, ASPN, CX3CL, and A2M) and 10 down-regulated DEGs (NEFM, FST, TNC, NPTX1, EDN1, FGF2, SEMA3A, NGF, NRG1, and GRIN2A). The qPCR results of miRNA reveals that the expression of 6 down-regulated DEMs (miR-25-5p, miR-3940-5p, miR-584-5p, miR-642b-3p, miR-663a, and miR-874-3p) and three up-regulated DEMs (miR-337-3p, miR-376c-3p, and miR-483-5p) was consistent with the analysis results of miRNA-seq (Figures 6A,B). The expression of miR-4288 and miR-654-3p were increased at day 3 and decreased at day 7 and day 14 during osteogenic differentiation. The possible change patterns of hub genes and miRNAs verified during PDLSCs osteogenic differentiation in vitro are shown in a schematic illustration (Figure 6C).

FIGURE 5

Figure 5. Hub genes change during osteogenic differentiation of PDLSCs in vitro. PDLSCs cells were cultured in osteogenic induction medium for 0, 3, 7, 14 days and tested by qPCR. (A) Representative microscope pictures show results of ALP staining at 3 and 7 days, and the results of Alizarin Red staining at 14 and 21 days. ALP and ARS staining were used to prove that the process of osteogenic differentiation was successful. (B) qPCR results shows up-regulated hub genes APOB, ASPN, ALPL, COL3A1, A2M, COL4A5, OMD, COL11A1, ZBTB16, and CX3CL1 (C) down-regulated hub genes FGF2, NGF, NRG1, FST, EDN1, SEMA3A, TNC, NEFM, GRIN2A, and NPTX1 during osteogenic differentiation of PDLSCs. Bar = 25 μm. *p < 0.05.

FIGURE 6

Figure 6. Candidate miRNAs change during osteogenic differentiation of PDLSCs in vitro. PDLSCs cells were cultured in osteogenic induction medium for 0, 3, 7, 14 days and reverse transcribed with Stem-loop Method and then tested by qPCR. The histogram shows qPCR results of five up-regulated miRNAs (A) and 6 down-regulated miRNAs (B). Demonstration diagram of hub genes and miRNAs change patterns verified during PDLSCs osteogenic differentiation in vitro (C). *p < 0.05.

miR-584-5p Regulates the Osteogenic Differentiation of Periodontal Ligament Stem Cells

We predict that miR584-5p can bind to ALPL and may regulate the expression of ALPL. ALPL is a very important osteogenic regulatory gene. Tian et al. (2020) found that inhibiting miR-584-5p in dental pulp stem cells can promote cell proliferation. Therefore, we chose miR-584-5p for follow-up experiments. In order to verify the regulation of miR-584-5p on the osteogenic differentiation of PDLSCs. PDLSCs were transfected with inhibitor NC, miR584-5p inhibitor, mimics NC and miR584-5p mimics. MiR-584-5p inhibitor or miR-584-5p mimics effectively reduced or increased the expression of miR-584-5p, respectively (Figure 7A). Cells were cultured in the osteogenic induction medium for 7, 14 days. ALP activity and the formation of mineralized nodules, respectively, were detected by ALP staining and Alizarin red staining (Figure 7B). Osteogenesis-related indicators ALPL, BMP2, RUNX2, and OCN were used to further prove the osteogenic regulation of miR-584-5p (Figures 7C,D). The osteogenic differentiation-related protein expression was detected by WB assay in different groups at day 7 (Figure 8A). The semi-quantitative results of related proteins are calculated and displayed, ALPL, SP7, RUNX2, and COL1A1 (Figures 8B–E).

FIGURE 7

Figure 7. MiR-584-5p regulates the osteogenic differentiation of PDLSCs in vitro. PDLSCs were transfected with inhibitor NC, miR584-5p inhibitor, mimics NC and miR584-5p mimics. After 24 h of transfection, the relative expression of miR584-5p was detected by qPCR (A). Cells were cultured in osteogenic induction medium for 7, 14 days. ALP staining (7 days) and Alizarin Red staining (14 days) were used to detect ALP activity and the formation of mineralized nodules, respectively (B). Osteogenesis related indicators ALPL, BMP2, RUNX2 and OCN were used to further prove the osteogenic regulation of miR-584-5p (C,D). **p < 0.01 and ***p < 0.001.

FIGURE 8

Figure 8. MiR-584-5p regulates osteogenic differentiation-related protein expression in vitro. (A) The osteogenic differentiation-related protein expression was detected by WB assay in different groups at day 7. The semi-quantitative results of related proteins are calculated and displayed, ALPL (B), SP7 (C), RUNX2 (D), and COL1A1 (E). ***p < 0.001.

Discussion

In this study, through the DEGs analysis of multiple datasets, the gene changes of PDLSCs after osteoinduction were obtained, and the GO enrichment analysis of DEGs was performed. Fan reported (Fan et al., 2020) that the main enrichment pathways of BMSCs osteogenic differentiation were extracellular matrix organization, ossification, negative regulation of cell proliferation, vasculature development, cell division and spindle. Wang reported (Yang et al., 2019) that regulation of the cellular amino acid metabolic process, DNA metabolic process and regulation of apoptotic process were mostly enriched during the osteogenic differentiation of BMSCs. Some pathways obtained in this study overlap with the enrichment pathways of BMSCs. Our study found that extracellular matrix organization, extracellular structure organization and ossification played an important role in the osteogenic differentiation of PDLSCs. The process of growth and osteogenic differentiation of pluripotent stem cells mainly includes cell adhesion, proliferation, extracellular matrix secretion and maturation, and extracellular matrix mineralization (Stein et al., 2004). The enriched function pathways derived from the analysis of these DEGs can well meet the functional requirements needed in the process of osteogenic differentiation.

A total of 114 DEGs were screened in the two datasets, and 62 up-regulated genes and 52 down-regulated genes were identified. Then, 10 up-regulated hub genes and 10 down-regulated hub genes were screened out through Cytoscape. Among these 10 up-regulated hub genes, most of them such as ALPL, COL4A5 (Sun J. et al., 2021), COL3A1 (Minaříková et al., 2015; Williams et al., 2018), COL11A1 and OMD (Lin et al., 2019, 2021; Ustriyana et al., 2021) are reported to be related to bone mineralization. This result partly proves that DEGs analysis is effective. Other up-regulated hub genes are APOB, ZBTB16, ASPN, CX3CL1, and A2M. ZBTB16 has been reported to be significantly elevated in the osteogenic differentiation of human mesenchymal stem cells (MSCs) (Ikeda et al., 2005), and promote the osteogenic differentiation of dental follicle–derived precursor cells (DFCs), another type of odontogenic pluripotent stem cells (Felthaus et al., 2014). APOB gene is related to the regulation of lipid metabolism (Seo et al., 2004), and there is no literature report that it is directly involved in the regulation of osteogenic differentiation. In addition, APOB was also strongly upregulated in the osteogenic differentiation of canine dental pulp stem cells (Sirirat et al., 2020).

Several studies point out that ASPN was increased in the osteogenic differentiation of PDLSCs (Yamada et al., 2006; Ueda et al., 2016; Garna et al., 2020), and ASPN is a specific molecule marker of PDLSCs. However, overexpression of ASPN in PDLSCs alone cannot promote the osteogenic differentiation of PDLSCs, but inhibits the osteogenic differentiation of PDLSCs partly (Yamada et al., 2007). In addition, ASPN accelerated bone resorption at the orthodontic tension side (Zhang et al., 2020). This contradictory result may be related to the different subgroups of PDLSCs, or to the balance of bone formation and resorption. The molecular mechanism and effects of ASPN in the osteogenic differentiation of PDLSCs needs more research.

The expression of osteoblasts CX3CL1 is critical for osteoclast differentiation. CX3CL1 plays an important role in bone formation and resorption balance by playing dual functions as a chemotactic factor and adhesion molecule for osteoclast precursors expressing CX3CR1 (Koizumi et al., 2009; Xiao et al., 2009; Hoshino et al., 2013). This study found that the expression of CX3CL1 increased during the osteogenic differentiation of PDLSCs. There are two possible explanations, PDLSCs express CX3CL1 that may allow chemotactic migration of osteoclasts to bone remodeling area. Another possibility is PDLSCs which also express the receptor CX3CR1, and secreted CX3CL1 can homing more PDLSCs even BMSCs (Zhang et al., 2015) to the position of bone remodeling through paracrine. The specific effects of CX3CL1 during the osteogenic differentiation of PDLSCs has not been elucidated, and the CX3CL1-CX3CR1 axis may serve as a potential chemokine pathway in order to improve the therapeutic efficacy of these cells.

A recent study (Ghale-Noie et al., 2018) reported that serum protein level of A2M was raised in avascular necrosis of the femoral head patients. A2M also has several roles in fibrinolysis and the coagulation cascade (de Boer et al., 1993). Sadeghi showed that ACTH promoted the osteogenic differentiation of MSCs by up-regulation of A2M (Sadeghi et al., 2020).

Down-regulated 10 hub genes was studied and analyzed by searching the published literature. NEFM, GRIN2A, NPTX1, NRG1, and NGF genes are related to nerve cell development and nerve signal transmission. NEFM is associated with neural maturation (Urrutia et al., 2019; Kawase-Koga et al., 2020) and is down-regulated after early adipogenic differentiation (Urrutia et al., 2019; Marcon et al., 2020; Sadeghi et al., 2020; Nie et al., 2021). GRIN2A is found in nerve cells (neurons) in the brain and spinal cord and one component of a subset of NMDA receptors. NPTX1 may be involved in mediating the uptake of synaptic material during synapse remodeling. Chorion Stromal Cells express higher NPTX1 than Amnion Stromal Cells (Jones et al., 2021). NRG1 is one of four proteins in the neuregulin family that act on the EGFR family of receptors. It was decreased during BMSCs chondrocyte differentiation (Maumus et al., 2020). Nerve growth factor (NGF) is a neurotrophic factor and neuropeptide primarily and was reported to mediate osteogenic differentiation of BMSCs (Peyton, 2017). Interestingly, PDLSCs were able to induce neural progenitor differentiation.

FST was decreased during DPSCs osteogenic differentiation [62] and FST did not enhance osteogenic differentiation of MSCs, but increased committed osteoblast mineralization (Fahmy-Garcia et al., 2019). TNC, as an Extracellular protein, helps osteoblast adhesion to the ECM by creating cell-matrix adhesion sites (Baroncelli et al., 2018). TNC is regarded as one of the markers of the tendon extracellular matrix (Liu et al., 2018). EDN1 is a potent vasoconstrictor and produced by vascular endothelial cells. EDN-1 overexpressed BMSCs showed increased proliferation and significantly increased osteogenesis potential ability (Hu et al., 2017). The overexpression of EDN1 in ADSCs significantly promoted the proliferation and migration of co-cultured HUVECs (Wang et al., 2021). FGF2 inhibited the osteogenic differentiation of PDLSCs. For BMSCs, FGF2 has no obvious function in improving the osteogenic-related genes, but it can ameliorate the impaired osteogenesis (Hao et al., 2020). However, some authors have reported different results. FGF2 enhances the proliferation and osteogenic differentiation of MSCs to promote bone formation (Li et al., 2017; Poudel et al., 2019). Autoregulation of osteocyte SEMA3A promotes bone formation and counteracts bone aging (Hayashi et al., 2019) and is a factor that accelerates osteogenic differentiation (Qiao et al., 2018; Sun Z. et al., 2021). In summary, the effects of FST, TNC, EDN1, FGF2, and SEMA3A in the osteogenic differentiation of PDLSCs need more research and work.

It has been reported that miRNA can regulate the osteogenic differentiation of BMSCs and PDLSCs. For BMSC, MiRNAs control gene expression in osteogenic differentiation by regulating two crucial signaling cascades in osteogenesis: the transforming growth factor-beta (TGF-β)/bone morphogenic protein (BMP) and the Wingless/Int-1(Wnt)/β- catenin signaling pathways (Mazziotta et al., 2021). In previous literature, miR-17, miR-21, miR23a, miR-24, miR-132, miR-138, miR-214, miR-218, miR-543, and miR1305 have been reported to regulate PDLSCs Osteogenic differentiation (Ahmad et al., 2021). In these documents, the target genes of miRNA are mostly concentrated in the WNT signaling pathway, such as TCF3 (Liu et al., 2013), CTNNB1 (Cao et al., 2017), or the TGF-β pathway (Wei et al., 2015). In addition, miR-1305 (Chen and Liu, 2017) and miR-218 (Gay et al., 2014) have been reported to directly regulate the expression of RUNX2. The research results of previous scholars provide us with ideas for further searching for the target genes of DEMs. The miRNA-seq analysis obtained 11 DEMs. The increased miR-654-3p, miR-4288 and decreased miR-663a and miR-874-3p have the same trends as reported in the literature during the osteogenic differentiation of PDLSCs (Hao et al., 2017) (GSE106887). In addition, miR-584-5p was also reduced after osteogenic differentiation of DPSCs (Liu et al., 2020). MiR-337-3p upregulated in DPSCs while downregulated in BMSCs during osteogenic differentiation (Gaus et al., 2021). The role and expression trend of these miRNAs are consistent with our analysis results. However, MiR-483–5p is involved in the pathogenesis of osteoporosis by promoting osteoclast differentiation (Li et al., 2020). MiR-376c-3p inhibits BMSCs osteogenesis (Camp et al., 2018; Kureel et al., 2018). The role of these two miRNAs needs further study. There is no literature describing miR-25-5p, miR-642b-3p, and miR-3940-5p directly related to osteogenesis.

The sample size of the datasets we included is not very large. In addition, we only selected miR-584-5p as an in-depth study and did not prove the exact target gene of miR-584. This is the shortcoming of our experiment. In this study, through the analysis of multiple datasets, the expression changes of nine miRNAs and 18 hub genes were obtained. This provides new preliminary target genes for studying the osteogenic differentiation of PDLSCs.

Data Availability Statement

The datasets GSE99958, GSE159507, and GSE159508 for this study can be found in the (Gene Expression Omnibus) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi).

Author Contributions

CW, LD, ZJ, JZ, and YW contributed in the material preparation, collected and analyzed the data. CW and GY wrote the first draft of the manuscript. All authors contributed to the study conception and design, commented on previous versions of the manuscript, and read and approved the final manuscript.

Funding

This work was supported by grants from the National Natural Science Foundation of China (81801026), the Natural Science Foundation of Zhejiang Province (Q16H140002), and the Health Department Foundation of Zhejiang Province (2018KY499).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank Obwte and Jimmy for their help in R language and bioinformatics analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.750827/full#supplementary-material

Supplementary Figure 1 | Isolation and characterization of human PDLSCs. PDLSCs at P0 showed a typical fibroblast-like spindle appearance (A). Mineralized nodules were stained by ARS at 14d (B), cartilage tissue was stained by Alcian Blue at 21d (C), Clone formation (D), and lipid droplets were stained by Oil Red O staining at 14d (E).

Supplementary Figure 2 | The GOCord circle graph and KEGG enrichment analysis of DEGs. The GOCord circle graph shows the relationship between the GO terms and the up-regulated 62 DEGs (A) or the 52 down-regulated DEGs (B). Each GO term is indicated with a specific color. The (C) and (D) show KEGG enrichment analysis of DEGs. The KEGG enrichment analysis of up-regulated genes (C) is mainly concentrated in “Protein digestion and absorption” and “Thiamine metabolism.” The KEGG enrichment of down-regulated expressed genes (D) mainly includes the “Calcium signaling pathway” and “TGF-beta signaling pathway.”

Supplementary Table 1 | Detailed GO enrichment analyses of DEGs.

Supplementary Table 2 | MiRNA primers used in this study.

Supplementary Table 3 | Hub genes primers used in this study.

Abbreviations

BMSCs, bone marrow-derived mesenchymal stem cells; MSC, mesenchymal stem cell; IBMX, 3-isobutyl-1-methyl-xanthine; α-MEM, alpha Modification Minimum Essential Medium Eagle; FBS, fetal bovine serum; ARS, alizarin red staining; PDLSCs, periodontal ligament stem cells; ADSCs, adipose-derived stem cells; DEMs, differentially expressed miRNAs; DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; PDLCs, periodontal ligament cells; GO, gene ontology; PPI, Protein-protein interaction network; qPCR, quantitative Polymerase Chain Reaction; DFCs, dental follicle derived precursor cells.

Footnotes

References

Abdan, A., Mustafa, A., Kathrin, A. A., Opeyemi, A. R., Filipe, A. d. C., Collins, A., et al., (2020) Analyse der Genexpression von Humanen Stro-1-Positiven Zahnkeim- und Beckenkammzellen in DME-Medium und Osteogenem Differenzierungsmedium. Doctoral thesis. Available online at: https://ediss.uni-goettingen.de/handle/21.11130/00-1735-0000-0005-145D-D?show=full (accessed August 18, 2020).

Ahmad, P., Stoddart, M. J., and Della Bella, E. (2021). The role of noncoding RNAs in osteogenic differentiation of human periodontal ligament stem cells. Craniomaxillof. Traum. Reconstr. Open 6:247275122199922.