A Single-cell Interactome of Human Tooth Germ Elucidates Signaling Networks Regulating Dental Development

Background: Development of dental tissue is regulated by extensive cell crosstalk based on various signaling molecules, such as BMP, FGF and SHH pathways. However, an intact network of the intercellular regulation is still lacking. Result: To gain an unbiased and comprehensive view of this dental cell interactome, we applied single-cell RNA-seq on immature human tooth germ of the third molar, discovered rened cell subtypes, and applied multiple network analysis to identify the central signaling pathways. We found that immune cells made up over 80% of all tooth germ cells, which exhibited profound regulation on dental cells via TGF-β, TNF and IL-1. During osteoblast differentiation, expression of genes related to extracellular matrix and mineralization was continuously elevated by signals from BMP and FGF family. As for the self-renewal of apical papilla stem cell, BMP-FGFR1-MSX1 pathway directly regulated the G0-to-S cell cycle transition. We also conrmed that CSF1 secreted from pericyte and TNFSF11 secreted from osteoblast regulated a large proportion of genes related to osteoclast transformation from macrophage and monocyte. Conclusion: We constructed the intercellular signaling networks that regulated the essential developmental process of human tooth, which served as a foundation for future dental regeneration engineering and the understanding of oral pathology.


Background
The human tooth and periodontal tissue emerge from the neural crest-derived ectomesenchyme of frontonasal, maxillary, and mandibular protrusions [1]. Tooth development is a long-term and complex biological process involving cell-cell and epithelial-mesenchymal interaction, cell differentiation, morphogenesis, tissue mineralization, and tooth eruption [2]. At the initial stage, cell proliferation activates in speci c areas of the dental lamina. The proliferative epithelium then extends to the deep connective tissue, and the terminal cells proliferate and further develop into the enamel organ. At the same time, the ectomesenchyme cells under the proliferative epithelium also proliferate rapidly and gather around the epithelium. These locally proliferated epithelia and mesenchyme together form the tooth germ [3]. The tooth germ consists of three parts: 1. Enamel organ originated from oral ectoderm and forms enamel; 2.
Dental papilla originated from ectomesenchyme, forms pulp and dentin; 3. Dental follicle originated from ectomesenchyme, forms cementum, periodontal ligament and alveolar bone [4]. Undoubtedly, a comprehensive understanding of tooth development requires dissection of these tooth germ substructures.
Like all developmental processes, tooth development is regulated by a series of complex gene cascades, which drive the cells to enter a predetermined location and differentiate in a speci c direction [2]. The formation of dental arch depends on the coordinated regulation of a variety of signaling molecules and location signals, which jointly regulate the development process of cell division rate, trend and direction of cell migration, cell differentiation and apoptosis. In the above process, a series of genes and signaling pathways play an important role in regulation, including Shh (Sonic hedgehog) pathway [5], Wnt pathway[6], FGF ( broblast growth factors) pathway [7], TGF-β (Transforming growth factor-beta) [8] pathway and BMPs (bone morphogenesis proteins) family [9][10][11]. These signaling molecules bind to corresponding receptors and regulate the expression of speci c genes. Speci c functions elicited by activation of these pathways are noted during distinct phases of dental tissue differentiation, some of which are bene cial for cell stemness and proliferation (FGF, Shh) while others such as Wnt, TGF-β, and BMPs act in postnatal differentiation phases and promote polarization, migration, and calci cation [1,2].
Understanding the intact signaling network in tooth development is essential to dental regeneration engineering and clinical dentistry. So far, many functional studies have elucidated the components and processes of speci c pathways [1,2,[11][12][13][14][15], and the application of these insights has promoted translational medicine and the understanding of oral disease. The rapid advancement of single-cell RNA sequencing (scRNA) and corresponding data analysis algorithms has provided a chance to draw this entire cell interactome. Scientists have applied scRNA to draw the cell atlas of mouth dental development at various stages [16][17][18] and human periodontal tissue [19,20]. However, the comprehensive cell interactome of human tooth is still lacking.
To achieve this goal, we applied scRNA to tooth germs isolated from developing third molar of healthy volunteers who planned to have orthodontic treatment (Fig. 1A). In human, the third molar has the latest development time scale of all teeth, and properly retains immature tooth germ structure before one's twenties, making it an ideal object for dental development research. In this study, we rst identi ed re ned dental cell subtypes in human tooth germ, discovered essential genes involved in dental cell differentiation and transformation, then constructed ligand-receptor-transcription factor networks that regulate these essential genes.

Method
Sample collection and pre-processing This study was approved and supervised by Ethical committee of Shanghai Ninth People's Hospital.
Written informed consent was provided by all participants. Both the surgical procedures were performed by the same surgeon and assistant with patients under local anesthesia. After a full-thickness mucoperiosteal ap elevation, the buccal bone of the third molar was removed with piezosurgery handpiece to expose the tooth germ, which was carefully enucleation by curette immediately. Where necessary, the crown sectioning were performed with a high-speed handpiece and ssure burs. After removing the mineralized part of it, the tooth germ was rinsed with normal saline and stored in MACS Tissue Storage Solution (Miltenyi, German), and was immediately delivered to Single-cell RNA-seq platform.

Single-cell Dissociation
Single-cell RNA-seq experiment was performed by experimental personnel in the laboratory of NovelBio Bio-Pharm Technology Co Ltd. The tissues were surgically removed and kept in MACS Tissue Storage Solution (Miltenyi Biotec) until processing. The tissue samples were processed as described below.
Brie y, samples were rst washed with phosphate-buffered saline (PBS), minced into small pieces (approximately 1mm3) on ice and enzymatically digested with 0.5 mg/mL collagenase I / (Worthington) and 50 U/mL DNase I (Worthington) for 45 min at 37°C, with agitation. After digestion, samples were sieved through a 70µm cell strainer, and centrifuged at 300g for 5 min. After the supernatant was removed, the pelleted cells were suspended in red blood cell lysis buffer (Miltenyi Biotec) to lyse red blood cells. After washing with PBS containing 0.04% BSA, the cell pellets were re-suspended in PBS containing 0.04% BSA and re-ltered through a 35μm cell strainer. Dissociated single cells were then stained for viability assessment using Calcein-AM (Thermo Fisher Scienti c) and Draq7 (BD Biosciences). The single-cell suspension was further enriched with a MACS dead cell removal kit (Miltenyi Biotec).
Single-cell RNA Sequencing BD Rhapsody system was used to capture the transcriptomic information of the (two sample-derived) single cells. Single-cell capture was achieved by random distribution of a single-cell suspension across >200,000 microwells through a limited dilution approach. Beads with oligonucleotide barcodes were added to saturation so that a bead was paired with a cell in a microwell. The cells were lysed in the microwell to hybridize mRNA molecules to barcoded capture oligos on the beads. Beads were collected into a single tube for reverse transcription and ExoI digestion. Upon cDNA synthesis, each cDNA molecule was tagged on the 5′ end (that is, the 3′ end of a mRNA transcript) with a unique molecular identi er (UMI) and cell barcode indicating its cell of origin. Whole transcriptome libraries were prepared using the BD Rhapsody single-cell whole-transcriptome ampli cation (WTA) work ow including random priming and extension (RPE), RPE ampli cation PCR and WTA index PCR. The libraries were quanti ed using a High Sensitivity DNA chip (Agilent) on a Bioanalyzer 2200 and the Qubit High Sensitivity DNA assay (Thermo Fisher Scienti c). Sequencing was performed by illumina sequencer (Illumina, San Diego, CA) on a 150 bp paired-end run.
Single-cell RNA Analysis scRNA-seq data analysis was performed by NovelBio Bio-Pharm Technology Co., Ltd. with NovelBrain Cloud Analysis Platform. We applied fastp [47] with default parameter ltering the adaptor sequence and removed the low-quality reads to achieve the clean data. UMI-tools[48] were applied for Single Cell Transcriptome Analysis to identify the cell barcode whitelist. The UMI-based clean data was mapped to human genome (Ensemble version 91) utilizing STAR [49] mapping with customized parameter from UMItools standard pipeline to obtain the UMIs counts of each sample. Cells contained over 200 expressed genes and mitochondria UMI rate below 20% passed the cell quality ltering and mitochondria genes were removed in the expression table. Seurat package [22] was used for cell normalization and regression based on the expression table according to the UMI counts of each sample and percent of mitochondria rate to obtain the scaled data. PCA was constructed based on the scaled data with top 2000 high variable genes and top 10 principals were used for tSNE construction and UMAP construction. We calculated cell cycle score using cell cycle gene lists from Tirosh et al. [50] and included this score in the normalization.
We applied graph-based Louvain cluster method with resolution = 0.8, we acquired the unsupervised cell cluster result based the PCA top 10 principal and we calculated the marker genes by FindAllMarkers function with wilcox rank sum test algorithm under default criteria. We also applied the same function to nd genes signi cantly differed between osteoclast, macrophage and monocyte. For T cell subtypes, we reran FindVariableGenes and PCA, reran clustering analysis with top 13 PCA and resolution=0.6. For SOX9 + cells, osteoblasts and Neutrophils, we also carried out similar second-level clustering analysis.

Pseudotime analysis
We applied the Single-Cell Trajectories analysis utilizing Monocle2 [51] (http://cole-trapnelllab.github.io/monocle-release) using DDR-Tree and default parameter. For osteoblast pseudotime analysis, we rst obtained a list of high variation genes by FindVariableGenes function in Seurat, then applied differential expression analysis using differentialGeneTest function in Monocle to nd high variation genes that differentially expressed between immature and mature osteoblasts. Single-Cell Trajectories analysis was applied on these differential expression genes. For APSC, we used differentialGeneTest to nd high variation genes that highly correlated with cell cycle scores and applied trajectory analysis on them.
After we calculated the pseudotime for each cell, we applied differentialGeneTest function to nd genes that signi cantly altered along pseudotime. According to the trends of alteration (i.e., sign of Pearson Correlation Coe cient between gene expression and pseudotime), we separated these genes into ascending and descending genes. We then applied Gene Ontology Biological Process (GO-BP) analysis by clusterPro ler [52] R package separately on these genes to elucidate their functions.

NicheNet analysis
Having identi ed the key gene lists in osteoblast maturation, APSC self-renewal and osteoclast transformation, we applied NicheNet [40] analysis to nd the regulation network upstream of these gene sets. NicheNet has constructed a priori networks consisting of ligands, receptors and targets. Given a set of targets and the range of expressed ligands and receptors, NicheNet nds the ligands and receptors showing the highest regulation potentiality on them. In the current study, NicheNet was applied with the following parameter: threshold of expression =25% in receptor cell and 10% in sender cell, number of ligand-target pairs=100, ligand-target activity threshold=0. 33. For uniformity, we reported top 20 ligands for all analysis.
For each of the highlighted ligand and receptor, we calculated their average expression in the corresponding cell clusters. We denoted the origin of each ligand as the cell cluster showing the highest average expression. For the ease of visualization, we aligned each target to only one receptor with the highest regulation potentiality in the a priori network. All targets that were predicted to be regulated by at least one of the ligands, but did not had potential upstream receptors, were labelled "other receptor" in the circus plot.

SCENIC analysis
To assess transcription factor regulation strength, we applied the Single-cell regulatory network inference and clustering (pySCENIC, v0.9.5) [41] work ow, using the 20-thousand motifs database for RcisTarget and GRNboost. A regulation score>1 calculated by SCENIC was taken as evidence of the regulation activity of the corresponding TF and target genes. To select TF of interest, we only included those TF downstream of the NicheNet ligands or receptors, as denoted by the a priori network of NicheNet.

Statistical analysis
Statistical analysis was carried out in R 4.0 (R Core Team). All p values were two-tailed unless otherwise speci ed. For the comparison of gene expression between cell clusters, we applied permutation test by coin R package [53]. For the enrichment of TF targets in speci c gene lists, we applied Fisher exact test with background gene list de ned as all genes with regulation score>1.

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 (calci cation of cusp tips without coalescence of other calci cations). 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 quanti cation 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 identi ed SPARC + RUNX2 + osteoblast [23], ACP5 + osteoclast [24], RGS5 + pericytes and VWF + endothelium. Lastly, we identi ed 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 immuno uorescence 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 colocalization 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 classi cation 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 subclusters ( Fig. 2A-D). We rst de ned 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 de ned 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, re ected 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-speci c intercellular signal transduction from T cell to other tooth germ cells. By summarizing the ligand-receptor pairs reaching signi cant threshold (Method), all T cell subpopulations showed strongest association with ameloblast (Number of signi cant 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 signi cant communication with ameloblast, and they exhibited subpopulation-speci c pathways like CXCR6-to-CXCL16 and CRTAM-to-CADM1. T cell also exhibited strong communication with osteoclast (Number of signi cant 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 signi cant 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 rst applied clustering analysis to dissect neutrophil subpopulation. As shown in Fig. 3A-C, we obtained eight subpopulations of neutrophil. They were classi ed by unique expression of genes related to neutrophil functions (Fig. 3C). We rst identi ed PGLYRP1 + neutrophil with speci c 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 signi cant 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 signi cant pairs = 10) and showed subpopulation-speci c 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 signi cant 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 signi cant differential expression along pseudotime, we found a clear temporal cascade  Fig. 4B and Table S3), etc. These results con rmed that the pseudotime analysis successfully reconstruct the process of osteoblast differentiation, and we therefore highlighted all genes signi cantly 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 ligandtarget 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 broblast 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 identi ed 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-speci c 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 rst 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 identi ed 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 de ned 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 nd the signaling pathways that regulate this transformation, we rst identi ed transformation-related genes by differential expression analysis (Fig. 6A). At the signi cance 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 in ammation. Since CellPhoneDB [33] analysis (Fig. 3D) found that neutrophil subtypes had dense connection with macrophage but few interactions with monocyte, we studied the neutrophil-tomacrophage signaling that could regulate the transformation-related genes. As shown in Fig. 6C, we found ve ligand-receptor pathways that were signi cantly activated between neutrophil and macrophage. Among them, TNFSF14-LTBR pathway was shared by all seven neutrophil subtypes, whereas CEACAM1-CD209 was speci c 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 signi cantly 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 signi cantly enriched in functions related to osteolysis, such as phagosome acidi cation (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 signi cant 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 signi cant functional enrichment.

Discussion
In the current study, we characterized the single-cell transcriptome of human tooth germ to decipher the cell subtype-speci c signaling pathways that regulate the biological process of tooth development. We deciphered the subtypes of resident immune cells, re ned the network of known tooth development regulators like BMP, FGF and MSX1, and discovered novel signaling pathways like.
The role of BMP family in the development and regulation of dental cells has long been highlighted by researchers. BMP family encodes various bone morphogenetic proteins, which consist of large subdivision of transforming growth factor-β ligand family [9]. In skeletal tissue, BMP regulates the osteoblastogenesis and extra-cellular matrix formation, whereas in dental tissue, BMP also regulate functions of dental pulp cells [10] and osteoclasts [14]. Following these observations, our hypothesis-free signaling network analysis further discovered that downstream pathways of these regulation were distinct. For the maturation of osteoblast, BMP regulated the osteoblast expression of ID1 and VCAN via BMP2-BMPR1A signaling, in line with their known activities during osteoblastogenesis [9]. In the selfrenewal of APSC, BMP4 and BMP5 secreted from osteoblast activated FGFR1 and downstream MSX1 to regulate a large number of renewal-related genes like SOX9 and ID3. concordantly, knock-down study on mice have demonstrated that BMP4 and MSX1 are essential in tooth organogenesis [11]. In the transformation of osteoclast from monocyte and macrophage, the role of BMP is less signi cant than major regulators CSF1, TNF and TNFRSF11, but the BMP4/2-BMPR2 pathway still showed regulation potentiality on osteoclast genes like SPP1 and GREM1. This result suggested that while BMP family took part in various biological process of tooth, the mechanism of each process is distinct and should be analyzed separately. Concordantly, each member of BMP family also showed distinct roles in different process of dental development, and their functions might be from multi-aspect. For example, BMP7 secreted from APSC could serve as ligand to act on osteoblast, whereas BMP7 expressed in osteoblast could also serve as transcription factors and regulate osteoblast maturation.
We also identi ed the complex regulation network initiated by TGF and FGF, in accordance with their known roles in tooth development. Aside from BMP family, the transforming growth factor-β (TGFB) ligand family includes various genes including TGFB1, TGFB2 and TGFB3 [13]. In dental tissue, TGFB ligands regulate the pulpal repair and dentinogenesis, possibly through the SMAD2 and ERK pathways in pulp cell [12]. In our single-cell analysis, we further found that TGFB1, which was mainly originated from T cells, activated receptor TGFBR1 and TGFBR2 to regulate a large number of genes involved in osteoblast maturation. This regulation was in part mediated by transcription factor FOXQ1. As for APSC renewal, both TGFB1 and TGFB3 (mainly originated from osteoblast) showed high regulation potentiality.
Similarly, we also found the cell-type speci c networks of FGF signaling, which was known to play a role in tooth development but the intact pathway remained to be elucidated [45]. Speci cally, receptor FGFR3 was involved in osteoblast maturation, whereas FGFR1 was mainly responsible for ASPC self-renewal.
The downstream transcription factor (ID1 and MSX1, respectively) was also distinct in these two processes. Another interesting fact is that aside from FGF2, ligands BMP2, 4 and 5 also showed a nity with FGF receptor and exhibited even higher regulation potentiality, highlighting the importance of crosspathway signaling transduction.
Our analysis also discovered the role of immune cell in dental development, which was not strengthened by previous study. For example, IL1B, which was mainly expressed in monocyte, regulated 34 APSC renewal-related genes, more than other non-immune ligands could regulate. It also regulated the maturation of osteoblast together with IL1A. It is known that lymphocytes could inhibit the dental pulp development by secreting cytokines like IL1B and IL6, but this inhibition was only found in in ammatory status [46]. Alternatively, since the tooth germ sample in the current study was collected at normal status, our result suggested that IL1B could also regulate the normal developmental process of human tooth in the absence of in ammation. Similarly, TGFB1 also regulated 39 renewal-related genes, and was mainly expressed in T cells. As for osteoclast transformation, the top ligand was TNF from macrophage. These results indicated that dental immune cells not only defense against pathogens, but also regulate the dental development via secretion of ligands that act on other dental cell types. Under the situation of in ammation or stress, such regulation might lean towards suppression of osteogenesis and activation of osteoporosis.
It should be noted that since dental germ developmental stages exhibit high heterogeneity, and that samples collected from different individual may represent distinct status and could not cover the entire process of dental development. For example, the samples analyzed in the current study contained relatively small number of ameloblasts and odontoblasts, which are mainly activated at the early stage of dental development. Since it is relatively di cult to obtain highly immature dental tissue in clinical scenario, the signaling networks of these cell types might be better analyzed using prenatal tissues or by translational studies of mouse model. In the future, more efforts should be devoted to the single-cell analysis of human embryo dental tissues.
In conclusion, we provided a cell interactome landscape for postnatal human germ and discovered the key signaling pathways regulating the development of dental cells, which provided novel insights into the mechanism of dental development and highlighted potential targets for disease intervention and dental regeneration.

Declarations
Ethics approval and consent to participate This study was approved and supervised by Ethical committee of Shanghai Ninth People's Hospital.
Written informed consent was provided by all participants.
Availability of data and materials Expression data will be made available at https://github.com/WeiCSong upon publication.