Time series clustering of mRNA and lncRNA expression during osteogenic differentiation of periodontal ligament stem cells

Cell isolation and culture

We isolated human PDLSCs using the process described previously (Seo et al., 2004). Briefly, we collected PDL tissue from the middle third root of healthy premolars, cut it into pieces, and digested it in equal volumes of type I collagenase and dispase. Next, we cultured the isolated cells using α-Modified Eagle Medium (α-MEM) supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin (Gibco, Grand Island, NY, USA). We passaged the cells using 0.25% trypsin and further expanded them to passage 4 before they were used in the experiments. For osteogenic induction, we cultured PDLSCs in α-MEM supplemented with 10% FBS, 100 nM dexamethasone, 200 μM L-ascorbic acid, and 10 mM β-glycerophosphate. We changed the medium every 2 days and harvested the cells at the indicated time points. The study protocol was approved by the Ethics Committee of Peking University School of Stomatology (PKUSSIRB-2011007).

RNA sequencing

We performed sequencing of the mRNA and lncRNA as described previously (Zheng et al., 2016; Zheng & Jia, 2016). Briefly, we extracted the RNA of the PDLSCs using the Trizol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA). After digestion with DNase, we depleted rRNA using a Ribo-Zero magnetic kit and prepared a library according to the manufacturer’s instructions (Epicentre Biotechnologies, Madison, WI, USA). We performed paired-end sequencing using the Illumina Hiseq X Ten. After removing low-quality reads, we trimmed the raw reads and mapped them to the human genome (hg19) using TopHat2. We used HTseq to count the genes and calculated the reads per kilobase transcriptome per million mapped reads (RPKM) to evaluate the gene expression level. We used the EB-Seq package for the differential gene expression analysis among PDLSCs at each time point. We defined the differentially expressed genes (DEGs) based on fold changes greater than or equal to 2.0 and a false discovery rate of less than 0.05. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the Database for Annotation Visualization and Integrated Discovery. Raw data of the performed RNA-seq have been recorded in the GEO public database (accession number: GSE99958).

STEM analysis and Gene Set Enrichment Analysis

We used the STEM (Ernst, Nau & Bar-Joseph, 2005) clustering algorithm to identify temporal gene expression profiles during the osteogenic differentiation of PDLSCs. Representative temporal expression profiles were first defined as model profiles, which were independent from the sequencing data. The values of the gene expression were transformed to log ratios relative to the expression at Day 0 (D0). Then each gene was assigned to the filtering criteria of the model profiles, and the correlation coefficient was determined. We performed standard hypothesis testing using the true ordering of time points, and determined the p-value using the number of genes assigned to the model profile and the expected number of assigned genes (adjusted p-value, 0.05 by Bonferroni correction). The boxes of figures were colored if the profiles were statistically significant. We also performed the Gene Set Enrichment Analysis (GSEA) to identify the GO biological processes involving genes with significant profiles as previously described (Subramanian et al., 2007).

Alkaline phosphatase staining

We undertook a histochemical analysis of alkaline phosphatase (ALP) activity as previously described (Huang et al., 2015). For the ALP staining, we seeded cells in 24-well plates and cultured them in an osteogenic medium for 3, 7, or 14 days. Subsequently, we fixed the cells in 4% paraformaldehyde for 30 min. After washing them with PBS, we incubated the cells for 20 min at 37 °C in an alkaline solution according to the directions for the NBT/BCIP staining kit (CoWin Biotech, Beijing, China). We took images of PDLSCs using a Zeiss system.

Mineralization assay

To visualize the mineral deposition in the extracellular matrix (ECM), we cultured the PDLSCs in an osteogenic medium for 3, 7, or 14 days. We then fixed the cells in 4% paraformaldehyde and stained them with 0.1% Alizarin Red S (ARS; Sigma-Aldrich, St. Louis, MO, USA) for 20 min at room temperature. We took images of this samples using a Zeiss system.

Quantitative realtime PCR

We extracted total RNA using the Trizol reagent and carried out quantitative real-time PCR (qRT-PCR) in triplicate using the Power SYBR Green® PCR Mastermix (Applied Biosystems, Foster City, CA, USA) as previously described (Zheng et al., 2016). The mean expression values were calculated relative to the mean expression level of the housekeeping gene GAPDH. The primers that we used are listed in Table 1. We conducted a statistical analysis using one-way ANOVA, followed by Student–Newman–Keuls post hoc tests using SPSS 16.0 software.

