Single-cell RNA-seq analysis of rat molars reveals cell identity and driver genes associated with dental mesenchymal cell differentiation

Abstract

Background

The molecular mechanisms and signaling pathways involved in tooth morphogenesis have been the research focus in the fields of tooth and bone development. However, the cell population in molars at the late bell stage and the mechanisms of hard tissue formation and mineralization remain limited knowledge.

Results

Here, we used the rat mandibular first and second molars as models to perform single-cell RNA sequencing (scRNA-seq) analysis to investigate cell identity and driver genes related to dental mesenchymal cell differentiation during the late bell hard tissue formation stage. We identified seven main cell types and investigated the heterogeneity of mesenchymal cells. Subsequently, we identified novel cell marker genes, including Pclo in dental follicle cells, Wnt10a in pre-odontoblasts, Fst and Igfbp2 in periodontal ligament cells, and validated the expression of Igfbp3 in the apical pulp. The dynamic model revealed three differentiation trajectories within mesenchymal cells, originating from two types of dental follicle cells and apical pulp cells. Apical pulp cell differentiation is associated with the genes Ptn and Satb2, while dental follicle cell differentiation is associated with the genes Tnc, Vim, Slc26a7, and Fgfr1. Cluster-specific regulons were analyzed by pySCENIC. In addition, the odontogenic function of driver gene TNC was verified in the odontoblastic differentiation of human dental pulp stem cells. The expression of osteoclast differentiation factors was found to be increased in macrophages of the mandibular first molar.

Conclusions

Our results revealed the cell heterogeneity of molars in the late bell stage and identified driver genes associated with dental mesenchymal cell differentiation. These findings provide potential targets for diagnosing dental hard tissue diseases and tooth regeneration.

Background

Tooth organogenesis is an extremely complex biological process. It depends on a series of reciprocal interactions between the epithelium and neural crest-derived mesenchyme, which develop into the enamel organ, dental papilla, and dental follicle, respectively [1]. Besides the interaction between epithelium and mesenchyme, tooth development involves cell differentiation, morphogenesis, tissue mineralization, maturation, tooth eruption, and integration with its surrounding tissues [2]. At the late bell stage, dental papilla cells differentiate into dentin-forming odontoblasts, and inner enamel epithelia differentiate into enamel-forming ameloblasts. These processes are regulated by a complex network of cell signaling pathways. A series of signaling molecules and their receptors play an important role in tooth formation, including wingless-type (WNT) pathway [3], sonic hedgehog (SHH) pathway [4], fibroblast growth factors (FGF) pathway [5], bone morphogenesis proteins (BMPs) family [6], ectodysplasin A (EDA) pathway [7], and transforming growth factor-beta (TGF-β) pathway [8, 9]. In addition to these signals, mutations in some genes regulated by these pathways have been shown to cause dental defects [10].

To date, previous studies have provided some information about tooth organogenesis, genetic mechanisms, and signaling pathways involved in dental diseases [1, 10, 11]. More than 300 genes have been reported to be associated with odontogenesis, but most are involved in tooth germs at relatively early morphogenetic stages, from initiation through the cap/early bell stage [12]. Analyses of molecular signaling networks and new insights into cellular heterogeneity have greatly improved our knowledge of tooth development and homeostasis. However, these studies primarily focus on the regulation of organ, tissue, or cell population levels.

Recently, the technology of single-cell RNA sequencing (scRNA-seq) can perform unbiased transcriptional profiling at the single-cell level, reveal complex and diverse cell populations, and delineate the trajectories of cell lineages during development [13]. Studies on continuously growing mouse incisors with scRNA-seq have revealed the complexity of cellular composition and provided unprecedented cell type annotations in mammalian teeth at various stages [14]. The transcriptomic characteristics of enamel-forming dental epithelial cells in whole mouse incisors have been identified [15]. PRX1+ cells have been shown to be involved in the development of mouse first molar and the angiogenesis in periodontal ligament cell development and repair [16]. The incisor has long been used as a model for continuous tooth development. However, there is a lack of information about the cell population in molars during the late bell stage, as well as the molecules that regulate the functional differentiation of cells and the mineralization of dental hard tissues.

Therefore, we performed scRNA-seq transcriptome analysis on rat mandibular molars at postnatal days 5 (PN5) to investigate the cellular composition and gene expression of different cell populations in the tooth germs at the late bell stage of hard tissue formation and to identify driver genes associated with mesenchymal cell differentiation.

Results

Characterization of rat mandibular molars single-cell atlas

To explore the cell identity and genes related to the cell fate of molars during the late bell stage, the mandibular first and second molars (M1 and M2) of rat at PN5 were dissected for enzymatic digestion and subjected to scRNA-seq (Additional file 1: Fig. S1). A total number of 18,998 cells and 21,714 cells were obtained for M1 and M2, respectively. After quality control filtering out low-quality and doublet cells, 13,007 and 14,976 cells for M1 and M2, respectively, were retained for subsequent analysis (Additional file 1: Table S1). The mean and median numbers of detected genes per cell were 2090 and 2123 for M1 and 2018 and 2054 for M2, respectively.

