Single-cell composition of human tooth germ
We isolated tooth germ tissue from two patients with different developmental statuses of left mandibular third molars (Fig. 1A and Figure S1). One patient's left mandibular third molar was at developmental stage A (calcification of cusp tips without coalescence of other calcifications). The other patient’s was at stage D (complete crown formation up to cementoenamel junction). The developmental status of the third molars was assessed using eight-stage developmental scoring (from A to H) proposed by Demirjian et al.[21]. We applied BD-seq on the dissected tooth germ and obtained transcriptome data for 9,855 cells, which on average contained about 28,000 mapped reads per cell, after RNA quantification and quality control. Using Louvain method embedded in Seurat 3.0 R package[22], we partitioned all cells into 11 clusters (Fig. 1B). Various immune cells, including T cell (CD3E+), neutrophil (S100A9+), macrophage (CCL3+), monocytes (FCN1+) and dendritic cell (CD1C+) (Fig. 1B-D and Table S1), etc., consisted of nearly 83% of all cells. In the remaining cells, we identified SPARC+RUNX2+ osteoblast[23], ACP5+ osteoclast[24], RGS5+ pericytes and VWF+ endothelium. Lastly, we identified a population of SOX9+ cells that exhibited heterogeneous transcriptome characteristics (Table S1).
To resolve this heterogeneity, we carried out further clustering analysis on the SOX9+ cells. As shown in Fig. 1E, SOX9+ cells consisted of two subpopulations: one expressed apical papilla stem cell (APSC) marker CD24[25], and another expressed ameloblast marker AMBN and epithelium-associated gene CLU[26]. We therefore separated SOX9+ cells into APSC and ameloblast. Similarly, osteoblast (Fig. 1F) also consisted of two subpopulations, namely, immature and differentiated osteoblast (iOsteoblast and dOsteoblast, respectively). iOsteoblast expressed higher level of VIM that inhibited osteoblast differentiation[27], whereas dOsteoblast highly expressed SPARC that took part in osteogenesis[23] and GJA1 that took part in osteoblast differentiation[28].
Since SPARC and GJA1 expression is not limited to osteoblast in dental tissue, we further applied immunofluorescence to elucidate their spatial and cellular distribution in tooth germ. On the tip of dental sac (Fig. 1F and Figure S2), which is proximal to the bone interface, SPARC and GJA1 showed co-localization in the osteoblast. On the outer surface of dental papilla (Fig. 1G and Figure S2), SPARC expressed in the odontoblast at the root end while GJA1 expressed in the odontalblast at the crown end, with little co-localization at the middle. This result supported our classification of SPARC+GJA1+ cells as the mature osteoblast.
T cell subpopulations and their intercellular interaction patterns
Our single-cell data revealed that more than 42% of tooth germ cells were T cell (Fig. 1C), which presumably contained diverse subpopulations with important roles in tooth structure[29]. To analyze the role of these subpopulations, we applied Seurat cluster analysis on T cells and obtained eight sub-clusters (Fig. 2A-D). We first defined the cytotoxic T cells by NKG7 and GNLY expression, and separated them into natural killer (NK) T and CD8T according to CD8A expression (Fig. 2D). We then defined memory T and Th17 by IL7R, CCR6 and CCR7 expression[30, 31]. Finally, we found small subpopulation of CRTAM+ activated CD8T, SELL+ naïve NK, and MKI67+ proliferation T cell. Separated by developmental period, we found that stage A tooth germ contained more Th17 (74% of Th17 came from stage A tooth germ) but less cytotoxic CD8T (17%, Fig. 2B). The proportion of tooth-residence memory cells was also higher in stage D tooth germ, reflected by higher expression of residence T cell marker CD69[32] (permutation p < 4.1×10− 5, Fig. 2C).
We then applied CellPhoneDB[33] analysis to explore the subpopulation-specific intercellular signal transduction from T cell to other tooth germ cells. By summarizing the ligand-receptor pairs reaching significant threshold (Method), all T cell subpopulations showed strongest association with ameloblast (Number of significant pairs = 2 to 9, association strength > 12. Figure 2E), especially by CD74-to-APP and HLADRB1-to-OGN signaling pathways (Fig. 2F). Naïve NK and activated CD8T showed the most significant communication with ameloblast, and they exhibited subpopulation-specific pathways like CXCR6-to-CXCL16 and CRTAM-to-CADM1. T cell also exhibited strong communication with osteoclast (Number of significant pairs = 4 to 11, association strength = 2.5 to 6.1. Figure 2E), especially by signals from CCL3/CCL4/CCL5 to CCR1/CCR5 (Fig. 2G). Interestingly, we observed CTLA4-to-CD86 signaling which was unique to Th17 cell, where CD86 was known to suppress osteoclast differentiation[34]. For other immune cells, T cell subpopulations showed divergent association with monocyte (Number of significant pairs = 2 to 8, association strength = 1.6 to 12.3. Figure 2E). These results suggested that human tooth germ contained diverse T cell subpopulation with distinct cell interaction patterns.
Non-T immune cell subpopulations and their intercellular interaction patterns
Despite T cell subpopulations, other immune cells also play an important role in tooth development[35]. Since the proportion of neutrophil was profoundly larger than remaining immune cells (Fig. 1C), we first applied clustering analysis to dissect neutrophil subpopulation. As shown in Fig. 3A-C, we obtained eight subpopulations of neutrophil. They were classified by unique expression of genes related to neutrophil functions (Fig. 3C). We first identified PGLYRP1+ neutrophil with specific roles in innate immunity[36], as well as MX1+ antiviral neutrophil[37], and SLPI+ inhibitory neutrophil[38]. Another three subpopulations were labeled by P2RY13, PRRG4 and S100P, all of which took part in neutrophil functions. Separated by developmental period, stage A tooth germ contained more S100P+ neutrophils (71%) and less P2RY13+ (17%) and antiviral (9%) neutrophils.
We then applied CellPhoneDB[33] to analyze the signal transduction from non-T immune cells to dental cells. For neutrophil subpopulation (Fig. 3D and E), they were overwhelmingly associated with endothelium (Number of significant pairs = 6 to 10, association strength > 30; strength between neutrophil and other cells < 18), especially via ACKR1 which guided neutrophil migration[39]. Antiviral neutrophil showed the strongest association with endothelium (Number of significant pairs = 10) and showed subpopulation-specific signaling of CEACAM1-to-SELE, which mediated neutrophil activation and migration[39]. For other immune cells, we found that B cells showed the strongest signaling transduction to dental cells (Number of significant pairs = 2 to 7, association strength > 25), especially the CD74-to-APP pathway.
Intercellular signaling network regulated osteoblast maturation
Having resolved all subpopulation for tooth germ cell types, we now had the opportunity to analyze the functional implication of intercellular signaling network. We started by delineating the process of osteoblast maturation, then analyzed whether this process was regulated by signals from other cells by using ligand-target[40] and transcription factor regulation network[41].
As shown in Fig. 4A, we applied monocle pseudotime analysis on the osteoblast and observed a clear immature-to-differentiated lineage. We found that the expression of ALPL, BGN and OGN, genes related to mature osteoblast function[42], increased rapidly at the beginning of this lineage (Fig. 4A and Figure S3). Other genes related to extracellular matrix formation such as MMP2 and COL12A1 also showed gradual increment throughout the lineage (Figure S3). On the other hand, genes regulating the proliferation of osteoblast, such as SCUBE1 and PTCH1 (Fig. 4A), gradually decreased during this lineage. Taken all genes showing significant differential expression along pseudotime, we found a clear temporal cascade (Fig. 4B). By Gene Ontology (GO) analysis, we found that genes showing decreasing expression in the cascade mainly took part in Wnt signaling pathway (adjusted p value of GO analysis [GO P] = 4.41×10− 11, odds ratio [OR] = 5.21), mesenchymal development (GO P = 1.91×10− 8, OR = 5.08),BMP signaling pathway (GO P = 1.16×10− 7, OR = 6.12, Fig. 4B and Table S2), etc. Inversely, ascending genes took part in extracellular matrix organization (GO P = 5.24×10− 31, OR = 8.46), bone growth (GO P = 6.65×10− 8, OR = 13.49) and biomineralization (GO P = 2.49×10− 9, OR = 7.07, Fig. 4B and Table S3), etc. These results confirmed that the pseudotime analysis successfully reconstruct the process of osteoblast differentiation, and we therefore highlighted all genes significantly altered along pseudotime as key osteoblast lineage genes (Table S4).
We hypothesis that these key lineage genes might be downstream targets of intercellular signaling networks, which regulated the process of osteoblast differentiation. Thus, we applied nichenetr network analysis[40] to prioritize the intercellular ligands and pathways that might be upstream of these key genes. Nichenetr prioritized 20 potential ligands that could regulate the expression of these key genes (Fig. 4C and Figure S3), such as IL1A and IL1B that mainly secreted from macrophage, TNF and APOE that mainly secreted from monocyte, as well as IFNG and TFGB1 that mainly secreted from T cell, etc. BMP2 and BMP7, top prioritized ligands that are known to regulate osteoblast activity[43], were mainly expressed in APSC. We then inferred the potential receptors on osteoblast that mediated these ligand-target associations (Fig. 4C), and found that TGFBR1, FGFR3 and TGFBR2 were linked with most key genes. Interestingly, ascending key genes were also enriched in GO term "response to fibroblast growth factor" (GO P = 1.82×10− 7, OR = 6.53, Fig. 4B). We managed all nichenetr-inferred regulation relations into a circus plot, showing the intact network underlying osteoblast maturation (Fig. 4C).
Taking one step further, we asked whether this regulation network involved any key transcription factors (TF). We applied SCENIC[41] to infer the TF-gene regulation network and found that eight key osteoblast lineage genes (PTCH1, BMP7, EGR2, etc.) functioned as TF that could regulate other key osteoblast lineage genes (Fig. 4D). They were downstream to receptor such as LTBR, FRFR3, TGFBR2, BMPR1A and RARG. By combination of ligand-receptor-TF-target regulation results, we identified several important pathways such as FGF2-FGFR3-ID1-MMP2/MMP13/ALPL/TNC, TGFB1-TGFBR2-FOXQ1-HGF/FST/WNT5A, etc. Interestingly, we also observed that BMP could serve as both ligand and TF in the regulation of osteoblast differentiation, and regulated the expression of PTCH1, LTBP1, EDNRA, etc. Taken together, our result generated a gene regulation network originated from cell type-specific ligands, which could regulate the osteoblast differentiation.
BMP-FGFR1-MSX1 pathway had central role in APSC renewal regulation
Apical papilla stem cell (APSC) resides in human tooth germ and retains multipotent and proliferation capacity via consistent self-renewal. To identify essential signaling pathways that regulate APSC renewal, we first applied pseudotime analysis to resolve the cell cycle alteration of APSC. As shown in Fig. 5A, APSC exhibited gradual transmission along cell cycle: on the right branch, APSC within G2/M phase gradually left cell cycle and transmitted into resting state, whereas in the left branch, G0 resting APSC transmitted into G1/S phase and entered cell cycle again. We took the left branch as the APSC renewal process and applied differential expression analysis to identified genes that altered along this process, such as SERPING1, DKK3 and HSPA1B (Figure S4). As shown in Fig. 5B and Table S5-S6, genes that were up-regulated during renewal took part in mesenchymal cell proliferation (GO P = 3.67×10− 4, OR = 30.3), cell growth (GO P = 8.82×10− 5, OR = 6.89), etc. On the other hand, genes related to tumor necrosis factor bio-synthesis (GO P = 9.98×10− 3, OR = 21.4) and glycosaminoglycan catabolic (GO P = 8.77×10− 4, OR = 19.9) were down regulated during APSC renewal. We took all differential expression genes together and defined them as APSC renewal genes (Table S7).
We next applied Nichenetr to identify upstream intercellular signals that regulated the APSC renewal genes. As shown in Fig. 5C and Figure S4, ligands released from monocytes (IL1B, etc.), Osteoblast (BMP4, TGFB3, BMP5) and T cells (TGFB1) had the strongest regulation potentiality on APSC renewal genes. The autocrine of BMP2 and BMP7 on APSC also regulated a large number of renewal genes. These ligands mostly acted on APSC receptor FGFR1, which regulated 28 downstream renewal genes (Fig. 5C and Figure S4), including MSX1, PTCH1, SOX9, etc. Similar to the analysis of Osteoblast, we applied SCENIC analysis to highlight key TF in this ligand-target network (Fig. 5D). We found a transcription factor MSX1, which was a downstream target of FGFR1, regulated 47 renewal genes (Fig. 5D), especially VCAN, C1S and TIMP2 (regulation score > 1.5). By taking all MSX1 targets as a whole (so-called "regulon" [41]), we found that they were generally elevated during APSC transition to G1/S phase. These results highlighted the role of BMP-FGFR1-MSX1 pathway in APSC renewal. In addition, we also found other TF such as KLF4, ID3, JUN and EGR1, etc., that regulated other renewal genes like HSPA1B, CALR, SGK1, etc.
Transformation of Osteoclast is regulated by signals from Osteoblast and macrophage
In tooth and other bone tissues, monocytes and macrophage continuously transform into osteoclasts, the dysregulation of which might disrupt the bone remodeling balance. To find the signaling pathways that regulate this transformation, we first identified transformation-related genes by differential expression analysis (Fig. 6A). At the significance threshold of FDR-adjusted P < 0.01 and log fold change > 1, we found 111 genes that were elevated during monocyte-to-osteoblast transformation and 149 genes that were elevated during macrophage-to-osteoblast transformation (Table S8-S9). We merged these two gene lists into 183 unique transformation-related genes and applied nichenetr[40] to discover the upstream signaling pathways that regulated them (Fig. 6B). In accordance with previous studies, ligands from osteoblasts (BMP4, BMP5, TNFSF11, etc.) and macrophage (CCL3 and TNF) regulated the largest number of transformation-related genes (40 and 35, respectively, Figure S5). CSF1 secreted from pericyte also regulated 23 transformation-related genes via receptor CSF1R. Interestingly, the cellular distribution of receptors of these ligands were different: BPMR1A and CSF1R mainly expressed on macrophage and osteoclast, whereas TNFRSF1B and NOTCH1 mainly expressed on monocytes (Fig. 6B). This result indicated that monocyte-to-osteoblast and macrophage-to-osteoblast transformation was regulated by different signaling pathways. Nonetheless, many key genes involved in osteolysis, such as ACP5, NFATC1, MMP9, were regulated by multiple signaling pathways (Figure S5).
It has been revealed that neutrophil could activate osteoclasts and trigger osteonecrosis during inflammation. Since CellPhoneDB[33] analysis (Fig. 3D) found that neutrophil subtypes had dense connection with macrophage but few interactions with monocyte, we studied the neutrophil-to-macrophage signaling that could regulate the transformation-related genes. As shown in Fig. 6C, we found five ligand-receptor pathways that were significantly activated between neutrophil and macrophage. Among them, TNFSF14-LTBR pathway was shared by all seven neutrophil subtypes, whereas CEACAM1-CD209 was specific to MX1+ antiviral neutrophils. We further explored the downstream TF regulation networks of these pathways by nichenetr and SCENIC, and found three TF whose target genes significantly enriched in transformation-related genes: EGR1 (Fisher P = 8.90×10− 16, OR = 5.32), TCF4 (Fisher P = 5.94×10− 15, OR = 4.70) and SMARCA1 (Fisher P = 0.02, OR = 2.11, Fig. 6D). Interestingly, they were all downstream to TNFSF14-LTBR signaling pathway, supporting its essential role in osteoclast transformation. These TFs regulated multiple osteoclast maturation markers[44], such as CALR, MMP9, NFATC1, etc. We further applied GO analysis and found that the target genes of EGR1 significantly enriched in functions related to osteolysis, such as phagosome acidification (GO P = 3.11×10− 5, OR = 19.7) and receptor-mediated endocytosis (GO P = 3.26×10− 4, OR = 4.10). For SMARCA1, the most significant enrichment was found for regulation of osteoclast development (GO P = 0.01, OR = 39.5), which suggested that SMARCA1 may be one of the regulators of osteoclast transformation. For TCF4 targets, we did not observe significant functional enrichment.