Skip to content

Advertisement

BMC Genetics

What do you think about BMC? Take part in

Open Access

Identification and integrated analysis of differentially expressed lncRNAs and circRNAs reveal the potential ceRNA networks during PDLSC osteogenic differentiation

BMC GeneticsBMC series – open, inclusive and trusted201718:100

https://doi.org/10.1186/s12863-017-0569-4

Received: 19 August 2017

Accepted: 16 November 2017

Published: 2 December 2017

Abstract

Background

Researchers have been exploring the molecular mechanisms underlying the control of periodontal ligament stem cell (PDLSC) osteogenic differentiation. Recently, long noncoding RNAs (lncRNAs) and circular RNAs (circRNAs) were shown to function as competitive endogenous RNAs (ceRNAs) to regulate the effect of microRNAs (miRNAs) on their target genes during cell differentiation. However, comprehensive identification and integrated analysis of lncRNAs and circRNAs acting as ceRNAs during PDLSC osteogenic differentiation have not been performed.

Results

PDLSCs were derived from healthy human periodontal ligament and cultured separately with osteogenic induction and normal media for 7 days. Cultured PDLSCs were positive for STRO-1 and CD146 and negative for CD31 and CD45. Osteo-induced PDLSCs showed increased ALP (alkaline phosphatase) activity and up-regulated expression levels of the osteogenesis-related markers ALP, Runt-related transcription factor 2 and osteocalcin. Then, a total of 960 lncRNAs and 1456 circRNAs were found to be differentially expressed by RNA sequencing. The expression profiles of eight lncRNAs and eight circRNAs were measured with quantitative real-time polymerase chain reaction and were shown to agree with the RNA-seq results. Furthermore, the potential functions of lncRNAs and circRNAs as ceRNAs were predicted based on miRanda and were investigated using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis. In total, 147 lncRNAs and 1382 circRNAs were predicted to combine with 148 common miRNAs and compete for miRNA binding sites with 744 messenger RNAs. These mRNAs were predicted to significantly participate in osteoblast differentiation, the MAPK pathway, the Wnt pathway and the signaling pathways regulating pluripotency of stem cells. Among them, lncRNAs coded as TCONS_00212979 and TCONS_00212984, as well as circRNA BANP and circRNA ITCH, might interact with miRNA34a and miRNA146a to regulate PDLSC osteogenic differentiation via the MAPK pathway.

Conclusions

This study comprehensively identified lncRNAs/circRNAs and first integrated their potential ceRNA function during PDLSC osteogenic differentiation. These findings suggest that specific lncRNAs and circRNAs might function as ceRNAs to promote PDLSC osteogenic differentiation and periodontal regeneration.

Keywords

Long non-coding RNAsCircular RNACompetitive endogenous RNAPeriodontal ligament stem cellsOsteogenic differentiationRNA sequencing

Background

Periodontitis is a slowly progressive disease characterized by the loss of periodontal tissue, which is the main cause of tooth loss [1]. With the development of stem-cell delivery therapeutics, periodontal ligament stem cells (PDLSCs) have been shown to produce typical periodontal ligament-like tissue to regenerate tissue damaged by periodontitis [2, 3]. Our previous studies also have shown the great potential of PDLSCs to regenerate periodontal tissue and form a bioengineered tooth root in miniature pigs [46]. The regeneration potency of PDLSCs contributes to their self-renewal and multi-differentiation capacity, especially osteogenic differentiation [7]. Exploring the molecular mechanisms of PDLSC osteogenic differentiation might provide new genetic strategies for periodontal regenerative medicine.

Recently, through microarray analysis, microRNAs (miRNAs) were identified and predicted to be involved in PDLSC osteogenic differentiation [8]. In addition, miRNA146a, miRNA17 and miRNA22 have been demonstrated to regulate PDLSC osteogenic differentiation by modulating the expression of target genes at the post-transcriptional level [911]. However, recent studies have revealed that a new player, competing endogenous RNA (ceRNA), is essential for the circuitry of miRNAs and target genes [12]. By competing for common miRNA response elements (MREs), ceRNAs can break the balance between miRNAs and target genes to regulate the physiological and pathophysiological process [13]. These ceRNAs include various types of RNAs, such as long non-coding RNAs (lncRNAs), circular RNAs (circRNAs), messenger RNAs (mRNAs) and pseudogenes.

LncRNA is a class of non-coding RNA (ncRNA) transcripts longer than 200 nucleotides [14]. Recent studies have reported that lncRNAs were involved in the osteogenic differentiation process [15]. For example, up-regulation of lncRNA HIF 1α-anti-sense 1 induced by TGFβ-mediated targeting of sirtuin 1 promotes the osteogenic differentiation of human bone marrow stromal cells [16]. Furthermore, lncRNA ANCR was proved to be critical in regulating the PDLSC osteogenic differentiation via the Wnt signaling pathway [17].

CircRNA is another new class of RNA composed of more than one exon with a covalently closed loop [18]. Compared to linear RNA, this circular structure is more stable and resistant to RNase R [19]. Emerging evidence has revealed that circRNAs participate in osteogenic differentiation [20, 21]. For instance, 154 differentially expressed circRNAs were found to associate with bone morphogenetic protein 2-induced osteogenic differentiation of MC3T3-E1 cells [22]. Moreover, circRNAs were predicted to have potential roles in osteogenesis of PDLSCs [23].