Using the 2000 most variable genes, Uniform Manifold Approximation and Projection (UMAP) dimensional reduction identified 20 cell clusters (Fig. 1a). We fitted the data for M1 and M2 and indicated that M1 and M2 cell distributions were largely consistent, but the proportion of cells in each cluster are different (Additional file 1: Fig. S2-4). The clusters were then annotated into seven main cell types based on the classic cell markers, including Sox9+Msx2+ mesenchymal cells, Pitx2+Krt14+ dental epithelium, C1qa+Aif1+ macrophages, Napsa+ lymphocytes, Cdh5+ endothelial cells, Rgs5+ perivascular cells, and Cdk1+ cycling cells (Fig. 1a–c). Among them, mesenchymal cells contained multiple clusters. Five differentially expressed genes (DEGs) for each cluster were shown in the heatmap (Fig. 1d). We found that cluster 15 exhibited gene expression patterns similar to macrophages (Additional file 2). Functional enrichment analysis of DEGs in cluster 15 indicated associations with immune regulation, osteoclast differentiation, and negative regulation of ossification (Additional file 3). However, it did not express the osteoclast marker gene Acp5. Therefore, we defined cluster 15 as a type of macrophage that is in the process of differentiating into osteoclasts.

[IMAGE OMITTED: SEE PDF]

Heterogeneity of the mesenchymal compartment in rat molar

Mesenchymal cells in teeth constitute the cementum, dentin, and soft tissue of the dental pulp cavity. Although multiple clusters were divided into mesenchymal cells during cell annotation, heterogeneity, and differential gene expression still existed among each cluster. We further annotated the molar mesenchymal cell population (Fig. 2a). We identified two types of Kit+ dental follicle cells, cluster 0 and 2, which exhibited high expression of osteogenic-related gene Igfbp5 and a novel marker Pclo (Fig. 2b). Functional enrichment analysis indicated that cluster 0 is associated with the regulation of osteoblast differentiation, while cluster 2 is related to sensory organ development and sodium ion transmembrane transport (Additional file 1: Fig. S5). Cluster 5 was annotated as Smpd3+ odontoblasts, and cluster 6 was annotated as pre-odontoblasts by the marker gene Gsc (Fig. 2c, d). Wnt10a was intensively expressed in pre-odontoblasts (Fig. 2d). Cluster 4 was defined as osteoblasts through gene set enrichment analysis. Postn+ periodontal ligament cells were identified in cluster 18. Fst and Igfbp2 were specifically expressed in this cluster, which could be used as markers for the isolation of periodontal ligament cells at this stage (Fig. 2e). Cluster 16 was annotated as Smoc2+/Sfrp2+ apical pulp. Igfbp3 was mainly expressed in this cluster (Fig. 2f). The expression of Igfbp3 in the apical pulp was validated with immunofluorescence (Fig. 2g, h). Clusters 1, 3, 7, and 10 were defined as dental pulp. Enrichment analysis of DEGs in cluster 3 revealed associations with skeletal system development, extracellular matrix organization, and osteoblast differentiation. Cluster 7 was related to extracellular matrix organization, skeletal system development, and the regulation of Igf transport and uptake by Igfbps. Therefore, we defined cluster 7 as the distal pulp involved in dentin formation [17]. Cluster 1 was associated with ECM proteoglycans, collagen fiber organization, and response to growth factor stimulus. The function of cluster 10 was associated with translation, response to stress, and immune response, suggesting it might be a type of dental pulp cell involved in immune regulation. The function of cluster 13 was related to response to growth factors, skeletal system development, osteoclast differentiation, and odontogenesis (Additional file 3). It was defined as cementoblasts through gene set enrichment analysis.

[IMAGE OMITTED: SEE PDF]

Maturation and differentiation trajectory of rat molar cells

RNA velocity was used to construct cell development and differentiation trajectories of rat molar cells and show how gene expression correlates with cell fate to identify putative driver genes. Velocities derived from the dynamical model for PN5 M1 and M2 were visualized as streamlines in a UMAP-based embedding (Fig. 3a, Additional file 1: Fig. S6a). In the dynamic model, we identified three differentiation trajectories within the mesenchyme. The first trajectory involves the differentiation of apical pulp cells into pre-odontoblasts, which subsequently differentiate into odontoblasts. The second trajectory consists of apical pulp cells differentiating into other functional dental pulp cells. The third trajectory includes two parts: one part involves the differentiation of dental follicle cells (cluster 0) into osteoblasts and the other part involves the differentiation of dental follicle cells (cluster 2) into cementoblasts (Fig. 3a). These differentiation trajectories are consistent with our current knowledge of dental mesenchymal cell differentiation, indicating that our dynamic model is reliable. Based on scVelo, genes with pronounced dynamic behavior often contain high likelihood and these genes are defined as driver genes in dynamic model. The dynamical model allowed us to systematically identify putative driver genes with high likelihoods. Gene expression dynamics resolved along latent time revealed a clear cascade of transcription in the top 300 likelihood-ranked genes in two molars (Fig. 3b, Additional file 1: Fig. S6b), 26.3% of the putative driver genes between M1 and M2 overlapping, including Nudt4, Vcan, Tnc, Pak1, Anxa5, Me1, and Col12a1 (Additional file 1: Fig. S6c). The top 300 driver genes (Additional file 4) in the first molar are functionally enriched in tissue morphogenesis, extracellular matrix organization, and skeletal system development (Additional file 1: Fig. S7a). For the second molar, they are enriched in skeletal system development, response to growth factors, and biomineral tissue development (Additional file 1: Fig. S7b).

[IMAGE OMITTED: SEE PDF]

The cluster-specific top five genes of the two dental follicle cell clusters (clusters 0 and 2) both included Kcnt2, Pak1, and Nebl. In cluster 0, the increased expression of Tnc and the decreased expression of Vim would promote the differentiation of dental follicle cells into osteoblasts (Fig. 3c). As for cluster 2, the increased expression of Slc26a7 and Fgfr1 would promote the differentiation of dental follicle cells into cementoblasts (Fig. 3d). The top five cluster-specific genes of the apical pulp showed that Ptn and Satbs are associated with apical pulp cells differentiation (Fig. 3e).

