Erf Affects Commitment and Differentiation of Osteoprogenitor Cells in Cranial Sutures via the Retinoic Acid Pathway

ETS2 repressor factor (ERF) haploinsufficiency causes late-onset craniosynostosis (CRS) (OMIM entry 600775; CRS4) in humans, while in mice Erf insufficiency also leads to a similar multisuture synostosis phenotype preceded by mildly reduced calvarium ossification. However, neither the cell types affected nor the effects per se have been identified so far. ABSTRACT ETS2 repressor factor (ERF) haploinsufficiency causes late-onset craniosynostosis (CRS) (OMIM entry 600775; CRS4) in humans, while in mice Erf insufficiency also leads to a similar multisuture synostosis phenotype preceded by mildly reduced calvarium ossification. However, neither the cell types affected nor the effects per se have been identified so far. Here, we establish an ex vivo system for the expansion of suture-derived mesenchymal stem and progenitor cells (sdMSCs) and analyze the role of Erf levels in their differentiation. Cellular data suggest that Erf insufficiency specifically decreases osteogenic differentiation of sdMSCs, resulting in the initially delayed mineralization of the calvarium. Transcriptome analysis indicates that Erf is required for efficient osteogenic lineage commitment of sdMSCs. Elevated retinoic acid catabolism due to increased levels of the cytochrome P450 superfamily member Cyp26b1 as a result of decreased Erf levels appears to be the underlying mechanism leading to defective differentiation. Exogenous addition of retinoic acid can rescue the osteogenic differentiation defect, suggesting that Erf affects cranial bone mineralization during skull development through retinoic acid gradient regulation.

C ranial sutures comprise the connective tissues between the bones of the skull and are considered to be major centers of bone growth during development (1,2). Mesenchymal stem cells that reside in the suture mesenchyme commit and enter the intramembranous ossification pathway, giving rise to proliferating populations of osteoprogenitor and preosteoblast cells that eventually differentiate to osteoblasts at the edges of the developing bones (3,4). The balance between the different cell types seems to be crucial for suture patency and consequently for the coordination of skull and brain development (5,6). Craniosynostosis is a developmental disorder in which one or more of the cranial sutures close prematurely, leading to skull and facial deformities, along with complications that can affect vision, hearing, breathing, and learning ability. During the last 25 years, considerable effort has been put into unraveling the mechanisms that underlie craniosynostosis (7). However, the presence of multiple cell populations in sutures, the paucity of specific cellular markers, and difficulties in the identification and isolation of suture stem cells have hindered these efforts.
Genetic analysis has identified an increasing number of genes that, when mutated, cause craniosynostosis. Activating mutations in three members of the fibroblast growth factor receptor family (FGFR1, FGFR2, and FGFR3) and alterations in downstream signaling cascades such as the p38 mitogen-activated protein kinase (MAPK), extracellular signal-regulated kinase 1/2 (ERK1/2), phosphatidylinositol 3-kinase (PI3K)/AKT, and phospholipase Cg (PLCg)/protein kinase C (PKC) pathways have been commonly reported to be involved in syndromic