Although several lncRNAs, circRNAs and miRNAs are suggested to be involved in PDLSC osteogenic differentiation, there is not much data published on their potential networks and functions. To fully understand the impact of ceRNA crosstalk on PDLSC osteogenic differentiation, it will be crucial to integrate the lncRNA/circRNA-miRNA-mRNA competitive regulatory networks. In this study, we developed RNA sequencing (RNA-seq) with Illumina HiSeq2000 to comprehensively identify differentially expressed lncRNAs and circRNAs in normal and osteogenic inductive PDLSCs. Subsequently, the representative lncRNAs and circRNAs were further confirmed using quantitative real-time polymerase chain reaction (qRT-PCR). Finally, we predicted the ceRNA networks of lncRNAs, circRNAs, miRNAs and mRNAs based on miRanda and investigated their potential regulatory roles via gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Our findings might provide new evidence for exploring the molecular mechanism of PDLSC osteogenic differentiation.

Methods

All protocols for the handling of human tissues were approved by the Research Ethics Committee of Stomatology Hospital of Shandong University, China (G201401601). Informed consent was obtained from all donors.

Cell culture and osteogenic induction

In this study, PDLSCs were derived from the middle third of the root surface of normal human impacted third molars, which were collected from 18- to 20-year-old patients at the Department of Oral Maxillofacial Surgery, Stomatology Hospital of Shandong University, using the explant culture method. Then, they were cultured with normal media, consisting of α-modification of Eagle’s media (HyClone, South Logan, UT, USA), 10% fetal calf serum (Gibco BRL, Grand Island, NY, USA), 100 U/ml penicillin and 100 μg/ml streptomycin (Invitrogen, Carlsbad, CA, USA) at 37 °C in 5% carbon dioxide. The stemness of PDLSCs was characterized by the scanning of cell surface markers (STRO-1, CD146, CD31, and CD45) through flow cytometric analysis (Becton, Dickinson and Company, Franklin Lakes, NJ, USA). For osteogenic differentiation, PDLSCs were cultured with osteogenic inductive media supplemented with 10 nM dexamethasone, 10 mM β-glycerophosphate and 50 μg/ml vitamin C (Sigma-Aldrich, St. Louis, MO, USA). Through separately culture with osteogenic induction and normal media for 7 days, PDLSCs were divided into two groups: induced and non-induced groups. All cells in this study were used at passage number 3.

Alkaline Phosphatase (ALP) staining and ALP activity assay

ALP staining and activity assay were performed using an Alkaline Phosphatase Kit (Sigma-Aldrich) as described previously [24]. Briefly, PDLSCs were fixed with 70% ethanol for 30 min and stained with a solution of sodium nitrite, FRV alkaline and naphthol AS-BI alkaline for 15 min away from light. For the ALP activity assays, total protein was isolated from two groups of PDLSCs, incubated for 15 min with citrate buffer and phosphatase substrate (Sigma-Aldrich), and then quantified by spectrophotometric absorbance at 405 nm.

Quantitative real-time PCR (qRT-PCR)

Total RNA was isolated from two groups of PDLSCs. For analysis of linear transcripts, 1 μg of RNA per sample was reverse transcribed into cDNA using a cDNA Reverse Transcription Kit (Takara, Tokyo, Japan). Convergent primers were designed to detect lncRNAs and mRNAs. Divergent primers were designed to detect the circular form of circRNAs. Relative transcript levels were measured with quantitative PCR using a Roche LightCycler®480 sequence detection system (Roche, Basel, Switzerland) following the manufacturer’s protocol. Each 20 μl reaction volume contained 10 μl SYBR® Premix Ex Taq™ (Takara), 0.4 μl 10 μM forward primer (0.4 μM final), 0.4 μl 10 μM reverse primer (0.4 μM final), 200 ng of template cDNA and DEPC-treated water. GAPDH was used as an internal control to quantify and normalize the results. The primer pairs are listed in Additional file 1. The specificity of the reaction was determined by detection of the Tms of the amplification products immediately after the last reaction cycle. The 2-CT value was used for comparative quantitation. All qRT-PCRs were performed in triplicate.

Construction of cDNA libraries and high-throughput sequencing

Total RNA was extracted from two groups using Trizol (Invitrogen) according to the manufacturer’s protocols. Strand-specific cDNA libraries were constructed following a previously described protocol [25] and were sequenced using an Illumina HiSeq2000 sequencer (LC Biotech, Hangzhou, China) by performing a paired-end run with a 100 bp read length. The raw reads were processed by removing the adaptor reads and low-quality tags. All subsequent analyses were performed using clean reads.

Identification of differentially expressed lncRNAs, circRNAs and mRNAs

The clean reads from two cDNA libraries were mapped to the human genome sequence in GENCODE Release 19 using TopHat version 2.0.9 [26]. The transcripts were then assembled and annotated using Cufflinks [27]. According to the annotation of the human genome sequence, the known lncRNAs and mRNAs were identified. The coding potential of the remaining unknown transcripts was calculated with a coding potential calculator based on the quality, completeness, and sequence similarity of the open reading frame to the proteins in the protein databases [28]. The remaining unknown transcripts of more than 200 base pairs (bp) were considered novel lncRNAs with a coding potential score of less than −1, which suggested that a transcript had no protein-coding capacity. The expression levels of lncRNAs and mRNAs were quantified as FPKM (fragments per kilobase of exon per million fragments mapped) using the Cuffdiff program [27]. The differential expression of lncRNAs and mRNAs was determined using DESeq, with P < 0.05 and fold change (FC) ≥ 2 [29].

Compared with linear RNAs, circRNAs exhibit distinct patterns of alternative back-splicing and alternative splicing. An upgraded computational pipeline (CIRCexplorer2) was used to systematically identify and annotate circRNAs [30]. The expression levels of circRNAs were quantified as RPM (mapped backsplicing junction reads per million mapped reads) using the EBSeq package [31]. The differential expression of circRNAs was determined with P < 0.05, FC ≥ 2, and circRNA junction reads ≥5 [31].

Functional analysis