Cell cluster-specific regulation by transcription factors in molars

Using pySCENIC, we identified the cluster-specific regulons for each cluster and determined the enriched target genes corresponding to each transcription factor. We found that the top five regulons differ among different cell types (Additional file 5). In the first molar, the top five regulons in osteoblasts were Klf6, Egr1, Nfil3, Klf4, and Fos, while in odontoblasts, they were Creb5, Tcl4, Snai2, Osr2, and Hivep. In the second molar, the top five regulons in osteoblasts were Egr1, Nr2f2, Nfil3, Crem, and Lhx8, among which two regulons were shared with osteoblasts in the first molar. Furthermore, we focused on the regulons of two types of dental follicle cells (clusters 0 and 2) in the second molar. We found that the top three regulons were the same: Tfap2b, Thrb, and Msx1. These results indicated that the regulons identified are cell cluster-specific and could be used to reflect differences in cell types or cell functions.

Validation of the function of driver gene TNC in odontoblastic differentiation

In the results of RNA velocity analysis, we found that the Tnc is an odontoblasts-specific driver gene in the second molar, but not in the first molar (Additional file 4). We hypothesized that Tnc could function in odontoblastic differentiation during dentin formation. Immunohistochemical analysis showed that the highest expression level of Tnc was at PN5, compared to PN1 and PN10 (Fig. 4). Tnc was mainly expressed in odontoblasts and dental pulp cells (Fig. 4d). The mRNA and protein expression levels of TNC increased in human dental pulp stem cells (hDPSCs) after 7 and 14 days of odontogenic induction (Fig. 5a, b). Knockdown experiments were performed to investigate whether TNC affects the function of hDPSCs. Compared to cells transduced with the shCtrl lentivirus, the expression levels of TNC mRNA and protein were both significantly lower in the shTNC1 and shTNC2 cells, indicating efficient downregulation of TNC by the specific shRNA (Fig. 5c, d). ShTNC2 cells with higher knockdown efficiency were selected for subsequent experiments, and shTNC1 cells were used to validate the results of 14 days of odontogenic induction. To further characterize the potential role of TNC in odontogenic differentiation, the expression of odontoblast differentiation markers was examined. After 7 and 14 days of induction, the mRNA expression of DSPP and RUNX2 were lower in the cells transduced with the shTNC2 lentivirus compared to the cells transduced with the shCtrl lentivirus (Fig. 5e). Furthermore, the protein levels of DSPP and RUNX2 were also decreased by the TNC knockdown (Fig. 5f). Lastly, alizarin red S (ARS) staining revealed that TNC downregulation inhibited the mineral deposition on day 14 of odontogenic induction (Fig. 5g). The results of shTNC1 cells after 14 days of odontogenic induction were consistent with those of shTNC2 cells (Additional file 1: Fig. S8).

[IMAGE OMITTED: SEE PDF]

[IMAGE OMITTED: SEE PDF]

Expression of osteoclast genes is increased in the first molar macrophages

At PN5, the hard tissue development in the rat’s first molar tooth crown extended from the cusp to the neck of the tooth but did not reach the tooth cervix. Meanwhile, the hard tissue of the second molar was forming in the tooth cusp. Hence, we analyzed the DEGs in the two tooth germs to gain further insights into the gene regulation involved in tooth germ development. A study on the human tooth germ from growing third molar revealed that immune cells make up 83% of all tooth germ cells, demonstrating their regulatory responsibilities [18]. In our study, we observed a higher proportion of macrophages in the first molar compared to the second molar (Additional file 1: Fig. S2b). Thus, we focused on the DEGs in macrophages (Fig. 6ab). Notably, the expression of chemokine ligands Ccl7, Ccl2, and Ccl24 was significantly higher in the first molar than in the second molar (Fig. 6c). These chemokines bind to CCR, participate in NF-kB and TGF-B signaling pathway, serving as crucial mediators of osteoclastogenesis [19, 20]. Additionally, the expression of Ifrd1 [21] and Slfn2 [22], which are essential for activating osteoclast differentiation, was also significantly higher in the first molar (Fig. 6c). Conversely, the expression of Stmn1 [23] and Ptn [24], which promote osteoblast differentiation, was increased in the second molar (Fig. 6d). These results suggest that osteoclast differentiation of is enhanced in the first molar, whereas osteoblast differentiation is enhanced in the second molar.

[IMAGE OMITTED: SEE PDF]

Discussion

A tooth includes three hard tissues, enamel, dentin, and cementum. Dental hard tissue hypoplasia diseases cause tooth structure defects or even loss, which affect the masticatory function and maxillofacial esthetics, and psychological status [25, 26]. Genetic factor plays an important role in these diseases. By constructing the gene expression profile of the dental hard tissue formation process, the pathogenesis of diseases could be well elucidated and assisted with clinical diagnosis and treatment [27, 28]. In this study, we performed single-cell transcriptome of the mandibular first and second molars tooth germ of rats at postnatal days 5. We deciphered the cell atlas in two molars and revealed specific genes related to cell fate and associated with the process of dental mesenchymal cell differentiation in the hard tissue formation stage.

