2.6.10.Statistical analyses
The version of R for statistical analysis and visualization was version 4.1.0. The The statistical methods used for protein differential analyses included one-way analysis of variance (ANOVA) and students’ t-test and the correlation analysis method was Spearman’s correlation coefficient.
Results
The quantitative protein profile of human deciduous molar germ from cap-to-early bell stage
During the cap stage (after PCW 10), the size and shape of the tooth crown become apparent[13]. To determine whether the gestational age was consistent with the developmental stage of tooth germ, H&E staining was performed. At PCW12 stage, the enamel organ was cap-like structure, the primary enamel knot was clearly visible. (Figure 1A). At PCW15 stage, the area of dental papilla (DP) significantly enlarged, the shape of the enamel organ changes (Figure 1B). When tooth germ developed to PCW18, the odontoblast cell layer could be visualized and the dentin cells begin to secrete a matrix, which differentiates into dentin in the future (Figure 1C).
To identify the protein compositions of tooth development from cap-to-early bell stage, 15 samples, grouped into PCW12, PCW15, and PCW18 (five samples per group), were examined for label-free mass-spectrometry proteomics. A total of 6327 proteins were determined in human deciduous molar germ after filtration by FDR <1% and containing at least one unique peptide. To ensure the data repeatability, unsupervised hierarchical clustering, Pearson’s correlation coefficient, and PCA analyses were performed (Figure S1A-C). The protein expression patterns of 15 samples showed a temporal partition of data points. Furthermore, only proteins detected in at least 3 replicates in each stage were selected for further analyses. After data filtrating, 6084 proteins were finally adopted. The overlapping proteins present in the three stages were summarized (Figure 2A). The CV-value of each protein for each stage was calculated, and the median CV for each stage was compared. The intra-group variability of the overlapping proteins exhibited a time-dependent increase (Figures 2B and 2C). Interestingly, the dental cell markers such as VCAN and COL1A1 showed high heterogeneity during tooth development, while the classical housekeeping proteins (GAPDH and TUBB2B) did not change significantly (Figure S1D).
Next, the PANTHER Classification System (version 17.0) was applied to classify the overlapping proteins (Figure 2D). The top three categories of proteins detected were metabolite interconversion enzyme (20.8%), protein modifying enzyme (13.3%), and RNA metabolism proteins (9.3%). Proteins that were expressed in only one stage were defined as stage-specific proteins (Table S1). Among the stage-specific proteins, 149 of PCW12 were involved in histone deacetylation, mitochondrial translation, and chromatin remodeling; 29 stage-specific proteins in PCW15 were involved in the establishment of skin barrier; PCW18 stage-specific proteins were involved in stem cell proliferation.
To identify the molar tissue-specific proteins, 6084 proteins in our deciduous molar proteome were compared with the proteins in Human Protein Atlas (HPA) and GTEx databases (Figure 2E). Comparison with HPA led to the identification of 450 tissue-specific proteins (Table S2); tissue-specific proteins mainly function in ossification and collagen metabolic processes. Comparison with GTEx led to the identification of 20 tissue-specific proteins. Cross-comparison revealed WASH2P and C11orf96 as unique proteins in the deciduous molar proteome, but their functions need to be further elucidated.
Next, the scores of protein signatures in representative signaling pathways were calculated using the ssGSEA algorithm (Figure 2F). A detailed list of proteins for signaling pathways is provided in Table S3. Pathways involved in epithelial cell development, mesenchyme development, osteoblast development, and tooth mineralization gradually increased with developmental stage. Besides, protein signatures involved in WNT and BMP signaling pathways also increased significantly.
The temporal expression patterns of human deciduous molar proteome
The overall temporal expression patterns of all proteins were firstly analyzed by fuzzy c-means algorithm. The numbers of clusters were determined from the elbow plot (Figure S2A). A total of five clusters of temporal patterns were summarized (Figure S3, left). GO analysis of proteins from C3 and C5 revealed that the up-regulated proteins had functions, including hydrogen peroxide catabolism, glycolysis, negative regulation of endopeptidase activity, and integrin-mediated signaling pathway. These proteins included more than 50% of detected extracellular matrix (ECM) related proteins (Table S4).
The down-regulated proteins functioned mainly in translation, DNA replication, intracellular protein transport, and mitotic cell cycle. These proteins mainly included factors that play critical roles in regulating cell proliferation, such as telomere replication-related proteins, NOP2 and MCM4 (Table S5). Meanwhile, the bipolar proteins were associated with the biological process of cellular catabolic and endomembrane system organization. These proteins included defined epithelial cell markers KRT14, KRT15, SPRR1B, CAVIN1, etc. All identified molar transcription factors (TFs) have been listed to comprehensively understand their expression pattern (Figure S3, right).
Next, the co-expression network was constructed using 6084 filtered proteins by WGCNA analysis. Sixteen modules were obtained through hierarchical clustering and dynamic tree clipping (Figure 3A). Among these modules, the MEpaleturquoise and MEgreenyellow modules were highly and positively associated with the fetal ultrasonic age (Figures 3B and 3C). In total of 1368 proteins in these two modules were selected. These proteins executed functions of protein transport and catabolic and metabolic processes (Figure 3D). In addition, the GLI family involving Hedgehog signaling transduction and VEGFA-VEGFR2 signaling pathway were significantly enriched. After ranking the connectivity of selected modules, 24 hub proteins were obtained (Figure 3E). U2AF2, HMGB2, and PPP1R12A are associated with embryonic stem cell proliferation and differentiation[14-16]. WLS, the key regulator for Wnt trafficking, was continuously up-regulated from the cap-to-bell stage[17].
In order to explore cell origins of the hub proteins, the single cell transcriptomic data from published GEO datasets were analyzed. After data quality control, the expression data for 23530 single cells were obtained (Figure 3F). Epithelial origin cells were defined based on dental epithelial markers KRT14, SPRR3, PITX2, SCUBE3 andMPPED1 . Mesenchymal origin cells were characterized by the expression of PAX9, VIM, FN1, VCAN and PRRX1 . Endothelial origin cells were defined based on CDH5, EMCN, PECAM1, RGS5 andMCAM . Schwann cells were defined based on PLP1, SOX10 andCDH19 . Immune cells were characterized by the expression ofPTPRC, CD4, CD163 and CD86 . Similar to the previous literature, the main cell types of embryonic deciduous tooth germ were dental mesenchymal and epithelial cells. And a relative proportion of endothelial cells, a small number of immune cells and nerve cells were also identified.
Then, the relative expression of hub protein-encoding genes was compared among different cell populations. Since the number of Schwann cells was too small, this cell population was excluded from analyses. As shown in Figure 3G, the epithelial cells expressed higher level of WLS, RBP1, CLSTN1, GSS, ANP32A and the mesenchymal cells expressed higher levels of SPON1, DDAH1, HMGB1, GPALPP1, FREM1, PRPS2, RAP1GDS1,and AKAP12 . The mRNA expression of WLS in dental epithelial cells and the expression of AKAP12, RAP1GDS1, DDAH1, SPON1 and FREM1 in dental mesenchymal cells were also gradually increased from PCW12 to PCW18. The hub protein-encoding genes that were highly expressed in each cell populations showed similar up-regulated temporal expression pattern. These results indicated that the hub proteins screened by proteomics WGCNA analysis were worthy of further exploration. It should be emphasized that HMGB1 may play important role in osteogenic development. It has been recently reported that variants in HMGB1 could cause brachyphalangy, polydactyly and tibial aplasia syndrome, a rare complex malformation syndrome[18].