RESULTS
LIF-selected long-term expanded suture derived cells possess in vitro characteristics of mesenchymal stem/progenitor cells. Cranial sutures constitute niches of highly heterogeneous cell populations related to bone growth (37). We thus focused our efforts on mesenchymal stem cell (MSC)-derived populations and, based on previous studies, established a new protocol utilizing leukemia inhibitory factor (LIF) for the selective expansion and maintenance of mesenchymal stem/progenitor cells from cranial sutures. Suture explants from postnatal day 5 (P5) mice and the resulting suture-derived cells were cultured in the presence of leukemia inhibitory factor, which is known for its role in sustaining the stem cell state while inhibiting differentiation (41,42). Cultivation of suturederived cells in the presence of LIF for a minimum of 8 population doublings (PDs) during a period of 50 to 60 days resulted in a population of cells that were plastic adherent, fibroblast-like in shape ( Fig. 1A and B), and expressed increased levels of the MSC marker Axin2 (43) and reduced levels of the osteogenic differentiation marker Sp7 compared to those in the initial population (Fig. 1C). The majority of these cells expressed the MSC-associated surface antigens CD44, CD90, CD29, and Sca1 (4,(44)(45)(46), while neither hematopoietic nor endothelial cell markers could be detected (Fig. 1D). This cell population growing in culture for more than 8 PDs could effectively undergo differentiation toward the chondrogenic, osteogenic, and adipogenic lineages (Fig. 1E), a hallmark of mesenchymal stem cells. These cells, which we label suture-derived mesenchymal stem/progenitor cells (sdMSCs), can be routinely maintained in culture for more than 20 PDs (Fig. 1F) and sustain their characteristics for at least 3 freeze-thaw cycles.
Utilizing this approach, we established sdMSCs from Erf loxP/1 and Erf loxP/2 P5 littermates, in at least 5 independent experiments, to study the effect of limited Erf levels on MSC growth and differentiation. At this time point, the mice have not yet developed the phenotype of synostosis.
Erf insufficiency compromises the commitment of suture mesenchymal stem/ progenitor cells toward the osteogenic lineage. Although Erf is known to affect cellular proliferation (16,47), cell cycle phase analysis of Erf loxP/1 and Erf loxP/2 sdMSCs showed no significant difference in the cell distribution profiles ( Fig. 2A). There was also no difference in the cell doubling time throughout the life of the cultures (Fig. 2B), suggesting that Erf insufficiency does not affect sdMSC self-renewal rate. We then (C to E) sdMSCs were induced to differentiate along the chondrogenic lineage for 21 days (C), the adipogenic lineage for 21 days (D), and the osteogenic lineage for 28 days (E) and stained with alcian blue and hematoxylin, oil red O, and alizarin red S, respectively. Bars, 10 mm, 50 mm, and 100 mm, respectively. (F) Measurements of the alizarin red S dye extracted from the cells after 28 days of osteogenic differentiation. Three independent biological experiments were conducted, and the statistical analysis was performed using an unpaired t test with two-tailed distribution. *, P , 0.05. examined the impact of Erf levels on sdMSC differentiation. Erf loxP/1 and Erf loxP/2 cells showed comparable efficiency in in vitro chondrogenic and adipogenic commitment ( Fig. 2C and D). However, Erf loxP/2 cells displayed decreased ability to mineralize ( Fig. 2E and F), implying an impairment in the osteogenic differentiation of these cells. The reduced osteogenic differentiation was also apparent in the initial heterogeneous suture-derived cell population, in which Erf loxP/2 cells displayed initially comparable but later decreased capacity to mineralize ( Fig. 3A and B). Chondrogenic differentiation appeared unaffected, while the limited adipogenic differentiation was increased in Erf loxP/2 cells (Fig. 3C to E), indicating the possible existence of a greater proportion of precursor populations among Erf-insufficient suture cells. These results indicate that the insufficiency of Erf reduces the osteogenic ability of suture-derived cells in vitro.
To gain further insight into the apparent differences in sdMSC osteogenic differentiation and to understand the possible role of Erf levels in this process, we analyzed transcriptome changes in suture-derived cells from Erf loxP/1 and Erf loxP/2 mice in 3 conditions. sdMSCs at PD8 in self-renewal or osteogenic differentiation medium for 24 h were used for early differentiation events, while the initial heterogeneous suture cell population growing in expansion medium, without LIF was used as a mixed differentiation state specimen (Fig. 4A). Consistent with a decrease rather than the elimination of Erf expression, a very limited number of genes were found to differ between the two genotypes. This was evident by the clustering of transcriptome samples according to culture conditions but not genotypes ( Fig. 4B and Table 1). However, all 3 growth conditions produced distinct patterns that clustered together, assuring the consistency of the ex vivo cultures (Fig. 4C). The limited number of differentially expressed genes indicated a deficit of core matrisome genes in Erf loxP/2 cells when tested against the GSEA database (48) (Fig. 5A and B and Table 1), supporting the hypothesis that Erf insufficiency leads to a defect in the osteogenic differentiation. This was consistent with the  ). At least 3 or 4 mice were used per genotype for each of the above conditions tested per experiment, and at least 4 independent experiments were conducted. (B) Unsupervised clustering of gene expression experiment from Erfcompetent Erf loxp/1 (P/1) and Erf-insufficient Erf loxp/2 (P/2) cells, indicating no genotype-specific associations. Self-renewal, sdMSCs in self renewal medium; osteo differentiation, sdMSCs 24 h in osteogenic differentiation medium; fresh cells, freshly isolated suture-derived cells. (C) Heatmaps based on the 60 most differentially expressed genes of the indicated comparison, showing a clear clustering between conditions but not between genotypes. P/1, Erfcompetent (Erf loxP/1 ) cells; P/2, Erf-insufficient (Erf loxP/2 ) cells. transcriptome differences of Erf-competent and Erf-insufficient cells upon osteogenic induction (see Table S1 in the supplemental material). Erf loxP/2 cells exhibited considerably fewer genes associated with ossification and extracellular matrix organization than Erf loxP/1 cells during induction of sdMSCs ( Fig. 5C, L-O_minus and L-O_plus). Consistently, Erf loxP/2 sdMSCs that either self-renewed or differentiated for 24 h required many more ossification-related changes to reach the differentiation state of the initial heterogeneous cell population than the Erf loxP/1 cells (Fig. 5C, O-F_minus and O-F_plus).
We further examined the apparent contribution of Erf expression in the effective osteogenic differentiation, interrogating single cell RNA-sequencing data from mouse sutures available through the FaceBase Consortium (49,50). Given the ubiquitous expression of Erf and its posttranscriptional regulation, we developed an approach to examine gene coexpression rather than cell cluster expression. We determined any expression correlation between the gene of interest and the rest of the cellular transcripts for each informative cell in the data and evaluated the known function of the correlated genes. Erf expression appeared to correlate with genes involved in ossification and extracellular matrix organization ( Fig. 6A; see also Table S2 in the supplemental material). Compared to other suture ossification landmark genes subjected to the same correlation analysis, Erf clustered closely with Sp7 in cells at E16.5 and with Fgfr1, Runx2, Twist1, and Alpl in cells at E18.5 and P10 ( Fig. 6B and Table S2).
These data suggest that appropriate Erf expression level is required for proper differentiation of cranial suture cells toward the osteogenic pathway and are consistent with the decreased mineralization pattern observed previously in vivo (20), which could account for the late onset of Erf-related synostosis phenotype.
Erf insufficiency-induced osteogenesis defect can be rescued by retinoic acid. In spite of the limited number of genes found to differ between Erf-competent and Erf-insufficient cells in all growth conditions, a group of genes associated with the retinoic acid (RA) pathway could be identified (Table 1). Characteristically, Cyp26b1, a gene coding for an RAcatabolizing enzyme known to affect suture development leading to craniosynostosis (32), was elevated upon Erf insufficiency in both proliferating and differentiating sdMSCs and in the initial heterogeneous suture cell population (Table 1 and Fig. 7A). Cyp26b1 was drastically reduced upon normal sdMSC differentiation but remained in relatively high levels in Erf-insufficient cells (Fig. 7B), suggesting decreased levels of retinoic acid. Comparison of the genes differentially expressed during the differentiation of Erf-competent and Erf-insufficient sdMSCs with genes regulated by retinoic acid in mouse embryonic stem cells (mESCs) (51, 52) revealed a 4-to 5-order-of-magnitude-higher significance for the Erf-competent (Erf loxP/1 ) cell gene set ( Fig. 7C; see also Tables S3 and S4 in the supplemental material). Analysis of these common genes (Table S4) through Metascape (53) indicated that indeed, in the case of Erf-competent (Erf loxP/1 ) cells, they are associated with genes known to play a   Table S1 in the supplemental material) of sdMSCs (L), sdMSCs growing in osteogenic differentiation medium for 24 h ("O"), and freshly isolated suture cells ("F"), from Erf-competent (Erf loxP/1 ; "plus") and Erf-insufficient (Erf loxP/2 ; "minus") animals, were clustered, based on their ontology, with the Metascape program. L-O_comm, genes found to be differentially expressed during the 24-h induction of sdMSCs of both Erf-competent and Erf-insufficient cells; L-O_plus, genes found to be (Continued on next page) major role is suture development, such as Runx, Twist, and BMP-Smad. In contrast, in the case of Erf-insufficient (Erf loxP/2 ) cells, they appear to associate with lymphoid and endothelial pathways (Fig. 7D). It would thus appear that Erf deficiency leads to a decreased retinoic acid response that affects suture development.   Table S2 in the supplemental material) to correlate with Erf expression were clustered based on their ontology with the Metascape program. P10, single-cell RNA sequencing (scRNA-seq) data from P10 animals; E16.5, scRNA-seq data from E16.5 embryos; E18.5, scRNA-seq data from E18.5 embryos. (B) As in panel A, but with the addition of data from genes correlated with Alpl, Fgfr2, Runx2, Sp7, and Twist1 as markers for osteogenic differentiation stages. Red and orange lettering indicate ossification-related and differentiation-related ontologies, respectively. Only the top 20 categories are shown. Retinoic acid affects multiple developmental processes and pathways, and it has been suggested that its homeostasis is crucial for normal skeletogenesis (33,54,55). We thus evaluated the effect of RA on Erf-competent and Erf-insufficient sdMSCs. A characteristic feature of the differentiating suture-related Erf loxP/2 sdMSCs is the higher initial proliferation and final cell numbers ( Fig. 7E and F), consistent with their decreased capacity to exit self-renewal and commit. Addition of retinoic acid at low concentrations does not appear to affect the growth/survival of Erf-competent sdMSCs (Fig. 7G). However, RA addition fully suppressed the increased cell numbers of the differentiating Erf-insufficient cells (Fig. 7G). More importantly, the decreased mineralization potential of Erf loxP/2 cells was fully alleviated in the presence of RA without affecting the potential of the Erf-competent cells (Fig. 7H).
These data strongly suggest that Erf deficiency decreases retinoic acid levels leading to increased cellular proliferation and decreased osteogenic differentiation. Such changes could be the underlying cause of the late onset Erf-related craniosynostosis phenotype.