Tooth development involves a series of cell proliferation, apoptosis, and differentiation. As tooth development progresses from crown to root formation, cell proliferation decreases in epithelial cell types while it is maintained or increases in mesenchymal cells [29]. In this study, transcriptome analysis divided tooth molar cells into seven main cell types, including Sox9+Msx2+ mesenchymal cells, Pitx2+Krt14+ dental epithelium, C1qa+Aif1+ macrophages, Napsa+ lymphocytes, Cdh5+ endothelial cells, Rgs5+ perivascular cells, and Cdk1+ cycling cells (Fig. 1a). Additionally, we annotated cell cluster 15 as macrophages undergoing differentiation into osteoclasts. Mesenchymal cells, the largest population, contained multiple clusters, and were further annotated as osteoblasts, odontoblasts, pre-odontoblasts, periodontal ligament cells, apical pulp, dental pulp, cementoblasts, and dental follicle cells (Fig. 2a). We identified two clusters of dental follicle cells, which shared a similar gene expression profile, and both expressed stem cell marker genes Kit [30] and Nes [31]. This result is consistent with the current understanding that the dental follicle is considered a stem niche in adults [32]. Cell markers of each cell type have been previously presented [14,15,16, 18, 33,34,35,36,37]. To date, numerous studies have sought to identify the genes involved in tooth development. To identify the specific factors in dentin and cementum formation, we identified previously uncharacterized and preferentially expressed genes in tooth germ: Pclo in dental follicle cells, Fst and Igfbp2 in periodontal ligament cells, and Igfbp3 in apical pulp (Fig. 2). Furthermore, we validated the expression of Igfbp3 in apical pulp using immunofluorescence. Compared to other organs, teeth lack comprehensive references to annotate the various cell types [14]. As more cell markers are identified, this will be addressed and a clear hierarchy established.

RNA velocity enables the identification of genes that display pronounced dynamic behavior, which are candidates for important drivers of the main process in the population [38]. We constructed cell developmental and differentiation trajectories of rat molar cells and analyzed driver genes that related to cell fate using RNA velocity (Fig. 3, Additional file 1: Fig. S6). In the dynamic model, we identified three cell populations of differentiation origin: two types of dental follicle cells and apical pulp cells. These populations eventually differentiate into osteoblasts, cementoblasts, odontoblasts, or functional dental pulp cells. These differentiation trajectories are consistent with the reported lineages of dental mesenchymal cells [39]. Cluster-specific driver genes, Ptn and Satb2, are crucial for apical pulp differentiation. These genes have been reported to be essential for odontogenic differentiation of mesenchymal stem cells [40, 41]. We found that the two types of dental follicle cells share some of the same cluster-specific driver genes (Fig. 3cd). The pySCENIC analysis results also revealed distinct regulons among different cell types, while functionally similar cell types shared similar regulons (Additional file 4). This suggests that regulons reflect the differences between cell clusters to some extent, and the functions of different cell clusters might be determined by different transcription factor regulatory networks. A deeper investigation into the differences between the two types of dental follicle cells might provide better insights into bone and cementum development. Combining the results of RNA velocity and pySCENIC for more in-depth analysis is crucial for refining the regulatory network of dental mesenchymal cell development and for discovering new genes associated with cell differentiation.

Tnc is one of the top 300 driver genes both in M1 and M2. Among PN1, PN5, and PN10, Tnc has the highest expression level on PN5, primarily in odontoblasts and dental mesenchyme (Fig. 4), consistent with previous studies [42, 43]. We took the Tnc gene as an example to verify the function of driver genes of rats in the hDPSCs mineralization process. In the present study, we first confirmed the expression of TNC was upregulated in the odontogenic differentiation of hDPSCs. Compared with shCtrl infection, shTNC lentivirus infection decreased calcium nodule formation and DSPP and RUNX2 expression. Our results suggest potential applicability to humans; however, more functional analyses of these cluster-specific driver genes should be performed. In future studies, these driver genes could be prioritized to provide a comprehensive understanding of the mechanisms underlying dental hard tissue hypoplasia diseases.

Finally, we addressed the heterogeneity of macrophages in the first and second molars. A previous study has strengthened the role of immune cells in tooth development. Besides defense against pathogens, dental immune cells regulate dental development by secreting ligands that act on other dental cell types [18]. In the first molar, which was in the more mature development stage, the proportion of macrophages was more than that in the second molar. Previous studies have shown that the maximal number of osteoclast formation in the rat mandibular first molar is on postnatal day 3, while a secondary minor burst of osteoclastogenesis occurs around postnatal day 9 [44]. In our study, the expression of osteoclast differentiation factors was increased in the macrophages of the first molar on postnatal day 5, compared to the second molar (Fig. 6c). Studying how these factors function could further understand the process of bone resorption and root formation during the hard tissue formation period.

Conclusions

Incisors have long been used as models for studying continuous tooth development, whereas much less is known about molar development. In this study, we provided a detailed cell atlas of the rat mandibular molars at the single-cell level and identified key genes involved in dental mesenchymal cell differentiation. Our findings provide potential targets for diagnosing dental hard tissue diseases like dentin hypoplasia diseases and tooth regeneration, but further detailed investigation is required.

Methods

Single-cell RNA sequencing

Animal experiments were conducted with approval from the Southern Medical University Laboratory Animal Welfare and Ethics Committee. Rats were anesthetized with sodium pentobarbital and sacrificed by cervical dislocation. Three PN5 wild-type SD rats (the Laboratory Animal Center, Southern Medical University) from the same litter with similar body weights were chosen to extract single-cell suspensions. Mandibles were dissected under a stereomicroscope, and the first and second molar tooth germs of the left mandibular were carefully isolated and cut into small pieces. Subsequently, the small pieces were digested in 3 mg/mL collagenase type I (Sigma-Aldrich, St. Louis, MO, USA) and 3 mg/mL dispase II (Sigma-Aldrich) and 20 U/ml Dnase I (Solarbio, Beijing, China) at 37 ℃ for 1 h, followed by filtration through a 40-μm cell strainer (BD Biosciences, New Jersey, USA). After washing with cold PBS, the cells were resuspended in PBS containing 0.04% BSA. Upon confirming the quality of the cell suspensions, reverse transcription sequencing was immediately proceeded. The 10 × Genomics Chromium single-cell v3.0 reagent was used to construct the cDNA library by following the manufacturer’s protocol. The resulting libraries were sequenced on an Illumina NovaSeq 6000 System. Raw sequencing data have been uploaded to the Gene Expression Omnibus (GEO) database [45] (accession code GSE217465).