siRNA knock down assay

Periodontal ligament stem cells were transfected with siRNA to knock down the expression of MEG8 and MIR22HG. Cells in six-well plates were transfected with 10 µl siRNA, 250 µl Opti-MEM medium, and three µl Lipofectamine® RNAiMAX (Life Technologies, Carlsbad, CA, USA). Total RNA was extracted 24 h after transfection and the expression levels of runt-related transcription factor 2 (RUNX2) and osteocalcin (OCN) were examined.

Summary of the RNA-seq data

The human PDLSCs showed osteogenic differentiation potential when they were cultured in the osteogenic medium. The calcified nodules and ARS staining of the PDLSCs were significantly increased on Days 7 and 14 (D7 and D14, respectively; Fig. 1A). Using the RNA-seq technique, we determined the mRNA and lncRNA expression profiles of the PDLSCs during osteogenic differentiation, which provided insight into the molecular mechanisms underlying osteoblast differentiation. We generated paired-end libraries and sequenced more than 2 × 10⁸ reads per sample. We filtered the raw reads (Fig. S1) and mapped them to the hg19 version of the human genome. Our analysis of the gene structure, as well as their distribution on the chromosomes of the mapped reads, is shown in Figs. 1B and 1C. Generally, the mapping numbers per chromosome correlated broadly with chromosome length (Fig. 1B). Chromosome Y received a relatively low number of read mappings. The majority of reads were localized in the exon and intron regions. Moreover, the intergenic regions produced a great number of transcripts (Fig. 1C). Pairwise correlation analyses using the gene-level expression data demonstrated that there was a high degree of similarity between the PDLSCs on Days 3 (D3) and D7 (Fig. 1D).

Characterization of the mRNA transcriptome of PDLSCs

Our analysis of the RNA-seq data for Gencode-annotated mRNAs showed that 24,881 mRNAs were detected. To gain insight into the roles of these mRNAs, we performed comparative expression analyses according to the differentiation stages of the osteoblasts (D3 vs. D0, D7 vs. D0, and D14 vs. D0). As shown in Fig. 2A, 664 mRNAs were upregulated and 477 mRNAs were downregulated on D3. On D7, 969 mRNAs were upregulated and 485 mRNAs were downregulated. On D14, 1,046 mRNAs were upregulated and 732 mRNAs were downregulated. A Venn analysis showed that 290 genes were simultaneously upregulated on D3, D7, and D14 (Fig. 2B). These genes were enriched in bone remodeling-associated processes, including ECM organization and focal adhesion. Genes associated with cell differentiation and cell apoptosis were also upregulated (Fig. 2C). In contrast, 194 genes were downregulated on D3, D7, and D14 (Fig. 2D). These genes were enriched mainly in processes related to the regulation of cellular proliferation, including mitotic nuclear division and sister chromatid cohesion (Fig. 2E).

Characterization of lncRNA transcriptomes

We organized the lncRNAs according to their expression profiles. The overall lncRNA expression profiles are presented in Fig. 3A. Small distinct groups of lncRNAs were expressed either early or late in osteogenic differentiation or predominately at one specific time point (e.g., D3 or D7). A much larger group was expressed mainly on D0 and D14. A pairwise comparison of lncRNA expression patterns during the osteogenic differentiation revealed that, from D0 to D3, 75 genes were upregulated and 148 genes were downregulated (Fig. 3B). Compared with the PDLSCs on D0, 145 genes were upregulated and 110 genes were downregulated on D7 (Fig. 3C), while 163 upregulated genes and 208 downregulated genes were detected on D14 (Fig. 3D). A Venn diagram showed that 17 genes were simultaneously upregulated and 31 lncRNAs were downregulated on D3, D7, and D14 (Figs. 3E and 3F).

Temporal gene expression patterns of mRNAs