DISCUSSION
Syndromic craniosynostosis due to ERF haploinsufficiency presents some unique challenges and opportunities for disease understanding and management. In contrast to FGFR/MAPK-driven craniosynostosis syndromes, it has a late-onset phenotype, variable severity, and, in the mouse model, an initially decreased calvarial ossification. It would thus appear that Erf could mediate effects at multiple stages during suture development.
Understanding the formation of the cranial sutures is a challenging problem that is hampered by the multiple origins of the involved cells. We thus established a reliable and reproducible system to derive mesenchymal stem cells from murine cranial sutures and address the contribution of Erf levels in the process. Our cellular data indicate that although Erf elimination can affect multiple developmental processes, Erf insufficiency specifically attenuates osteogenic differentiation. This would be consistent with the role of Erk1/2 in sustaining the undifferentiated state of mesenchymal stem cells (56) and conforms with the absence of other major developmental defects in Erf 1/2 and Erf loxP/2 mice (17,18,20). The transcriptional changes directly associated with Erf insufficiency were fairly limited, probably because the 50% reduction in Erf may not be sufficient to induce extensive quantitative changes that would evade our detection cutoff limits. However, the expression program during the differentiation of the sdMSCs was altered extensively as a result of the decreased Erf levels, suggesting a homeostatic effect or an effect on one or more morphogenetic factors.
One consistent effect of Erf insufficiency through all biological replicates and culture conditions was the elevated expression of Cyp26b1, a major retinoic acid catabolizing-enzyme, suggesting a possible decrease in retinoic acid levels affecting bone development (32-35, 54, 55, 57, 58). Consistently, addition of RA during differentiation of sdMSCs fully alleviated the Erf deficiency defect, with no apparent effect on Erf-competent cells, suggesting that Erf, through Cyp26b1 regulation, affects retinoic acid levels, thus controlling the differentiation outcome. The inhibitory effect of Erf on Cyp26b1 and the concomitant increase in retinoic acid levels would be consistent with the role of FGF, which has already been reported to activate Cyp26b1 and downregulate RA availability (39,40,59). Given that Erf is inactivated via nuclear export as a result of Fgf signaling, the decreased Erf levels would resemble the increased Fgf signaling state. It is unclear at this point if Erf regulates Cyp26b1 directly or indirectly. Chromatin occupancy experiments in four independent mouse and human cellular systems, from our laboratory and others (20,(60)(61)(62), fail to identify any ERF binding site in the vicinity of Cyp26b1. The closest Erf chromatin binding site appears to be 300 to 400 kb away in the proximity of the Exoc6b gene, rather than that of Cyp26b1. Neither reports of other Ets domain transcription factors directly affecting Cyp26b1 transcription nor any confirmed functional Ets binding motifs in its regulatory regions have been found so far. Thus, Erf could affect its expression, either through a distal control element or indirectly by inhibiting a Cyp26b1 activator (Fig. 8).
It is also conceivable that Erf levels affect MSC differentiation and homeostasis in synergy or in competition with other factors. We have previously shown a close proximity of Erf with Runx binding sites and a possible orchestrated regulation by the two factors (20). Erf could antagonize other Ets proteins, such as Ets1 and Ets2, that have been reported to affect skull shape and osteogenic differentiation (63,64). Interestingly, genes differentially affected during Erf-competent and Erf-insufficient cell differentiation that harbor sites identified in Erf chromatin immunoprecipitation sequencing (ChIP-seq) experiments have these sites almost exclusively in putative regulatory regions and away from the transcription start site, where the repressive effect of Erf is more prominent (see Table S5 in the supplemental material). It would thus appear that Erf, in an orchestrated manner with other transcription factors, could regulate the homeostasis and kinetics of the process rather than being instructive.
The delayed ossification observed in MSC cultures may directly explain the early ossification defect observed in Erf-insufficient mice but may also be responsible for the craniosynostosis phenotype. Decreased MSC differentiation may lead to decreased populations of committed progenitor cells within the suture, which in turn fail to keep the suture open, resulting in craniosynostosis. Stem cell depletion (4), osteogenic differentiation-proliferation imbalance (6,65,66), and growth dynamics (67) have been previously proposed as potential mechanisms of synostosis. Indeed, there are several systems where Erf regulates differentiation. Chorionic trophoblast differentiation depends on Erf presence (17), Erf blocks the differentiation of Ras-deficient mESCs (62), and the rate of embryonic hematopoietic differentiation depends on Erf level (18). Whether Erf has additional roles, either positive or inhibitory ones, at later stages of the osteogenic differentiation pathway remains to be explored. Such an inhibitory effect in suture preosteoblast differentiation would be consistent with the reported function of Erk in late stages of the osteogenic pathway (68)(69)(70), as well as with the role of Runx2 in osteogenic differentiation and the antagonistic role of Erf (20), which could account for the craniosynostosis phenotype. Interestingly, comparable concentrations of retinoic acid have opposite effects on the mineralization of mesenchymal stem cells and more differentiated osteogenic populations such as MC3T3 preosteoblasts, with the former being induced and the latter being suppressed by RA (58,(71)(72)(73)(74)(75). Therefore, it is plausible that Erf, either through retinoic acid signaling or independently, could also affect the osteogenic pathway at later stages. Upon Erf deficiency, the initial decreased MSC differentiation and osteoprogenitor expansion could lead to decreased ossification, while an accelerated differentiation of progenitors when they comprise the predominant suture cell population could lead to depletion and premature ossification of the suture. Further studies would be necessary to explore the spatiotemporal function of Erf and its effect on RA concentration gradients in cranial bone and suture development.
In conclusion, our work provides evidence that Erf, an effector within the FGF/ERK pathway, affects suture formation and calvarial ossification through the regulation of retinoic acid levels.

MATERIALS AND METHODS
Mouse lines. Mice were bred and maintained in the animal facility of the Institute of Molecular Biology and Biotechnology (IMBB) in Greece. All experimental protocols were conducted with ethical guidelines and in compliance with the 3Rs. Protocols were approved by the bioethics committee of the IMBB and licensed from the General Directorate of Veterinary Services, Region Crete (permit numbers EL 91BIO-02 and EL91-BIOexp-02; project license no. 27289 to G. Mavrothalassitis). Erf loxP/1 and Erf loxP/2 littermates were obtained from crossing Erf 1/2 mice with Erf loxP/loxP mice, both of which have been reported in the literature (20).
Suture cell cultures and differentiation assays. Freshly isolated cranial suture cells were obtained as described previously (37,58). Briefly, sagittal and coronal sutures containing approximately 0.5-mm bony margins were dissected from skulls of 5-day-old mice (P5). Following the removal of skin and dura mater, the suture explants were placed into 100-mm cell culture dishes, cut into tiny pieces, and cultured in the presence of Dulbecco's modified Eagle's medium (DMEM)-F-12 medium (catalog no. 31330038; Gibco, Thermo Scientific) supplemented with 10% fetal bovine serum (FBS, catalog no. 10270106; Gibco, Thermo Scientific) and 100 U/ml penicillin-streptomycin (catalog no. 15140122; Gibco, Thermo Scientific) for 7 days. To obtain suturederived long-term expanded mesenchymal stem/progenitor cells, suture explants were cultured into gelatinized dishes in the presence of KnockOut DMEM (catalog no. 10829018; Gibco, Thermo Scientific) supplemented with 12% FBS (catalog no. SV30160.03; HyClone, GE Healthcare), 2 mM L-glutamine (catalog no. 25030024; Gibco, Thermo Scientific), 1% nonessential amino acids (MEM-NEAA; catalog no. 11140050; Gibco, Thermo Scientific), 100 U/ml penicillin-streptomycin, 0.1 mM b-mercaptoethanol, and 0.3% vol/vol leukemia inhibitory factor (LIF; produced in our laboratory) for 7 days. Cells were then harvested and further cultured (without the explants) under the same conditions for a total of at least 8 population doublings (PDs) during a period of approximately 50 to 60 days (Fig. 1A). The cumulative population doubling level (CPD) was calculated as described in previous studies (76). At this stage, the cells were immunophenotypically evaluated by flow cytometric analysis and checked for osteogenic, chondrogenic, and adipogenic differentiation potential.
Osteogenic and chondrogenic differentiation was performed as described previously (77,78). Osteogenesis was induced in 60 to 70% confluent cultures by the addition of DMEM (low glucose, catalog no. 21885025; Gibco, Thermo Scientific) supplemented with 10% FBS, 0.1 mM dexamethasone, 50 mM ascorbate-2-phosphate, and 10 mM b-glycerophosphate. In the experiments with retinoic acid, freshly prepared osteogenesis medium was supplemented with all-trans retinoic acid at a final concentration of 0.5 mM before every medium change. Adipogenesis was induced in postconfluence cultures by switching between adipogenic induction and adipogenic maintenance medium (79). One cycle of induction-maintenance was performed for freshly isolated suture cells and 3 cycles for LIF selection-subjected long-term expanded mesenchymal stem/progenitor cells. To assess the extent of differentiation, staining of cultures with oil red O (catalog no. O-9755; Sigma) and alcian blue (catalog no. A-5268; Sigma) was performed for the detection of adipocytes and cartilage, respectively. The evaluation of osteogenic differentiation was performed by alizarin red S (Sigma A-5533) staining of the cultures, followed by acetic acid extraction and quantification of the dye at 405 nm as already described (80).
Cell growth and viability studies. Cell doubling time was estimated at specific population doubling levels of the culture by using the already described logarithmic equation (81). The cell cycle phase study was performed in isolated nuclei by propidium iodide (PI) staining of cells in hypotonic solution (82), followed by flow cytometric analysis. Cells of population doubling level 20 (20 PDs) at 60 to 70% culture confluence were used in all cell cycle experiments. Data were analyzed using ModFitLT software. In order to evaluate the viability of cells during the osteogenic differentiation, an MTT [3-(4,5-dimethyl-2-thiazolyl)-2,5-diphenyl-2H-tetrazolium bromide] assay was conducted, and formazan absorbance was measured at 600 nm (83).
For in vitro proliferation assays, bromodeoxyuridine (BrdU; catalog no. B5002; Sigma) was added to the cell culture medium at a final concentration of 10mM for 8 h, followed by fixation of the cells with 4% paraformaldehyde (PFA) solution. The detection of BrdU-positive cells was performed using the following antibodies and reagents: rat anti-BrdU antibody (catalog no. MCA2060GA; AbD Serotec) at a dilution of 1:800, biotin-conjugated anti-rat antibody (catalog no. B7139; Sigma) at 1:100, and fluorescein isothiocyanate (FITC)-conjugated streptavidin (catalog no. 405201; BioLegend) at 1:1,000. A TCS SP2 confocal microscope (Leica Microsystems) and an Operetta imaging system were used for signal visualization and analysis.
Flow cytometric analysis. LIF selection-subjected mesenchymal stem/progenitor cells of 8 PDs were harvested using 0.25% trypsin-EDTA solution (catalog no. 25200072; Gibco, Thermo Scientific) for 2 min at 37°C and stained with the following antibodies in 1% FBS-phosphate-buffered saline (PBS) solution for 30  Quantitative PCR. Total RNA was extracted from cultured suture cells using TRI reagent (catalog no. T9424; Sigma) according to company's instructions. The mRNA was reverse transcribed using the SuperScript first-strand synthesis kit (catalog no. 11904-018), and 5 ng of the total synthesized cDNA was added in each real-time qPCR using 2Â Brilliant III SYBR green quantitative PCR (qPCR) master mix (catalog no. 600882-51; Agilent) in an Applied Biosystems StepOne Plus real-time PCR machine. The expression levels of the following genes were detected using the following sets of primers: Axin2 FW, 59-AGCCTAAAGGTCTTATGTGG-39, and RV, 59-ATGGAATCGTCGGTCAGT-39; Osterix (Sp7) FW, 59-TCTGCTTGAGGAAGAAGCTC-39, and RV, 59-TCCATTGGTGC TTGAGAAGG-39; and Gapdh FW, 59-CCAGTATGACTCCACTCACG-39, and RV, 59-GACTCCACGACATACTCAGC-39. The expression levels of the genes of interest were normalized to Gapdh expression levels for each particular sample.
RNA sequencing. Total RNA was isolated using the Qiagen RNeasy minikit. Each biological replicate was created by pooling suture-derived cells of at least three mice. Three or four independent biological replicates were conducted for the conditions tested. Next-generation sequencing (NGS) libraries were generated from 500 ng input total RNA with the Lexogen-QuantSeq 39 mRNA-Seq library prep kit FWD for Illumina and run on an Illumina 500 instrument on 1 Â 150 FlowCells.
Single-cell correlation analysis. Count matrices of single-cell RNA sequencing (scRNA-seq) data were first filtered following the quality assessment suggested by Harvard Chan Bioinformatics Core (https://hbctraining.github.io/scRNA-seq/lessons/04_SC_quality_control.html) and normalized following Seurat's default method (88). Features that were not detected in at least 2% of the cells were also eliminated to improve reliability of a possible correlation. Gene correlations with the false discovery rate at 0.05 significance were calculated using the "corr.test function" (89) in the R statistical environment (90). The Wilcoxon rank sum test, as implemented in the "wilcox.test" function from the stats package (90), was used to further evaluate differences in the distribution of the correlated gene in cells expressing the target gene or not. Enrichment analysis sets for Mus musculus were performed with the gprofiler2 package (91), with a statistical domain size comprising genes that have at least one annotation and with the g:SCS multiple testing correction method. The whole workflow was implemented in R version 3.6.1 (5 July 2019). Clustering of correlated gene sets across different scRNA data sets and target genes was visualized with the Metascape web tool (53) and ComplexHeatmap (92). The analysis software and genelist files are deposited at https://github.com/mpaltsai/iRNA_project. Statistical analysis. All experiments had a minimum of two biological replicates and two experimental replicates for each. At least 3 or 4 sibling mice were used to derive each biological sdMSC sample. Variation among experimental replicates was minimal. Equal numbers of experimental replicates were used for each biological replicate in each experiment. The data, unless otherwise stated, were analyzed using SPSS software. An unpaired (two-sided) t test was conducted for comparisons between two groups, and one-way analysis of variance (ANOVA) was performed for multiple comparisons, followed by a post hoc Dunnett's two-sided test according to the experimental requirements.
Data availability. Sequencing data are deposited in the BioProject database under accession number PRJNA664970.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, XLSX file, 0.2 MB.