Data cleaning, normalization, and scaling

The rat (Rattus norvegicus) reference genome (Rnor_6.0) was downloaded from Ensemble. The concatenated gene-cell barcode matrix which was generated by using the official 10 × Genomics pipeline Cell Ranger v3.1.0 was imported into Seurat v3.1.1 [46], a toolkit for single-cell RNA-seq data analysis, for data processing. To exclude genes that might be detected from random noise, we filtered genes whose expression was detected in fewer than 10 cells. To exclude poor-quality cells that might result from multiplets or other technical noise, we filtered cells that were considered outliers (> third quartile + 1.5 × interquartile range or < first quartile − 1.5 × interquartile range) based on the number of expressed genes detected, the sum of UMI counts and the proportion of mitochondrial genes. In addition, we limited the proportion of mitochondrial genes to a maximum of 0.2 to further remove potential poor-quality data from broken cells. DoubletFinder v3 [47] was used to predict doublets in single-cell RNA sequencing data. After removing unwanted cells from the dataset, we employ a global-scaling normalization method “LogNormalize” that normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result.

Dimensional reduction, clustering, and visualization

Uniform Manifold Approximation and Projection (UMAP) [48] dimensional reduction was performed on the scaled matrix (with most variable genes only) using the first 20 components of principal component analysis (PCA) [49] to obtain a two-dimensional representation of the cell states. Cell clustering was performed using the FindClusters function. DEGs between two groups of cells were detected with a likelihood-ratio test implemented in “FindMarkers” function.

For cell annotation, we collected makers from Cell Marker database [50] and studies related to single-cell sequencing of the tooth [14,15,16, 18, 33,34,35,36,37]. We used two methods below to perform cell annotation. First, we conducted gene set enrichment analysis on DEGs in each identified cell cluster to detect if any marker genes of typical cell types were enriched. Gene set enrichment analysis was performed using ClusterProfiler (v3.14.3) with default parameters [51]. Second, we used AddModuleScore in Seurat to detect if one cell was high-scored with any marker genes of typical cell types. We combined these two methods. If the DEGs of a cell cluster were significantly enriched in the marker genes of a cell type in the CellMarker database, the cell cluster was defined as this cell type (BH-adjusted p-value < 0.05). Also, metascape was used to perform functional enrichment analysis [52].

RNA velocity analysis

RNA velocity aims to infer directed differentiation trajectories from snapshot single-cell transcriptomic data. It predicts the differentiation trajectories and state transitions of cells by analyzing the rates of change in gene expression. Considering the potential batch effects between samples, we performed RNA velocity analysis on samples M1 and M2 separately using the dynamic model in the scVelo 0.3.2 software [38]. In the analysis, cell cluster information and UMAP information were derived from the integrated analysis of M1 and M2. Based on scVelo, genes with pronounced dynamic behavior often contain high likelihood, and these genes are defined as driver genes. In the dynamic model, we calculated cluster-specific driver genes using the function rank_dynamical_genes. For the aforementioned analysis, all parameters were set to default.

Analysis of regulatory network transcription factors

Considering the potential batch effects between samples, we performed TF regulatory analysis on M1 and M2 using pySCENIC (v0.9.19) [53] as follows. We used “pyscenic grn” for network inference between TFs and targets, “pyscenic ctx” for regulon prediction, and “pyscenic auc” to evaluate the activity of each regulon across all cells. We used regulon_specificity_scores to identify cell cluster-specific regulons. All parameters were set to their default values.

Cell culture and odontoblastic differentiation

Third molars from healthy 18- to 22-year-old donors were collected at the Department of Stomatology, Nanfang Hospital, Guangzhou, China. The isolation of human dental pulp stem cells (hDPSCs) was performed as described elsewhere [54]. The cells were cultured at 37 ℃ and in a 5% CO2 atmosphere in high-glucose Dulbecco’s modified Eagle’s medium (DMEM; Gibco, Thermo Fisher Scientific, California, U.S.A.) supplemented with 10% fetal bovine serum (FBS; Gibco), hereafter referred to as the growth medium (GM). For the odontoblastic differentiation experiments, the cells were cultured in odontogenic medium (OM) for 7 and 14 days. The OM consisted of DMEM, 10% FBS, 50 mg/ml ascorbic acid (Sigma-Aldrich), 5 mM b-glycerophosphate (Sigma-Aldrich), and 10 nM dexamethasone (Sigma-Aldrich) [55].

Lentivirus infection

The TNC-specific shTNC1, shTNC2, and control shRNA (shCtrl) were designed, synthesized, and packaged into lentivirus by Ubigene (Guangzhou, China), with sequences provided in Additional file 1: Table S2. The hDPSCs were seeded in 12-well plates with a density of 1 × 105 cells/well and grown for 24 h. They reached approximately 60–70% confluence at the time of infection. The cells were infected at a multiplicity of infection (MOI) of 20 with 5 mg/mL polybrene for 12 h. The knockdown efficiency was evaluated using quantitative real-time polymerase chain reaction (RT-qPCR) and western blot analysis. To confirm the reproducibility of the results, infection, RT-qPCR, and western blot were repeated three times.

RT-qPCR

