Human periodontal ligament stem cells (hPDLSCs) demonstrate multilineage differentiation capacity both in vivo and in vitro, forming cementum/PDL-like tissues in vivo and differentiating into osteoblasts, chondrocytes, and adipocytes in vitro (Seo et al., 2004; Gay et al., 2007). To gain insight into the multilineage differentiation capacity of hPDLSCs, we analyzed a publicly available single-cell RNA-sequencing (scRNA-seq) dataset from pooled periodontal ligament tissues of third molar obtained from 3 healthy individuals (aged 21–27 years, 1 female and 2 males) (Yang et al., 2024). Unsupervised clustering using Seurat (Hao et al., 2024) identified nine distinct cell subpopulations (Figure 1A, clusters 0–8), each exhibiting unique gene expression patterns (Figure 1B) that revealed cellular heterogeneity within hPDLSCs.

Notably, cluster 4 was characterized by expression of IGFBP5 (Figure 1B, red rectangles), a stem cell marker for dental mesenchymal progenitors that contribute specifically to periodontal cell lineages (Seidel et al., 2017). Consistent with this finding, unsupervised inference of developmental directions using VECTOR (Zhang et al., 2020) revealed that all subpopulations derived from IGFBP5-expressing cluster 4 cells (Figure 1C, blue rectangle) differentiating along two major developmental directions. Given the observed developmental heterogeneity within hPDLSCs, we then investigated the expression patterns of key regulatory factors governing lineage specification. RUNX2 is the master transcription factor required for osteogenic lineage determination (Komori, 2010), while SCX serves as the key transcription factor for tenogenic differentiation and is highly specific for tendon and ligament tissues (Schweitzer et al., 2001). Since osteogenic and tenogenic differentiation represent critical developmental processes for periodontal tissue regeneration, we analyzed the expression of these lineage-specific transcription factors alongside BMP4, which is widely used to induce in vitro differentiation of human mesenchymal stem cells (Setiawan et al., 2022). Kernel gene-weighted density estimation revealed that within IGFBP5-expressing cluster 4 cells, higher BMP4 expression correlated with higher RUNX2 expression, while lower BMP4 expression correlated with higher SCX expression (Figure 1D), consistent with the established osteogenic function of BMP4 (Beederman et al., 2013). Interestingly, along the two major developmental directions, BMP4 upregulation correlated with decreased RUNX2 and increased SCX expression (Figure 1D, red arrows), implicating a BMP4 dose-dependent effect on hPDLSC lineage specification. Quantitative analysis confirmed that BMP4 upregulation inversely correlated with RUNX2 (Figure 1E, cluster 4→2→0) and positively correlated with SCX expression (Figure 1E, cluster 4→1→3). These results uncover previously unrecognized levels of cellular heterogeneity and developmental plasticity within hPDLSCs and indicate that BMP4 regulates the bias between osteogenic and tenogenic differentiation programs in a dose-dependent manner.

To further dissect the temporal dynamics of BMP4, RUNX2 and SCX during hPDLSC lineage specification, we reconstructed the differentiation trajectory using Monocle2 (Qiu et al., 2017). Trajectory analysis identified nine distinct cellular states (Figure 2A, State 1–9) that span four branch points (Figure 2B, number 1–4) over pseudotemporal differentiation of hPDLSCs. Originating from initial state 8, which expresses IGFBP5 (Figure 2D), hPDLSCs bifurcate into five terminal states (Figure 2C, States 7, 5, 3, 9, and 1). Analysis of transition relationships among cell subpopulations (clusters 0–8) confirmed that all subpopulations are derived from cluster 4 cells (Supplementary Figure S1A,C). Cluster 4 cells are defined by IGFBP5 expression (Supplementary Figure S1B), which aligns with the VECTOR inference results (Figure 1C). Notably, BMP4 was expressed throughout the entire hPDLSC differentiation trajectory (Supplementary Figure S1B). Characterization of BMP4, RUNX2 and SCX expression during hPDLSC differentiation revealed a consistent pattern: at the onset of each cellular state transition, BMP4 expression was low, coinciding with high expression of both RUNX2 and SCX (Figure 2E). During subsequent state transitions, BMP4 gradually increased then decreased at all branch points except the bifurcation from state 8 to state 5 (Figure 2E, State 8→5), where BMP4 maintained high expression (Figure 2E, branch point 3). Notably, following upregulation of BMP4 expression, RUNX2 expression immediately decreased at all branch points (Figure 2E). In contrast, SCX expression immediately increased and remained elevated longer than RUNX2 at all branch points except the bifurcation from state 2 to state 1 (Figure 2E, State 2→1). Interestingly, SCX expression peaked at two specific bifurcations: from state 8 to state 7 (Figure 2E, State 8→7) and from state 2 to state 9 (Figure 2E, State 2→9). These peaks occurred after BMP4 expression had declined. Visualization of BMP4, RUNX2, and SCX expression kinetics during hPDLSC cellular state transitions (Figure 2F) further confirmed these observations. Collectively, these results reveal biphasic expression dynamics of BMP4, RUNX2, and SCX expression during hPDLSC differentiation: low BMP4 expression correlates with high co-expression of RUNX2 and SCX, while BMP4 upregulation exhibited an inverse correlation with RUNX2 expression and a positive, dynamic relationship with SCX expression. These findings further support a dose-dependent regulatory mechanism of BMP4 in hPDLSC lineage specification.