The ceRNA networks among lncRNAs, circRNAs, miRNAs and mRNAs were predicted based on miRanda with a maximum binding free energy of −20 [32]. First, we predicted and selected miRNAs binding with both differentially expressed lncRNAs and circRNAs. Then, target mRNAs of these selected miRNAs were predicted and compared to the differentially expressed mRNAs that were identified in the RNA-seq results. Subsequently, we selected the intersecting elements of target mRNAs and differentially expressed mRNAs to analyze their potential functions through GO functional annotation and KEGG pathway analysis. GO terms were enriched using Blast2GO [33] by referring to GO databases. Meanwhile, KEGG pathway analysis was performed by referring to KEGG pathway databases. Cytoscape3.5.1 was used to display the lncRNA/circRNA-miRNA-mRNA networks.

Statistical analysis

Quantitative qRT-PCR datasets are presented as the means ± standard deviation (SD) of at least three independent experiments. The statistical calculations were performed with SPSS statistics software version 17.0. Student’s t-test was performed for normally distributed data to determine statistical significance. A P-value less than 0.05 was considered statistically significant.

Results

Identification and Osteogenic differentiation of PDLSCs

PDLSCs derived from periodontal ligament explants were cultured with normal media to passage number 3 (Fig. 1a, b). Cultured PDLSCs were positive for STRO-1 and CD146 and negative for CD31 and CD45 (Fig. 1c-f). Increased ALP activity identified via ALP staining and ALP activity assay indicated osteogenic differentiation of osteo-induced PDLSCs (Fig. 1g, h). Subsequently, the up-regulated expression levels of the osteogenesis-related markers ALP, Runt-related transcription factor 2 (Runx2) and osteocalcin (OCN) provided further proof for the occurrence of PDLSC osteogenic differentiation (Fig. 2i-k). These findings agreed with previous reports on PDLSC differentiation into osteoblasts [7].
Fig. 1

Identification and osteogenic differentiation of PDLSCs. a PDLSCs were derived from periodontal ligament explants. b PDLSCs were cultured with non-induced media at passage number 3. c-f PDLSCs were positive for STRO-1 and CD146 and negative for CD31 and CD45. g, h ALP activity was enhanced in osteo-induced PDLSCs (Induced), as evidenced by ALP staining and ALP activity assay. i-k Compared with the non-induced group, the induced group showed up-regulated expression of the osteogenic genes ALP, Runx2 and OCN by qRT-PCR. All PCRs were performed in triplicate. The data are represented as means ± SD. *, P < 0.05; **, P < 0.01

Fig. 2

The apparent variations of differentially expressed lncRNAs and circRNAs. (a, b) LncRNAs and circRNAs were broadly distributed across the 24 pairs of human chromosomes according to their locations. The inner blue ring corresponds to the non-induced group; the outer yellow ring corresponds to the induced group. c Among differentially expressed lncRNAs, 17 common lncRNAs and 180 specific lncRNAs in the non-induced group and 763 specific lncRNAs in the induced group were identified. d Among differentially expressed circRNAs, 95 common circRNAs and 642 specific circRNAs in the non-induced group and 719 specific circRNAs in the induced group were identified. e Differentially expressed lncRNAs, consisting of 777 up-regulated lncRNAs and 183 down-regulated lncRNAs, are displayed in the heatmap. f Differentially expressed circRNAs, consisting of 766 up-regulated circRNAs and 690 down-regulated circRNAs, are displayed in the heatmap

Differential expression of lncRNAs, circRNAs and mRNAs upon osteogenic differentiation of PDLSCs

Ribosome-depleted RNA-seq generated 91.4 million raw reads of the non-induced group and 113.8 million raw reads of the induced group (Table 1). After filtering the adaptor reads and low-quality tags, we separately obtained 90.9 million and 112.4 million clean reads. More than 98% of the raw reads per sample were clean reads. According to the annotation of the human genome sequence GENCODE Release 19, 11,529 lncRNAs (non-induced group) and 9511 lncRNAs (induced group) were identified from the two cDNA libraries (Table 2) (Additional file 2). We also identified 5913 circRNAs in the non-induced group and 3162 circRNAs in the induced group (Table 2) (Additional file 3). In addition, 66,134 mRNAs (non-induced group) and 52,227 mRNAs (induced group) were annotated (Table 2) (Additional file 4). Both lncRNAs and circRNAs were broadly distributed across the 24 pairs of human chromosomes (Fig. 2a, b). A total of 960 lncRNAs, 1161 circRNAs and 1887 mRNAs were found to be differentially expressed, with P-value <0.05 and FC ≥ 2 (Additional file 5, Additional file 6 and Additional file 7). Among the differentially expressed lncRNAs, 17 common lncRNAs and 180 specific lncRNAs in the non-induced group and 763 specific lncRNAs in the induced group were detected, with 777 up-regulated lncRNAs and 183 down-regulated lncRNAs (Fig. 2c). Meanwhile, we also identified 95 common circRNAs and 642 specific circRNAs in the non-induced group and 719 specific circRNAs in the induced group, with 766 up-regulated circRNAs and 690 down-regulated circRNAs (Fig. 2d). The apparent variations in transcripts between the two groups are visually displayed with heatmaps (Fig. 2e, f).
Table 1

Statistical data of RNA-Seq reads for two samples

Sample

Raw reads

Clean reads

Unique lncRNAs

Unique circRNAs

Unique mRNAs

Non-induced

91.4 million

90.9 million

11,529

5913

66,134

Induced

113.8 million

112.4 million

9511

3162

52,227

Table 2

Expression profiles of lncRNAs validated by RNA-seq

Test_id.

Induced (FPKM)

Non-induced (FPKM)

Regulation

TCONS_00019601

93.6452

0

up

TCONS_00227764

830.447

47.6825

up

