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].