Building on our findings (Figures 1, 2) and previous studies showing that BMPs at various concentrations differently modulate the proliferation and osteogenic differentiation of hPDLSCs (Hakki et al., 2014), we further characterized the osteogenic and tenogenic differentiation potential of hPDLSCs following treatment with low (10 ng/mL) or high (100 ng/mL) BMP4 concentrations. We selected these specific concentrations based on established literature demonstrating that 10 ng/mL represents a threshold for maintaining PDLSC multipotency (Liu et al., 2013), while 100 ng/mL induces robust differentiation responses (Hyun et al., 2017). These concentrations at opposite ends of the physiologically relevant range allowed us to directly test whether BMP4 exerts dose-dependent effects consistent with the biphasic expression patterns observed in our scRNA-seq analysis. Notably, low-concentration BMP4 treatment increased expression of both SCX and RUNX2 (Figure 3A), consistent with our pseudotemporal analysis (Figure 2E). To confirm that BMP4 directly mediates upregulation of these transcription factors, we treated hPDLSCs with Noggin, a well-known BMP4 antagonist. This intervention effectively attenuated the BMP4-induced upregulation of both SCX and RUNX2 (Figure 3A). A key observation emerged when comparing the effects of low- and high-concentration BMP4: relative to the low-dose group, high concentration BMP4 further enhanced SCX expression while simultaneously suppressing RUNX2 expression, which aligns with our scRNA-seq analysis (Figures 1, 2). Given the positive correlation between SCX expression and both exogenous BMP4 treatment (Figure 3B) and endogenous BMP4 expression levels (Figures 1, 2) in hPDLSCs, we next investigated whether SCX responds directly to elevated BMP4. Using a lentiviral overexpression system to upregulate BMP4 in hPDLSCs, we found that BMP4 overexpression increased nuclear accumulation of SCX (Figures 3F–H') by approximately 2.5-fold (Figure 3I) compared with the control group (Figures 3C–E'), directly confirming BMP4’s pivotal role in driving hPDLSCs commitment to the tenogenic lineage. Collectively, these data validate that BMP4 elicits dose-dependent effects on hPDLSC lineage specification: low BMP4 concentration preserves hPDLSC multipotency, whereas high BMP4 concentration promotes tenogenic lineage commitment while attenuating osteogenic differentiation.

Given that BMP4 initiates signal transduction through ligand-receptor interactions, we utilized CellChat (Jin et al., 2021) to infer and analyze the cell-cell communications among hPDLSCs subpopulations (cluster 0–8). This analysis aimed to understand how in vitro treatment with high-concentration BMP4 elicits opposing osteogenic and tenogenic differentiation responses (Figure 3). Quantitative assessment of incoming BMP signaling in osteogenic subpopulations (Figure 4A, cluster 0 and 2) and tenogenic subpopulations (Figure 4B, cluster 1 and 3) revealed that BMP4 signaling for these two lineages is transduced through different receptor pairs. In osteogenic hPDLSCs, BMP4 signals are transduced through pairs of BMPR1A or BMPR1B with one of three type II receptors: BMPR2, ACVR2A, and ACVR2B (Figure 4A).

These three type II receptors are functionally similar, as they interact with an overlapping group of TGF-β family growth factors (Chu et al., 2022). In contrast, tenogenic hPDLSCs utilize BMPR1B paired with the same three type II receptors (BMPR2, ACVR2A, and ACVR2B) for BMP4 signaling transduction (Figure 4B). Based on these findings, we propose a mechanistic model underlying the dose-dependent osteogenic and tenogenic differentiation of hPDLSCs in response to in vitro BMP4 treatment. Specifically, two key differences between osteogenic (RUNX2+) and tenogenic (SCX+) hPDLSCs collectively drive their divergent differentiation responses: (1) distinct cell surface BMP receptor pairs, and (2) inherently different patterns of RUNX2 and SCX expression changes upon BMP4 upregulation (Figure 4C).