TCONS_00085268

0

26.6246

down

TCONS_00254538

18,683.1

3962.3

up

TCONS_00198784

9.25995

0.0230801

up

TCONS_00136898

8.67356

0.0105046

up

TCONS_00125934

0

4.65289

down

TCONS_00115113

0

4.05314

down

Differentially expressed lncRNAs/circRNAs validated by qRT-PCR

To verify the results of the RNA-seq experiments, eight lncRNAs and eight circRNAs with P < 0.05, FC ≥ 2 and FPKM/RPM in at least one of the samples ≥4 were selected for qRT-PCR validation. The lncRNAs were amplified with convergent primers and the circRNAs were amplified with divergent primers (Fig. 3a, b). Compared to the non-induced group, the induced group showed increased expression of the lncRNAs coded as TCONS_00019601, TCONS_00227764, TCONS_00254538, TCONS_00198784 and TCONS_00136898 and decreased expression of TCONS_00085268, TCONS_00125934, and TCONS_00115113 (Fig. 3a). The circRNAs named CDR1-AS, NCOA3, and SKIL were up-regulated in the induced group compared to the non-induced group, and the circRNAs IFF01, NTNG1, PLOD2, SMO, and SMURF2 were down-regulated (Fig. 3b). All the results were consistent with the normalized expression of RNA-seq shown in Tables 2 and 3.
Fig. 3

Differentially expressed lncRNAs/circRNAs validated by qRT-PCR. a Convergent primers were designed to detect eight lncRNAs with P < 0.05, FC ≥ 2 and FPKM in at least one of the samples ≥4. LncRNAs coded as TCONS_00019601, TCONS_00227764, TCONS_00254538, TCONS_00198784 and TCONS_00136898 were up-regulated in the induced group compared with the non-induced group, and TCONS_00085268, TCONS_00125934, and TCONS_00115113 were down-regulated. b Divergent primers were designed to detect the circular form of circRNAs with P < 0.05, FC ≥ 2 and RPM in at least one of the samples ≥4. CircRNAs named CDR1-AS, NCOA3, and SKIL were up-regulated in the induced group compared to the non-induced group, and the circRNAs IFF01, NTNG1, PLOD2, SMO, and SMURF2 were down-regulated. The results agreed with the normalized expression of validated lncRNAs and circRNAs shown in Tables 2 and 3. All PCRs were performed in triplicate. The data are represented as means ± SD. *, P < 0.05; **, P < 0.01

Table 3

Expression profiles of circRNAs validated by RNA-seq

Name

Induced (RPM)

Non-induced (RPM)

Regulation

circRNA CDR1-AS

3038.695

1102.345

up

circRNA IFF01

0

337.2772

down

circRNA NCOA3

1421.615

0

up

circRNA NTNG1

0

1286.478

down

circRNA PLOD2

0

3742.131

down

circRNA SKIL

1956.098

0

up

circRNA SMO

0

783.6154

down

circRNA SMURF2

0

282.7267

down

Function of lncRNAs and circRNAs as ceRNAs via lncRNA/circRNA-miRNA-mRNA networks

Based on miRanda with the maximum binding free energy of −20, 430 lncRNAs were predicted to share at least two binding sites with 779 miRNAs (Fig. 4) (Additional file 8). We also predicted that 1401 circRNAs bind 855 miRNAs with at least two binding sites (Additional file 9). Considering that graphics cannot display the enormous amount of network information between 1401 circRNAs and 855 miRNAs, we selected circRNAs with more miRNA binding sites and less binding free energy to make the network diagram (Fig. 5). Through analysis of the common binding MREs of lncRNAs and circRNAs, 165 miRNAs were predicted to combine with both 158 lncRNAs and 1385 circRNAs (data not shown).
Fig. 4

CeRNA networks of 430 lncRNAs and 779 miRNAs based on miRanda

Fig. 5

CeRNA networks of selected circRNAs and miRNAs based on miRanda. Considering that graphics cannot display the enormous amount of network information between 1401 circRNAs and 855 miRNAs, we selected circRNAs with more miRNA binding sites and less binding free energy to make the network diagram using Cytoscape3.5.1

To reveal their potential function, we predicted target mRNAs of these miRNAs based on miRanda and compared these target mRNAs to the differentially expressed mRNAs that were identified in the RNA-seq results (Additional file 7). There were 744 differentially expressed mRNAs that were found to combine with 148 common miRNAs along with 147 lncRNAs and 1382 circRNAs (Additional file 10). The networks of 744 mRNAs and 148 common miRNAs are shown in Fig. 6.
Fig. 6

CeRNA networks of 744 mRNAs and 148 common miRNAs based on miRanda

The potential regulatory roles of the ceRNA networks were predicted by analyzing the functions of 744 mRNAs through GO and KEGG pathway analysis (Additional file 11 and Additional file 12). GO annotations (P < 0.05) involving the top 60 mRNAs are displayed in Fig. 7a and include multiple biological processes, cellular components and molecular functions. Among these GO terms, we obtained GO: 0001649 (osteoblast differentiation), which was significantly enriched by 21 mRNAs (Additional file 11). The complex mRNA networks involved in GO: 0001649 (osteoblast differentiation) and related miRNAs, lncRNAs, and circRNAs are displayed in Fig. 8. Through KEGG analysis, mRNAs were predicted to participte in 34 pathways (Fig. 7b). Among these KEGG pathways, the MAPK pathway, the Wnt pathway and the signaling pathways regulating pluripotency of stem cells were closely related to osteogenic differentiation.
Fig. 7

GO annotations and KEGG pathway analysis of 744 differentially expressed mRNAs. a GO annotations (P < 0.05) involving the top sixty numbers of mRNAs included multiple biological processes, cellular components and molecular functions. b These mRNAs were significantly enriched in 34 KEGG pathways, including the the MAPK pathway, the Wnt pathway and the signaling pathways regulating pluripotency of stem cells