We normalized the sequencing data to D0 and identified the temporal gene expression profiles using STEM. Within the 25 model profiles, 14 temporal mRNA profiles were statistically significant (Fig. 4A). The black lines in the profile boxes depict the gene expression patterns over the four time points. The profile number on the top left corner of each profile box was assigned by STEM, the number on the bottom left represents the adjusted p-value, and the number on the top right corner represents the cardinality of each cluster. We found continuous upregulation and downregulation patterns in profiles 0 (Fig. 5A) and 25. Gradual increases or decreases were distributed across profiles 1 (Fig. 5B), 9, 12 (Fig. 5C), 13, 16, 21, 22, and 24. Meanwhile, profiles 5, 17, 20, and 23 (Fig. 5D) showed biphasic responding expression patterns. Genes in profiles 0, 1, 12, and 23 were summarized in File S1. We then analyzed the genes associated with these profiles using GSEA analysis and retrieved the enriched GO terms.

Among the 14 profiles, the genes in profile 0 were associated with epigenetic regulation, including DNA methylation and histone modification (Fig. 5E). Profile 1 genes were enriched in cell mitosis (Fig. 5F). Profile 12 genes were associated with cell differentiation and osteoblast differentiation (Fig. 5G), while profile 23 genes were enriched in the collagen catabolic process, odontogenesis, and skeletal system development (Fig. 5H; p-value = 0.003). We found that 17 out of 112 genes in profile 12 were annotated as cell differentiation and osteoblast differentiation, and eight out of 118 genes in profile 23 were annotated as odontogenesis and skeletal system development. Specifically, COL6A1, VCAN, RRBP1, and CREB3L1 were annotated as osteoblast differentiation. Genes annotated as odontogenesis and skeletal system development included AMELX, BMP2, LAMA5, SCN5A, COL3A1, BMP1, ANKH, and COL12A1.

Furthermore, we used a KEGG pathway analysis to explore the signaling pathways of DEGs involved in profiles 0, 1, 12, and 23. Interestingly, the genes in profile 0 and profile 1 were enriched in signaling associated with cancer development (Figs. 6A and 6B), which is characterized by the misregulation of cell proliferation and apoptosis. Genes in the profile 12 and 23 were not only enriched in ECM-receptor interactions, but also in the Notch and PI3K-AKT signaling pathways (Figs. 6C and 6D).

Temporal gene expression patterns of lncRNAs

We classified the differentially expressed lncRNAs into 25 profiles according to the STEM analysis results and ranked the profiles according to their significance (p-value). The p-values of 8 clusters were statistically significant (Fig. 4B). The lncRNAs in clusters 0 and 25 represented lncRNAs whose expression was continuously increased or decreased. The lncRNAs in clusters 12 and 13 revealed differential expression at D7, while lncRNAs in clusters 16 and 13 revealed differential expression at both D7 and D14. lncRNAs in clusters 5, 17, 20, and 23 corresponded to lncRNAs with biphasic expression patterns. Since the differentially expressed mRNAs in gene expression profile 12 were significantly enriched in osteoblast differentiation, the lncRNAs included in profile 12 were potentially related to osteogenic differentiation of PDLSCs.

Validation of the expression of mRNAs and lncRNAs that associated with osteoblast differentiation

We validated the expression of osteogenic mRNAs and lncRNAs in the profile 12 using qRT-PCR. Though most lncRNAs were expressed at a low level, some lncRNAs were highly expressed, rivaling the expression of many mRNAs. Therefore, we applied an expression cutoff to remove the transcripts that were expressed at a low level. We found that two lncRNA genes in the profile 12 were expressed at a RPKM of more than 10. The expression patterns of mRNAs and lncRNAs showed similar trends. The levels of COL6A1 and RRBP1 were increased at D3 and stayed relatively stable afterward (Figs. 7A and 7C). The expression of VCAN was decreased significantly during the osteoblast differentiation of PDLSCs (Fig. 7B). CREB3L1 increased significantly at D3, but decreased from D7 to D14 (Fig. 7D). For the lncRNAs, the expression of MEG8 was significantly increased at D7, but reduced at D14 (Fig. 7E). Overall, the expression of MIR22HG was upregulated in the osteogenic medium, and its level peaked at D7 (Fig. 7F). Moreover, the expression of RUNX2 and OCN decreased once the siRNA targeting MEG8 and MIR22HG were applied to PDLSCs, indicating an osteogenic role for the two lncRNAs (Figs. 7G and 7H).