Publicly available hPDLSC scRNA-seq data (Yang et al., 2024) from third molars of three healthy donors (aged 21–27 were processed and subjected to unsupervised clustering with the standard Seurat (Hao et al., 2024) pipeline and parameters suggested by VECTOR software (Zhang et al., 2020). Since these hPDLSC were obtained from three individual donors, discontinuous clusters from unsupervised clustering were excluded from downstream analyses. Uniform Manifold Approximation and Projection (UMAP) and principal component analysis (PCA) embeddings from the resultant Seurat object were then used by VECTOR to infer developmental directions. To reconstruct hPDLSC differentiation trajectories, cell embeddings and clusters from the resultant Seurat object were converted to a CellDataSet object and ordered in pseudotime using Monocle2 (Qiu et al., 2017) with default parameters. For analysis of cell-cell communication between hPDLSC subpopulations, we employed CellChat (Jin et al., 2021). Gene expression data and cell group information from the resultant Seurat object were used to create a CellChat object for computing communication probability between hPDLSC subpopulations. Incoming BMP signaling to osteogenic/RUNX2+ and tenogenic/SCX + cells was then visualized using CellChat’s internal function netVisual_bubble. Pseudobulk analysis between SCX+ and RUNX2+ hPDLSC cells were performed on raw counts extracted from the processed Seurat object using DESeq2 (Love et al., 2014). Differentially expressed genes (DEGs) were visualized using the EnhancedVolcano package (Blighe et al., 2018). DEGs with p-value <0.05 were considered statistically significant and subjected to functional classification using clusterProfiler (Yu et al., 2012).

Human periodontal ligament stem cells (hPDLSCs) were isolated according to a previous study (Schweitzer et al., 2001). Normal third molars (n = 10) were collected at School and Hospital of Stomatology, Fujian Medical University. Briefly, periodontal ligament cells were scraped from third molars and digested with 3 mg/mL collagenase type I (Worthington Biochem) and 4 mg/mL dispase (Roche) for 1 h at 37 °C. Isolated hPDLSCs were then cultured in Dulbecco’s modified Eagle’s medium supplemented with 10% fetal bovine serum, 200 U/mL penicillin, 200 μg/mL streptomycin, and 100X GlutaMAX™ Supplement (all from Gibco). Cells were grown at 37 °C in a humidified atmosphere containing 5% CO₂. Passages 3-5 of hPDLSCs were used for experiments. After reaching 80%–90% confluence, cells were incubated with BMP4 (R&D Systems), NOG (R&D Systems), or control and BMP4 lentivirus (pNL-BMP4-IRES2-EGFP) at indicated concentrations for specified time periods in complete medium. The medium was changed every other day.

Total RNA was isolated using the Simply P Total RNA Extraction Kit (BioFlux) according to the manufacturer’s instructions. Complementary DNA was reverse-transcribed from total RNA using the PrimeScript RT Reagent Kit (Takara). RT-qPCR was performed in triplicate using SYBR Premix Ex Taq™ (Takara) on a Thermal Cycler Dice™ Real Time System TP800 (Takara). Five independent experiments were performed, and mean expression levels are shown. The data were analyzed with the Thermal Cycler Dice™ Real Time System analysis software (Takara). Data were analyzed using the Thermal Cycler Dice™ Real Time System analysis software (Takara). Gene expression levels were normalized to ACTB and expressed as fold changes relative to control. Primer sequences for all genes are listed in Supplementary Table S1.

For immunofluorescence analysis, hPDLSCs were sub-cultured on glass coverslips in 24-well cell culture dishes. After lentivirus treatment for 96 h, cells were washed with 1× PBS and fixed in 4% paraformaldehyde (PFA) for 10 min at room temperature. After removing PFA and washing twice with PBS containing 0.1% Tween-20, cells were blocked and incubated with SCX primary antibody (Abcam, ab58655, 1:200) for 2 h, followed by incubation with Alexa Fluor-488-conjugated secondary antibody. DAPI (Invitrogen) staining was used to label nuclei according to the manufacturer’s instructions. Images were acquired using a Leica DM 5500B microscope equipped with a Leica DFC 495 camera or a Zeiss LSM700 meta-confocal laser-scanning microscope. To quantify cytoplasmic and nuclear SCX fluorescence intensity, we employed CellProfiler (Stirling et al., 2021). Briefly, nuclei were identified using Global thresholding with the Otsu approach. Whole cells and cytoplasm were then segmented using the Identify Secondary and Tertiary Objects modules. Mean SCX intensity was measured for nuclei and cytoplasm, and nucleus-to-cytoplasm ratios were calculated. Data were exported to Excel for analysis in Prism.