Fig. 8

CeRNA networks of lncRNAs/circRNAs-miRNAs-mRNAs significantly participated in GO: 0001649 (osteoblast differentiation)

Based on the above results, we selected several lncRNAs, circRNAs, miRNAs and mRNAs associated with the MAPK pathway to further display the ceRNA networks (Fig. 9). The lncRNAs coded as TCONS_00212979 and TCONS_00212984, circRNA BANP and circRNA ITCH were predicted to combine with miRNA34a and miRNA146a. DUSP1, FAS and RAC1 were predicted to be target genes of miRNA34a. In addition, PDGFRA, TGFBR2 and MYC were predicted to be targeted by miRNA146a. These six mRNAs were the pivotal genes of the MAPK pathway according to the KEGG analysis. This complicated ceRNA network suggested that TCONS_00212979, TCONS_00212984, circRNA BANP and circRNA ITCH might play regulatory roles in the MAPK pathway through miRNA34a, miRNA146a and their targets during PDLSC osteogenic differentiation.
Fig. 9

CeRNA networks of lncRNAs/circRNAs-miRNAs-mRNAs significantly participated in the MAPK pathway. The lncRNAs coded as TCONS_00212979 and TCONS_00212984, as well as circRNA BANP and circRNA ITCH, were predicted to interact with miRNA34a and miRNA146a. DUSP1, FAS and RAC1 were predicted to be target genes of miRNA34a. PDGFRA, TGFBR2 and MYC were predicted to be targeted by miRNA146a. These six mRNAs are pivotal genes in the MAPK pathway according to KEGG analysis

Discussion

LncRNAs, circRNAs, miRNAs and mRNAs form large-scale ceRNA cross-talk networks through MREs, which has exciting implications for gene regulation at the post-transcriptional level during multiple physiological and pathophysiological processes [12, 34]. In recent years, studies have documented the functions and clinical implications of ceRNAs in cancer, systemic lupus erythematosus and other diseases, which may present opportunities for new therapeutic approaches for diseases [35, 36].

Recently, researchers have systematically constructed ceRNA networks through RNA-seq and bioinformatics in mouse germline stem cells to reveal functions and mechanisms of lncRNAs and circRNAs in mouse germline stem cell self-renewal and differentiation [37]. Moreover, lncRNA POIR was demonstrated to form a regulatory network with miRNA182 and FoxO1 to up-regulate PDLSC osteogenic differentiation in periodontitis patients [38]. However, the ceRNA networks were revealed to be intertwined [39]. To fully understand the impact of ceRNA crosstalk on PDLSC osteogenic differentiation, it will be crucial to integrate the competitive lncRNA/circRNA-miRNA–mRNA regulatory networks. In our study, 744 mRNAs were predicted to combine with 148 common miRNAs, along with 147 lncRNAs and 1382 circRNAs.

Through GO analysis, 21 mRNAs were predicted to significantly participate in osteoblast differentiation (GO: 0001649) (Fig. 8). Among them, ALPL, also called ALP, was reported to be an osteogenesis-related marker and was up-regulated during PDLSC osteogenic differentiation [7]. The up-regulated expression level of ALP was also detected by qRT-PCR in our study (Fig. 1i). SMAD3 and SMAD5, members of the SMAD family, were also predicted to form ceRNA networks and participate in osteoblast differentiation. Both SMAD3 and SMAD5 were demonstrated to regulate PDLSC osteogenic differentiation by modulating TGF-β signals [40, 41]. Additionally, Notch1, another participant in osteoblast differentiation (GO: 0001649), is part of the Notch signaling pathway, which is important for maintaining osteogenic differentiation of PDLSCs [42, 43]. These crucial osteogenic genes formed ceRNA networks with lncRNAs and circRNAs by targeting common miRNAs, and these networks might provide evidence of new regulatory mechanisms in PDLSC osteogenic differentiation.

Through KEGG pathway analysis, mRNAs of ceRNA networks were predicted to be involved in the Wnt pathway, MAPK pathway and signaling pathways regulating the pluripotency of stem cells. Previous studies have demonstrated that Wnt signaling contributes to the differentiation of periodontal ligament fibroblasts into osteoblasts [44]. In addition, the MAPK pathway was found to be critical for skeleton development and bone homeostasis [45]. Moreover, it plays significant roles in osteogenic differentiation of PDLSCs [46, 47]. Furthermore, we constructed a ceRNA network of TCONS_00212979, TCONS_00212984, circRNA BANP, circRNA ITCH, miRNA34a, miRNA146a, DUSP1, FAS, RAC1, PDGFRA, TGFBR2 and MYC. These mRNAs are important elements of the MAPK pathway based on KEGG analysis. Among them, DUSP1, FAS and RAC1 are targeted by miRNA34a, while PDGFRA, TGFBR2 and MYC are targeted by miRNA146a. Studies have illustrated that both miRNA34a and miRNA146a are closely related to osteogenic differentiation of human mesenchymal stem cells [48, 49]. In addition, miRNA146a was revealed to promote differentiation of periodontal ligament cells [9]. TCONS_00212979, TCONS_00212984, circRNA BANP and circRNA ITCH were predicted to bind miRNA34a and miRNA146a. TCONS_00212979, known as CARMEN, has been reported to be a cardiac mesoderm enhancer-associated lncRNA that modulates cardiac differentiation through miRNA143 and miRNA145 [50]. TCONS_00212984 is a novel lncRNA with a genomic origin similar to that of TCONS_00212979 according to RNA-seq. circRNA BANP and circRNA ITCH have both been reported to contribute to carcinogenesis and might serve as cancer biomarkers [51, 52]. However, the regulatory roles of these two lncRNAs and circRNAs in osteogenic differentiation remain unclear. In summary, the ceRNA network suggested that TCONS_00212979, TCONS_00212984, circRNA BANP and circRNA ITCH might interact with miRNA34a and miRNA146a to regulate PDLSC osteogenic differentiation via the MAPK pathway. However, their regulatory mechanisms need to be further investigated. Our future study plan will be to validate their differential expression profiles, verify their ceRNA networks and specify their effects on PDLSC osteogenic differentiation using knockdown and overexpression experiments.