Total RNA was extracted from hDPSCs using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and reverse transcribed into cDNA using the HiScript III RT SuperMix for qPCR (+ gDNA wiper) kit (Vazyme, Nanjing, China). RT-qPCR was carried out with ChamQ SYBR qPCR Master Mix (Vazyme) on a LightCycler 480 (Roche, Indianapolis, IN, USA). Gene expression was quantified using the 2–ΔΔCT method and normalized to GAPDH mRNA levels. The primer pairs used for quantitation of TNC, DSPP, RUNX2, and GAPDH mRNA expression are listed in Additional file 1: Table S2. RT-qPCR was repeated three times for each test.

Western blot analysis

Total protein (30 μg) was separated on a 10% SDS–polyacrylamide gel and transferred to a polyvinylidene difluoride (PVDF) membrane (Millipore, Massachusetts, USA). The membranes were blocked for 1 h with 5% skim milk and then incubated overnight at 4 °C with primary antibody (anti-GAPDH (Proteintech, Wuhan, China), anti-TNC (Abcam, Cambridge, UK), and anti-RUNX2 (Abcam); anti-DSPP (Solarbio). The next day, membranes were incubated with the corresponding secondary antibodies at room temperature for 1 h after washing with Tris-buffered saline-Tween 20 (TBST) three times. The immunoreactive proteins were visualized with an ECL Kit (Tanon Science & Technology, Shanghai, China) according to the manufacturer’s instructions. The intensities of the protein bands were quantified using ImageJ software. Western blot analysis was repeated three times for each test.

Alizarin Red S (ARS) staining

The number of calcium nodules formed by hDPSCs after transfected with shCtrl or shTNC lentiviruses was analyzed by ARS staining. When cells were 70% confluent, the GM was replaced with the OM to induce odontogenic differentiation. After 14 days, the induced cells were fixed for 30 min at room temperature in 4% paraformaldehyde and then stained for 30 min with 1% ARS (Leagene, Beijing, China). ARS staining was repeated three times for each test.

Immunohistochemical and immunofluorescence analysis

Immunohistochemical and immunofluorescence studies were performed on 4-µm-thick unstained sections generated from formalin-fixed, decalcified, paraffin-embedded rat mandible (n = 3). For immunohistochemical staining, the sections were incubated with anti-TNC antibody (1:200; Abcam, Cambridge, UK) and observed using an Olympus VS200 microscope (Olympus; Tokyo, Japan). For immunofluorescence, the sections were incubated with anti-Igfbp3 antibody (1:100; Absin, Shanghai, China) diluted in 3% BSA/PBS overnight at 4℃. The next day samples were incubated in Alexa-fluor 488 Goat anti-Rabbit secondary antibodies (Proteintech). Nuclei were counterstained DAPI. Images were acquired using the Olympus VS200 microscope.

Statistical analysis

Statistical analyses were performed with GraphPad Prism 9.0. All data were presented as mean ± standard deviation. A t-test was used to assess differences between two groups, while one-way ANOVA was utilized for comparing differences among multiple groups. p-values < 0.05 were considered statistically significant.

Availability of data and materials

The datasets generated during and analyzed during the current study are available in the Gene Expression Omnibus (GEO) database (accession code GSE217465) [45].

Abbreviations

scRNA-seq:

Single-cell RNA-sequencing

PN:

Postnatal days

UMAP:

Uniform Manifold Approximation and Projection

DEGs:

Differentially expressed genes

hDPSCs:

Human dental pulp stem cells

GM:

Growth medium

OM:

Odontogenic medium

ARS:

Alizarin Red S

References

Bei M. Molecular genetics of tooth development. Curr Opin Genet Dev. 2009;19:504–10. https://doi.org/10.1016/j.gde.2009.09.002.

Balic A, Thesleff I. Tissue interactions regulating tooth development and renewal. Curr Top Dev Biol. 2015;115:157–86. https://doi.org/10.1016/bs.ctdb.2015.07.006.

Duan P, Bonewald LF. The role of the wnt/beta-catenin signaling pathway in formation and maintenance of bone and teeth. Int J Biochem Cell Biol. 2016;77:23–9. https://doi.org/10.1016/j.biocel.2016.05.015.

Hosoya A, Shalehin N, Takebe H, Shimo T, Irie K. Sonic hedgehog signaling and tooth development. Int J Mol Sci. 2020;21:1587. https://doi.org/10.3390/ijms21051587.

Li C, Prochazka J, Goodwin AF, Klein OD. Fibroblast growth factor signaling in mammalian tooth development. Odontology. 2014;102:1–13. https://doi.org/10.1007/s10266-013-0142-1.

Shi C, Yuan Y, Guo Y, Jing J, Ho TV, Han X, et al. BMP signaling in regulating mesenchymal stem cells in incisor homeostasis. J Dent Res. 2019;98:904–11. https://doi.org/10.1177/0022034519850812.

Fons RJ, Star H, Lav R, Watkins S, Harrison M, Hovorakova M, et al. The impact of the Eda pathway on tooth root development. J Dent Res. 2017;96:1290–7. https://doi.org/10.1177/0022034517725692.

Yang G, Zhou J, Teng Y, Xie J, Lin J, Guo X, et al. Mesenchymal TGF-beta signaling orchestrates dental epithelial stem cell homeostasis through Wnt signaling. Stem Cells. 2014;32:2939–48. https://doi.org/10.1002/stem.1772.

Niwa T, Yamakoshi Y, Yamazaki H, Karakida T, Chiba R, Hu JC, et al. The dynamics of TGF-beta in dental pulp, odontoblasts and dentin. Sci Rep. 2018;8:4450. https://doi.org/10.1038/s41598-018-22823-7.

Thesleff I. The genetic basis of tooth development and dental defects. Am J Med Genet A. 2006;140A:2530–5. https://doi.org/10.1002/ajmg.a.31360.

Bailleul-Forestier I, Molla M, Verloes A, Berdal A. The genetic basis of inherited anomalies of the teeth. Part 1: clinical and molecular aspects of non-syndromic dental disorders. Eur J Med Genet. 2008;51:273–91. https://doi.org/10.1016/j.ejmg.2008.02.009.

Lan Y, Jia S, Jiang R. Molecular patterning of the mammalian dentition. Semin Cell Dev Biol. 2014;25–26:61–70. https://doi.org/10.1016/j.semcdb.2013.12.003.

Hwang B, Lee JH, Bang D. Single-cell RNA sequencing technologies and bioinformatics pipelines. Exp Mol Med. 2018;50:1–14. https://doi.org/10.1038/s12276-018-0071-8.

Fresia R, Marangoni P, Burstyn-Cohen T, Sharir A. From bite to byte: dental structures resolved at a single-cell resolution. J Dent Res. 2021;100:897–905. https://doi.org/10.1177/00220345211001848.

Chiba Y, Saito K, Martin D, Boger ET, Rhodes C, Yoshizaki K, et al. Single-cell RNA-sequencing from mouse incisor reveals dental epithelial cell-type specific genes. Front Cell Dev Biol. 2020;8:841. https://doi.org/10.3389/fcell.2020.00841.

Gong X, Zhang H, Xu X, Ding Y, Yang X, Cheng Z, et al. Tracing PRX1(+) cells during molar formation and periodontal ligament reconstruction. Int J Oral Sci. 2022;14:5. https://doi.org/10.1038/s41368-021-00155-z.

He P, Zheng L, Zhou X. IGFs in dentin formation and regeneration: progress and remaining challenges. Stem Cells Int. 2022;2022:3737346. https://doi.org/10.1155/2022/3737346.

Shi Y, Yu Y, Zhou Y, Zhao J, Zhang W, Zou D, et al. A single-cell interactome of human tooth germ from growing third molar elucidates signaling networks regulating dental development. Cell Biosci. 2021;11:178. https://doi.org/10.1186/s13578-021-00691-5.

Madureira DF, Taddei SA, Abreu MH, Pretti H, Lages EM, Da ST. Kinetics of interleukin-6 and chemokine ligands 2 and 3 expression of periodontal tissues during orthodontic tooth movement. Am J Orthod Dentofacial Orthop. 2012;142:494–500. https://doi.org/10.1016/j.ajodo.2012.05.012.

Yu X, Huang Y, Collin-Osdoby P, Osdoby P. CCR1 chemokines promote the chemotactic recruitment, RANKL development, and motility of osteoclasts and are induced by inflammatory cytokines in osteoblasts. J Bone Miner Res. 2004;19:2065–77. https://doi.org/10.1359/JBMR.040910.

Iezaki T, Fukasawa K, Park G, Horie T, Kanayama T, Ozaki K, et al. Transcriptional modulator Ifrd1 regulates osteoclast differentiation through enhancing the NF-kappaB/NFATc1 pathway. Mol Cell Biol. 2016;36:2451–63. https://doi.org/10.1128/MCB.01075-15.

Omar I, Guterman-Ram G, Rahat D, Tabach Y, Berger M, Levaot N. Schlafen2 mutation in mice causes an osteopetrotic phenotype due to a decrease in the number of osteoclast progenitors. Sci Rep. 2018;8:13005. https://doi.org/10.1038/s41598-018-31428-z.

Zhang X, Ning T, Wang H, Xu S, Yu H, Luo X, et al. Stathmin regulates the proliferation and odontoblastic/osteogenic differentiation of human dental pulp stem cells through Wnt/beta-catenin signaling pathway. J Proteomics. 2019;202:103364. https://doi.org/10.1016/j.jprot.2019.04.014.

Jin L, Gao F, Zhang L, Wang C, Hu L, Fan Z, et al. Pleiotropin enhances the osteo/dentinogenic differentiation potential of dental pulp stem cells. Connect Tissue Res. 2021;62:495–507. https://doi.org/10.1080/03008207.2020.1779238.

Song YL, Bian Z. Gene mutations and disorders of dental hard tissues. Zhonghua Kou Qiang Yi Xue Za Zhi. 2020;55:316–22. https://doi.org/10.3760/cma.j.cn112144-20200225-00086.

Cobourne MT, Sharpe PT. Diseases of the tooth: the genetic and molecular basis of inherited anomalies affecting the dentition. Wiley Interdiscip Rev Dev Biol. 2013;2:183–212. https://doi.org/10.1002/wdev.66.

Smith C, Poulter JA, Antanaviciute A, Kirkham J, Brookes SJ, Inglehearn CF, et al. Amelogenesis imperfecta; genes, proteins, and pathways. Front Physiol. 2017;8:435. https://doi.org/10.3389/fphys.2017.00435.

Chen D, Li X, Lu F, Wang Y, Xiong F, Li Q. Dentin dysplasia type I-A dental disease with genetic heterogeneity. Oral Dis. 2019;25:439–46. https://doi.org/10.1111/odi.12861.

Lungova V, Radlanski RJ, Tucker AS, Renz H, Misek I, Matalova E. Tooth-bone morphogenesis during postnatal stages of mouse first molar development. J Anat. 2011;218:699–716. https://doi.org/10.1111/j.1469-7580.2011.01367.x.

Foster BM, Zaidi D, Young TR, Mobley ME, Kerr BA. CD117/c-kit in cancer stem cell-mediated progression and therapeutic resistance. Biomedicines. 2018;6:31. https://doi.org/10.3390/biomedicines6010031.

Nouri P, Zimmer A, Bruggemann S, Friedrich R, Kuhn R, Prakash N. Generation of a NES-mScarlet red fluorescent reporter human iPSC line for live cell imaging and flow cytometric analysis and sorting using CRISPR-Cas9-mediated gene editing. Cells. 2022;11:268. https://doi.org/10.3390/cells11020268.

Rusu MC, Rascu A, Stoenescu MD. Endothelial expression of c-kit and CD68 in dental follicles of human impacted third molars. Folia Morphol (Warsz). 2018;77:485–8. https://doi.org/10.5603/FM.a2017.0108.

Chiba Y, Yoshizaki K, Tian T, Miyazaki K, Martin D, Saito K, et al. Integration of single-cell RNA- and CAGE-seq reveals tooth-enriched genes. J Dent Res. 2021:503863097. https://doi.org/10.1177/00220345211049785.

Krivanek J, Soldatov RA, Kastriti ME, Chontorotzea T, Herdina AN, Petersen J, et al. Dental cell type atlas reveals stem and differentiated cell types in mouse and human teeth. Nat Commun. 2020;11:4816. https://doi.org/10.1038/s41467-020-18512-7.

Takada K, Chiba T, Miyazaki T, Yagasaki L, Nakamichi R, Iwata T, et al. Single cell RNA sequencing reveals critical functions of Mkx in periodontal ligament homeostasis. Front Cell Dev Biol. 2022;10:795441. https://doi.org/10.3389/fcell.2022.795441.

Hu H, Duan Y, Wang K, Fu H, Liao Y, Wang T, et al. Dental niche cells directly contribute to tooth reconstitution and morphogenesis. Cell Rep. 2022;41:111737. https://doi.org/10.1016/j.celrep.2022.111737.

Wang Y, Zhao Y, Chen S, Chen X, Zhang Y, Chen H, et al. Single cell atlas of developing mouse dental germs reveals populations of CD24(+) and Plac8(+) odontogenic cells. Sci Bull (Beijing). 2022;67:1154–69. https://doi.org/10.1016/j.scib.2022.03.012.

Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat Biotechnol. 2020;38:1408–14. https://doi.org/10.1038/s41587-020-0591-3.

Jussila M, Thesleff I. Signaling networks regulating tooth organogenesis and regeneration, and the specification of dental mesenchymal and epithelial cell lineages. Cold Spring Harb Perspect Biol. 2012;4:a008425. https://doi.org/10.1101/cshperspect.a008425.

Zhang L, Xia D, Wang C, Gao F, Hu L, Li J, et al. Pleiotrophin attenuates the senescence of dental pulp stem cells. Oral Dis. 2023;29:195–205. https://doi.org/10.1111/odi.13929.

Chen Q, Zheng L, Zhang Y, Huang X, Wang F, Li S, et al. Special AT-rich sequence-binding protein 2 (Satb2) synergizes with Bmp9 and is essential for osteo/odontogenic differentiation of mouse incisor mesenchymal stem cells. Cell Prolif. 2021;54:e13016. https://doi.org/10.1111/cpr.13016.

Vainio S, Jalkanen M, Thesleff I. Syndecan and tenascin expression is induced by epithelial-mesenchymal interactions in embryonic tooth mesenchyme. J Cell Biol. 1989;108:1945–53. https://doi.org/10.1083/jcb.108.5.1945.

Midwood KS, Chiquet M, Tucker RP, Orend G. Tenascin-C at a glance. J Cell Sci. 2016;129:4321–7. https://doi.org/10.1242/jcs.190546.

Shiyan H, Nanquan R, Shuhao X, Xiaobing L. Research progress on the cellular and molecular mechanisms of tooth eruption. Hua Xi Kou Qiang Yi Xue Za Zhi. 2016;34:317–21.

Zheng Y, Lu T, Zhang L, Gan Z, Li A, He C, He F, He S, Zhang J, Xiong F. Single-cell RNA-seq analysis of rat molars reveals cell-identity and driver genes associated with dental mesenchymal cell differentiation. NCBI GEO accession: GSE217465. 2024. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE217465.

Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 2015;161:1202–14. https://doi.org/10.1016/j.cell.2015.05.002.

Mcginnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8:329–37. https://doi.org/10.1016/j.cels.2019.03.003.

Mcinnes L, Healy J. UMAP: uniform manifold approximation and projection for dimension reduction. J Open Source Softw. 2018;3:861. https://doi.org/10.21105/joss.00861.

Jolliffe IT. Principal component analysis. J Mark Res. 2002;87:513. https://doi.org/10.2307/3172953.

Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, et al. Cell Marker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res. 2019;47:D721-28. https://doi.org/10.1093/nar/gky900.

Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7. https://doi.org/10.1089/omi.2011.0118.

Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10:1523. https://doi.org/10.1038/s41467-019-09234-6.

Van de Sande B, Flerin C, Davie K, De Waegeneer M, Hulselmans G, Aibar S, et al. A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nat Protoc. 2020;15:2247–76. https://doi.org/10.1038/s41596-020-0336-2.

Ma D, Yu H, Xu S, Wang H, Zhang X, Ning T, et al. Stathmin inhibits proliferation and differentiation of dental pulp stem cells via sonic hedgehog/Gli. J Cell Mol Med. 2018;22:3442–51. https://doi.org/10.1111/jcmm.13621.

Fu T, Liu Y, Huang X, Guo Y, Shen J, Shen H. lncRNA SNHG1 regulates odontogenic differentiation of human dental pulp stem cells via miR-328-3p/Wnt/beta-catenin pathway. Stem Cell Res Ther. 2022;13:311. https://doi.org/10.1186/s13287-022-02979-w.