Conclusion

This study identified differentially expressed lncRNAs, circRNAs and mRNAs during osteogenic differentiation of PDLSCs. Competitive lncRNA/circRNA-miRNA–mRNA regulatory networks were comprehensively integrated and predicted to be involved in osteoblast differentiation by GO and KEGG pathway analysis. Moreover, the lncRNAs coded as TCONS_00212979 and TCONS_00212984, circRNA BANP and circRNA ITCH were predicted to interact with miRNA34a and miRNA146a to regulate PDLSC osteogenic differentiation via the MAPK pathway. Our study suggested that specific lncRNAs and circRNAs might function as ceRNAs to promote PDLSC osteogenic differentiation and periodontal regeneration.

Abbreviations

ALP: 

Alkaline phosphatase

bp: 

Base pairs

ceRNAs: 

Competitive endogenous RNAs

circRNAs: 

Circular RNAs

FC: 

Fold change

FPKM: 

Fragments per kilobase of exon per million fragments mapped

GO: 

Gene Ontology

KEGG: 

Kyoto Encyclopedia of Genes and Genomes

lncRNAs: 

Long non-coding RNAs

miRNAs: 

MicroRNAs

MREs: 

miRNA response elements

mRNAs: 

Messenger RNAs

ncRNA: 

Non-coding RNA

OCN: 

Osteocalcin

PDLSC: 

Periodontal ligament stem cell

qRT-PCR: 

Quantitative real-time polymerase chain reaction

RNA-seq: 

RNA sequencing

RPM: 

Mapped back-splice junction reads per million mapped reads

Runx2: 

Runt-related transcription factor 2

SD: 

Standard deviation

Declarations

Acknowledgments

We acknowledge the anonymous reviewers and academic editors for their constructive suggestions on the previous version of the manuscript. We acknowledge laboratory members for their contributions.

Funding

The performance of RNA-seq, analysis of the sequenced data and revision of the manuscript were supported by grant from the National Natural Science Foundation of China (81470709). All the cell biology experiments and molecular biology experiments were supported by grant from the Construction Engineering Special Fund of “Taishan Scholars” (tsqn201611068).

Availability of data and materials

All data analyzed during this study are included in this article and its supplementary information files. RNA-seq data sets are available in the Gene Expression Omnibus database under accession number GSM2856840 (for non-induced PDLSCs) and GSM2856841 (for induced PDLSCs).

Authors’ contributions

FW and DL conceived, designed and revised this study. XG analyzed the RNA-seq data and wrote the manuscript. ML detected the expression of genes using qRT-PCR. YJ performed the cell culture and osteogenic induction treatments. All authors read and approved the final manuscript.

Ethics approval and consent to participate

PDLSCs were derived from the extracted teeth, which were collected from patients at the Department of Oral Maxillofacial Surgery, Stomatology Hospital of Shandong University. All protocols for handling dental tissues were performed in accordance with relevant guidelines and regulations. All protocols for the handling of human tissues were approved by the Research Ethics Committee of Stomatology Hospital of Shandong University, China (G201401601). The participant informed consent was written.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Department of Orthodontics, Shandong Provincial Key Laboratory of Oral Tissue Regeneration, School of Stomatology, Shandong University
(2)
Shandong Provincial Key Laboratory of Oral Tissue Regeneration, School of Stomatology, Shandong University

References

  1. Damgaard C, Holmstrup P, Van Dyke TE, Nielsen CH. The complement system and its role in the pathogenesis of periodontitis: current concepts. J Periodontal Res. 2015;50(3):283–93.View ArticlePubMedGoogle Scholar
  2. Liu Y, Zheng Y, Ding G, Fang D, Zhang C, Bartold PM, Gronthos S, Shi S, Wang S. Periodontal ligament stem cell-mediated treatment for periodontitis in miniature swine. Stem Cells. 2008;26(4):1065–73.View ArticlePubMedPubMed CentralGoogle Scholar
  3. Park JY, Jeon SH, Choung PH. Efficacy of periodontal stem cell transplantation in the treatment of advanced periodontitis. Cell Transplant. 2011;20(2):271–85.View ArticlePubMedGoogle Scholar
  4. Wei F, Qu C, Song T, Ding G, Fan Z, Liu D, Liu Y, Zhang C, Shi S, Wang S, Vitamin C. Treatment promotes mesenchymal stem cell sheet formation and tissue regeneration by elevating telomerase activity. J Cell Physiol. 2012;227(9):3216–24.View ArticlePubMedPubMed CentralGoogle Scholar
  5. Wei F, Song T, Ding G, Xu J, Liu Y, Liu D, Fan Z, Zhang C, Shi S, Wang S. Functional tooth restoration by allogeneic mesenchymal stem cell-based bio-root regeneration in swine. Stem Cells Dev. 2013;22(12):1752–62.View ArticlePubMedPubMed CentralGoogle Scholar
  6. Ding G, Liu Y, Wang W, Wei F, Liu D, Fan Z, An Y, Zhang C, Wang S. Allogeneic periodontal ligament stem cell therapy for periodontitis in swine. Stem Cells. 2010;28(10):1829–38.View ArticlePubMedPubMed CentralGoogle Scholar
  7. Seo BM, Miura M, Gronthos S, Bartold PM, Batouli S, Brahim J, Young M, Robey PG, Wang CY, Shi S. Investigation of multipotent postnatal stem cells from human periodontal ligament. Lancet. 2004;364(9429):149–55.View ArticlePubMedGoogle Scholar
  8. Hao Y, Ge Y, Li J, Hu Y, Wu B, Fang F. Identification of MicroRNAs by microarray analysis and prediction of target genes involved in the Osteogenic differentiation of human periodontal ligament stem cells. J Periodontol. 2017:1–11.Google Scholar
  9. Hung PS, Chen FC, Kuang SH, Kao SY, Lin SC, Chang KW. miR-146a induces differentiation of periodontal ligament cells. J Dent Res. 2010;89(3):252–7.View ArticlePubMedGoogle Scholar
  10. Liu Y, Liu W, Hu C, Xue Z, Wang G, Ding B, Luo H, Tang L, Kong X, Chen X, et al. MiR-17 modulates osteogenic differentiation through a coherent feed-forward loop in mesenchymal stem cells isolated from periodontal ligaments of patients with periodontitis. Stem Cells. 2011;29(11):1804–16.View ArticlePubMedGoogle Scholar
  11. Yan GQ, Wang X, Yang F, Yang ML, Zhang GR, Wang GK. MicroRNA-22 promoted Osteogenic differentiation of human periodontal ligament stem cells by targeting HDAC6. J Cell Biochem. 2017;118(7):1653–8.View ArticlePubMedGoogle Scholar
  12. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146(3):353–8.View ArticlePubMedPubMed CentralGoogle Scholar
  13. Sen R, Ghosal S, Das S, Balti S, Chakrabarti J. Competing endogenous RNA: the key to posttranscriptional regulation. TheScientificWorldJOURNAL. 2014;2014:896206.View ArticlePubMedPubMed CentralGoogle Scholar
  14. Kapranov P, Cheng J, Dike S, Nix DA, Duttagupta R, Willingham AT, Stadler PF, Hertel J, Hackermuller J, Hofacker IL, et al. RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science. 2007;316(5830):1484–8.View ArticlePubMedGoogle Scholar
  15. Fatica A, Bozzoni I. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. 2014;15(1):7–21.View ArticlePubMedGoogle Scholar
  16. WS XY, Tang C, Chen W. Upregulation of long non-coding RNA HIF 1α-anti-sense 1 induced by transforming growth factor-β-mediated targeting of sirtuin 1 promotes osteoblastic differentiation of human bone marrow stromal cells. Mol Med Rep. 2015;12(5):7233–8.View ArticleGoogle Scholar
  17. Jia Q, Jiang W, Ni L. Down-regulated non-coding RNA (lncRNA-ANCR) promotes osteogenic differentiation of periodontal ligament stem cells. Arch Oral Biol. 2015;60(2):234–41.View ArticlePubMedGoogle Scholar
  18. Wilusz JE, Sharp PA. Molecular biology. A circuitous route to noncoding RNA. Science. 2013;340(6131):440–1.View ArticlePubMedPubMed CentralGoogle Scholar
  19. Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015;365(2):141–8.View ArticlePubMedGoogle Scholar
  20. Barrett SP, Salzman J. Circular RNAs: analysis, expression and potential functions. Development. 2016;143(11):1838–47.View ArticlePubMedPubMed CentralGoogle Scholar
  21. Lasda E, Parker R. Circular RNAs: diversity of form and function. RNA. 2014;20(12):1829–42.View ArticlePubMedPubMed CentralGoogle Scholar
  22. Qian DY, Yan GB, Bai B, Chen Y, Zhang SJ, Yao YC, Xia H. Differential circRNA expression profiles during the BMP2-induced osteogenic differentiation of MC3T3-E1 cells. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie. 2017;90:492–9.View ArticleGoogle Scholar
  23. Zheng Y, Li X, Huang Y, Jia L, Li W. The circular RNA landscape of periodontal ligament stem cells during Osteogenesis. J Periodontol. 2017;88(9):906–14.View ArticlePubMedGoogle Scholar
  24. Wei FL, Wang JH, Ding G, Yang SY, Li Y, YJ H, Wang SL. Mechanical force-induced specific MicroRNA expression in human periodontal ligament stem cells. Cells Tissues Organs. 2014;199(5–6):353–63.View ArticlePubMedGoogle Scholar
  25. Borodina T, Adjaye J, Sultan M. A strand-specific library preparation protocol for RNA sequencing. Methods Enzymol. 2011;500:79–98.View ArticlePubMedGoogle Scholar
  26. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.View ArticlePubMedPubMed CentralGoogle Scholar
  27. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.View ArticlePubMedPubMed CentralGoogle Scholar
  28. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.View ArticlePubMedPubMed CentralGoogle Scholar
  29. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.View ArticlePubMedPubMed CentralGoogle Scholar
  30. Zhang XO, Dong R, Zhang Y, Zhang JL, Luo Z, Zhang J, Chen LL, Yang L. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome Res. 2016;26(9):1277–87.View ArticlePubMedPubMed CentralGoogle Scholar
  31. Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, Haag JD, Gould MN, Stewart RM, Kendziorski C. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013;29(8):1035–43.View ArticlePubMedPubMed CentralGoogle Scholar
  32. miRanda: http://www.microrna.org/microrna/home.do Assessed August 2010.
  33. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.View ArticlePubMedGoogle Scholar
  34. Ala U, Karreth FA, Bosia C, Pagnani A, Taulli R, Leopold V, Tay Y, Provero P, Zecchina R, Pandolfi PP. Integrated transcriptional and competitive endogenous RNA networks are cross-regulated in permissive molecular environments. Proc Natl Acad Sci U S A. 2013;110(18):7154–9.View ArticlePubMedPubMed CentralGoogle Scholar
  35. Qi X, Zhang D-H, Wu N, Xiao J-H, Wang X, Ma W. ceRNA in cancer: possible functions and clinical implications. J Med Genet. 2015;52(10):710–8.View ArticlePubMedGoogle Scholar
  36. Li LJ, Zhao W, Tao SS, Leng RX, Fan YG, Pan HF, Ye DQ. Competitive endogenous RNA network: potential implication for systemic lupus erythematosus. Expert Opin Ther Targets. 2017;21(6):639–48.View ArticlePubMedGoogle Scholar
  37. Li X, Ao J, Wu J. Systematic identification and comparison of expressed profiles of lncRNAs and circRNAs with associated co-expression and ceRNA networks in mouse germline stem cells. Oncotarget. 2017;8(16):26573–90.PubMedPubMed CentralGoogle Scholar
  38. Wang L, Wu F, Song Y, Li X, Wu Q, Duan Y, Jin Z. Long noncoding RNA related to periodontitis interacts with miR-182 to upregulate osteogenic differentiation in periodontal mesenchymal stem cells of periodontitis patients. Cell Death Dis. 2016;7(8):e2327.View ArticlePubMedPubMed CentralGoogle Scholar
  39. Tay Y, Rinn J, Pandolfi PP. The multilayered complexity of ceRNA crosstalk and competition. Nature. 2014;505(7483):344–52.View ArticlePubMedPubMed CentralGoogle Scholar
  40. Huang Y, Zheng Y, Jia L, Li W. Long noncoding RNA H19 promotes Osteoblast differentiation via TGF-beta1/Smad3/HDAC signaling pathway by deriving miR-675. Stem Cells. 2015;33(12):3481–92.View ArticlePubMedGoogle Scholar
  41. Ishibashi O, Ikegame M, Takizawa F, Yoshizawa T, Moksed MA, Iizawa F, Mera H, Matsuda A, Kawashima H: Endoglin is involved in BMP-2-induced osteogenic differentiation of periodontal ligament cells through a pathway independent of Smad-1/5/8 phosphorylation. J Cell Physiol 2010, 222(2):465-473.Google Scholar
  42. Osathanon T, Ritprajak P, Nowwarote N, Manokawinchoke J, Giachelli C, Pavasant P. Surface-bound orientated Jagged-1 enhances osteogenic differentiation of human periodontal ligament-derived mesenchymal stem cells. J Biomed Mater Res A. 2013;101(2):358–67.View ArticlePubMedGoogle Scholar
  43. Li Y, Li SQ, Gao YM, Li J, Zhang B. Crucial role of notch signaling in osteogenic differentiation of periodontal ligament stem cells in osteoporotic rats. Cell Biol Int. 2014;38(6):729–36.View ArticlePubMedGoogle Scholar
  44. Heo JS, Lee SY, Lee JC. Wnt/beta-catenin signaling enhances osteoblastogenic differentiation from human periodontal ligament fibroblasts. Molecules and cells. 2010;30(5):449–54.View ArticlePubMedGoogle Scholar
  45. Thouverey C, Caverzasio J. Focus on the p38 MAPK signaling pathway in bone development and maintenance. BoneKEy reports. 2015;4:711.PubMedPubMed CentralGoogle Scholar
  46. Tang M, Peng Z, Mai Z, Chen L, Mao Q, Chen Z, Chen Q, Liu L, Wang Y, Ai H. Fluid shear stress stimulates osteogenic differentiation of human periodontal ligament cells via the extracellular signal-regulated kinase 1/2 and p38 mitogen-activated protein kinase signaling pathways. J Periodontol. 2014;85(12):1806–13.View ArticlePubMedGoogle Scholar
  47. Li L, Han MX, Li S, Xu Y, Wang L. Hypoxia regulates the proliferation and osteogenic differentiation of human periodontal ligament cells under cyclic tensile stress via mitogen-activated protein kinase pathways. J Periodontol. 2014;85(3):498–508.View ArticlePubMedGoogle Scholar
  48. Chen L, Holmstrom K, Qiu W, Ditzel N, Shi K, Hokland L, Kassem M. MicroRNA-34a inhibits osteoblast differentiation and in vivo bone formation of human stromal stem cells. Stem Cells. 2014;32(4):902–12.View ArticlePubMedGoogle Scholar
  49. Huszar JM, Payne CJ. MIR146A inhibits JMJD3 expression and osteogenic differentiation in human mesenchymal stem cells. FEBS Lett. 2014;588(9):1850–6.View ArticlePubMedPubMed CentralGoogle Scholar
  50. Ounzain S, Micheletti R, Arnan C, Plaisance I, Cecchi D, Schroen B, Reverter F, Alexanian M, Gonzales C, Ng SY, et al. CARMEN, a human super enhancer-associated long noncoding RNA controlling cardiac specification, differentiation and homeostasis. J Mol Cell Cardiol. 2015;89(Pt A):98–112.View ArticlePubMedGoogle Scholar
  51. Zhu M, Xu Y, Chen Y, Yan F, Circular BANP. An upregulated circular RNA that modulates cell proliferation in colorectal cancer. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie. 2017;88:138–44.View ArticleGoogle Scholar
  52. Guo W, Zhang J, Zhang D, Cao S, Li G, Zhang S, Wang Z, Wen P, Yang H, Shi X, et al. Polymorphisms and expression pattern of circular RNA circ-ITCH contributes to the carcinogenesis of hepatocellular carcinoma. Oncotarget. 2017;8(29):48169–77.PubMedPubMed CentralGoogle Scholar

Copyright

© The Author(s). 2017

Advertisement