Lineage-specific differentiation of osteogenic progenitors from pluripotent stem cells reveals the FGF1-RUNX2 association in neural crest-derived osteoprogenitors.

Human pluripotent stem cells (hPSCs) can provide a platform to model bone organogenesis and disease. To reflect the developmental process of the human skeleton, hPSC differentiation methods should include osteogenic progenitors (OPs) arising from three distinct embryonic lineages: the paraxial mesoderm, lateral plate mesoderm, and neural crest. Although OP differentiation protocols have been developed, the lineage from which they are derived, as well as characterization of their genetic and molecular differences, has not been well reported. Therefore, to generate lineage-specific OPs from human embryonic stem cells and human induced pluripotent stem cells, we employed stepwise differentiation of paraxial mesoderm-like cells, lateral plate mesoderm-like cells, and neural crest-like cells toward their respective OP subpopulation. Successful differentiation, confirmed through gene expression and in vivo assays, permitted the identification of transcriptomic signatures of all three cell populations. We also report, for the first time, high FGF1 levels in neural crest-derived OPs-a notable finding given the critical role of fibroblast growth factors (FGFs) in osteogenesis and mineral homeostasis. Our results indicate that FGF1 influences RUNX2 levels, with concomitant changes in ERK1/2 signaling. Overall, our study further validates hPSCs' power to model bone development and disease and reveals new, potentially important pathways influencing these processes. ©2020 The Authors. Stem Cells published by Wiley Periodicals, Inc. on behalf of AlphaMed Press 2020.

Given that the human skeleton arises from different embryonic origins, modeling early bone development with human pluripotent stem cells benefits from a method of lineage‐specific derivation. This study proposes a stepwise differentiation protocol toward paraxial mesoderm‐, lateral plate mesoderm‐, and neural crest‐derived osteogenic progenitors with characterization at each stage. This approach establishes the utility of pluripotent stem cells in recapitulating osteogenesis and potential application in disease modeling. Our study's identification of transcriptomic signatures of each subpopulation reveals, for the first time, high FGF1 levels in neural crest‐derived osteoprogenitors and its influence on RUNX2, a finding that suggests its potential role in craniofacial diseases.

INTRODUCTION

Human pluripotent stem cells (hPSCs), which include human embryonic stem cells (hESCs) and human induced pluripotent stem cells (hiPSCs), are powerful tools to study developmental biology and mechanisms underlying pathological processes. A pertinent application of hPSCs is in recapitulating the organogenesis of human bone. The human skeleton originates from three distinct embryonic lineages: the frontal skull and facial bones arise from the neural crest and the remaining skeleton from the paraxial and lateral plate mesoderm. These differences in the developmental origins of bone are captured by diseases that affect specific bone types. Examples include Robinow syndrome—which affects the vertebrae of the axial skeleton—and acheiropodia—which causes truncation of the upper and lower extremities of the appendicular skeleton. There are also diseases isolated to craniofacial bones, such as Muenke syndrome. In an in vivo study utilizing a Muenke syndrome mouse model, authors reported significant shortening of presphenoid and basisphenoid bones, whereas basioccipital bones were unaffected. Through fate mapping studies, McBratney‐Owen et al demonstrated that these affected bones are neural‐crest derived, while the basioccipital bone is paraxial mesoderm in origin. Further investigation of disorders affecting particular bone types would be greatly informative of bone development and disease mechanisms, making the derivation of lineage‐specific osteogenic progenitors (OPs) from hPSCs highly beneficial. This study aimed to generate three subpopulations of OPs derived from paraxial mesoderm‐like (PM) cells, lateral plate mesoderm‐like (LP) cells, and neural crest‐like (NC) cells.

Differentiation of hPSCs directly into OPs has been reported by several studies that employ various differentiation conditions and selection criteria for purifying cell populations. , , , However, the embryonic germ layer from which they originated is often not described. Instead, two‐step differentiation protocols broadly refer to osteogenic precursors as “mesenchymal stem cell‐like” populations. This is a major limitation in recapitulating early developmental stages of bone and reducing heterogeneity in cell populations. Previous work that attempts to address this shortcoming perform stepwise differentiation of hPSCs into osteoprogenitors in serum‐ and feeder‐free conditions , , but still do not generate all three OP lineages, perform in vivo transplantation to confirm true osteogenic capacity, and/or describe transcriptomic patterns. As a result, the characterization of the differences among hPSC‐derived OPs remains incomplete. To address these gaps, our study presents a method of stepwise differentiation of lineage‐specific OPs from hiPSC‐ and hESC‐derived PM, LP, and NC cells using chemically defined and serum‐free culture conditions. We performed cell sorting coupled with gene expression analysis to optimize induction purity and employed a novel hESC‐RUNX2‐YFP reporter cell line that allowed the identification of the earliest OPs. , With our differentiation system, we characterized differences among transcriptomic patterns and highlight key markers of the three cell populations.

We also report the presence of high levels of fibroblast growth factor 1 (FGF1), an important signaling molecule in bone‐related processes, in neural crest‐derived OPs. The 23 members of the fibroblast growth factor (FGF) family bind to fibroblast growth factor receptors (FGFRs), leading to receptor dimerization and trans‐autophosphorylation of the kinase domain. Together, FGF‐FGFRs play essential developmental and homeostatic roles in the skeleton by regulating chondrocyte and osteoblast differentiation and proliferation. Indeed, aberrancy in their signaling cascades causes various well‐established skeletal diseases, such as achondroplasia from FGFR3 gain‐of‐function mutations. FGFR signal transduction is comprised of four major pathways: phosphoinositide‐3‐kinase/AKT (PI3K/AKT), phospholipase Cγ (PLCγ), signal transducer and activator of transcription (STAT), and the RAS/mitogen‐activated protein kinase (MAPK) pathways. MAPK is the predominant downstream pathway of activated FGFRs, modulating cell proliferation and, in certain contexts, differentiation. Therefore, to explore the role of endogenous FGF1 in NC‐OPs and its potential influence on bone, we investigated its effects on Runt‐related transcription factor 2 (RUNX2), a master transcription factor for osteoblast differentiation. We hypothesized that FGF1 regulates RUNX2 at multiple levels, with evidence implicating MAPK involvement, specifically extracellular signal‐regulated kinases 1 and 2 (ERK1/2 or MAPK1/3). Taken together, our study validates hPSCs as a powerful tool to model bone development and draws attention to FGF1 as a protein of interest for future studies on disorders of neural crest‐derived structures.

MATERIALS AND METHODS

Cell lines

NCRM‐5 hiPSCs were derived from male CD34+ cord blood and reprogrammed by the NIH Center for Regenerative Medicine (https://commonfund.nih.gov/stemcells/lines#RMP-generated%20iPSC%20lines). hESC line H9 (WiCell, Madison, Wisconsin) transfected with a RUNX2‐YFP reporter (hESC‐RUNX2‐YFP) was produced by the Kaufman laboratory at the University of Minnesota and was maintained as undifferentiated cells as previously described. hESC‐RUNX2‐YFP cells used to report differentiation into osteoprogenitors (OPs) were previously described.

hPSC culture and differentiation

Before stepwise differentiation into OPs, hiPSCs and hESC‐RUNX2‐YFP were replated into human xeno‐free vitronectin XF (STEMCELL Technologies, Vancouver, BC, Canada) precoated wells at 10 μg/mL in Essential 8 Media (E8) (ThermoFisher Scientific, Waltham, Massachusetts). For primitive streak‐like (PS) cell differentiation, E8 media was replaced with basal differentiation media (STEMdiff APEL; STEMCELL Technologies) supplemented with 5 μM GSKi (CHIR99021) (Stemgent, Lexington, Massachusetts) for 24 hours as previously reported. For further differentiation into PM cells, PS cells were kept in basal differentiation media supplemented with 10 μM TGF‐β inhibitor (SB431542) (Sigma‐Aldrich, St. Louis, Missouri) and BMP inhibitor (LDN193189) (Axon MEDCHEM, Reston, Virginia) for 6 days. For LP cell differentiation, PS cells were kept in basal differentiation media supplemented with 25 ng/mL recombinant human bone morphogenetic protein 4 (rhBMP4) (Peprotech, Rocky Hill, New Jersey) and recombinant human vascular endothelial cell growth factor (rhVEGF) (ThermoFisher Scientific) for 6 days. NC differentiation was accomplished in basal differentiation medium with 10 μM SB431542 and 1 μM GSKi for 6 days as previously reported.

For osteogenic commitment, PM, LP, and NC cells were kept in the osteogenic basal medium (1% P/S, 1% MEM‐NEAA, 2 mM L‐Glutamine in α‐MEM [ThermoFisher Scientific], 10% Knockout Serum Replacer [KOSR], 50 μg/mL ascorbic acid, 10 mM β‐glycerophosphate and 100 nM dexamethasone) supplemented with osteogenic mediators: 100 ng/mL BMP2, 40 ng/mL FGF9, 4 nM rapamycin (all from ThermoFisher Scientific) and 0.5 μg/mL Wnt3a (Creative Biomart, Shirley, New York) for 6 days. Subsequently, differentiation was continued with osteogenic medium without osteogenic mediators until day 28 from the start of differentiation. The medium was changed every 3 days. Cells were cultured at 37°C in 5% CO2 at 95% humidity. Bone marrow stromal cells (BMSCs) grown in osteogenic medium supplemented with 20% FBS were taken as positive controls. Undifferentiated hPSCs were taken as the negative controls.

Real‐time reverse‐transcription polymerase chain reaction

Real‐time reverse‐transcription polymerase chain reaction (RT‐PCR) analysis was done as previously described. Briefly, total RNA was extracted using Qiagen RNeasy Mini Kit (Qiagen, Valencia, California) and 1 μg of RNA was reverse transcribed into cDNA using SuperScript II Reverse Transcriptase (ThermoFisher Scientific), based on the manufacturer's instructions. Quantitative real‐time RT‐PCR was performed using 150 ng cDNA product with SYBR Green PCR Master Mix (Qiagen) in 25 μL per PCR reaction according to the recommended conditions as previously described. The genes amplified are listed in Table S1. The level of the target genes was correlated with the standard concentrations and normalized by GAPDH levels as an endogenous reference.

Flow cytometric analysis and cell sorting

A single‐cell suspension of undifferentiated and differentiated cells was prepared as previously described and evaluated for RUNX2 and surface proteins using the fluorescence‐activated cell‐sorting facility (FACSCalibur, BD, San Jose, California) in the NIDCR Combined Technical Research Core. Flow cytometry data were analyzed with the FlowJo software (Tree Star, Ashland, Oregon). Antibodies are listed in Table S2.

Immunofluorescence staining, in situ hybridization, Immunohistochemistry, and staining

Immunofluorescence staining was performed as described previously. Briefly, cells were fixed with 2% formaldehyde, washed twice with PBS, and incubated with primary antibodies for 1 hour at room temperature. Then, cells were washed three times and incubated with secondary antibodies diluted 1:100 for 1 hour at room temperature and visualized by confocal microscopy. In situ hybridization was performed according to the manufacturer's recommendation (#A001K.9905, Rembrandt Universal Dish & AP Detection Kit). Briefly, the detection of human cells in the bone formed in vivo was assessed by in situ hybridizations for human‐specific ALU repetitive DNA sequences. Immunohistochemistry was performed as previously reported. Briefly, the sections were deparaffinized and antigens retrieved with Uni‐Trieve (Innovex Biosciences, Richmond, California). Sections were blocked for 20 minutes in blocking buffer (1% BSA, 2% donkey serum, 0.1% Triton X100 in PBS). Incubations with primary antibody were done overnight at 4°C in blocking buffer. Secondary antibodies were incubated at 1:400 dilution for 1.5 hours at room temperature. Nonimmune immunoglobulins of the same isotype were used as negative controls. For H&E staining, sections of in vivo transplants were stained with H&E or toluidine blue and imaged with bright‐field microscopy as reported previously.

Enzyme‐linked immunosorbent assay, siRNA knockdown, and western blot analysis

FGF1 was quantified using Human FGF1 SimpleStep ELISA kit (Abcam, Cambridge, Massachusetts) according to the manufacturer's instructions. Briefly, the culture medium was collected over 5 days for all OPs and centrifuged at 2000g for 10 minutes. The total protein concentration of the supernatant was quantified using the Pierce BCA Protein Assay Kit (ThermoFisher Scientific). Samples and standards were loaded in duplicate in a 96‐well plate coated with an anti‐tag antibody, along with capture and detector antibodies. After a 1‐hour incubation at room temperature, wells were washed three times and 3,3′,5,5′‐tetramethylbenzidine substrate was added for 10 minutes. Stop Solution was added, and optical density was measured at 450 nm using a Varioskan LUX microplate reader (ThermoFisher Scientific). FGF1 knockdown was performed with FGF1 Silencer Predesigned siRNA (ThermoFisher Scientific). siRNA and the negative control were diluted in Opti‐MEM I reduced serum medium (ThermoFisher Scientific) and added to Lipofectamine RNAiMAX transfection reagent (ThermoFisher Scientific). After incubation for 5 minutes at room temperature, the siRNA‐lipid complex was added to NC‐OPs cultured in 6‐well plates at 37°C in 5% CO2 and 95% humidity for 72 hours. The efficiency of knockdown was assessed through ELISA. For immunoblot analysis, the protein was extracted from NC‐OPs with Extraction Buffer 5× PTR (Abcam), and total protein was measured with BCA assay. 50 μg of total protein was used for SDS‐PAGE and transferred to nitrocellulose membrane. Erk1/2 was detected using p44/42 MAPK (Erk1/2) rabbit mAb (Cell Signaling Technology, Danvers, Massachusetts) diluted 1:1000 and β‐actin rabbit pAb (Cell Signaling Technology) diluted 1:5000 served as housekeeping.

Subcutaneous transplants in mice

The use of deidentified human samples was exempted by the NIH Office of Human Subjects Research Protection (exemptions #393 and #13255). For transplant experiments, mice were approximately 8 weeks old, 25~30 g in weight and immunodeficient (NSG, NOD.Cg‐Prkdc Il2rg/SzJ, The Jackson Laboratory, Farmington, Connecticut). Transplants were constructed that contained approximately 2 million cells attached to 40 mg of the ceramic scaffold (Attrax [ceramic only], Nuvasive, San Diego, California). The anesthetized mouse was placed in ventral recumbency and the surgical area (dorsal surface) was prepared by alternating wipes of betadine and 70% ethanol three times. Autoclaved scalpel blades and scissors were used to make a 3‐cm longitudinal incision in the skin. The tips of the scissors were used to make a pocket for the transplant via blunt dissection. Sterile scaffolds (40 mg) seeded with donor cells were placed into each subcutaneous pocket. The incision was closed with an autoclip and surgical tissue adhesive. The incision site was dried with sterile gauze.

cDNA/library preparation, RNA sequencing, and analysis

Total RNA was reverse transcribed by Superscript IV (Invitrogen, Carlsbad, California) using template switching oligo and oligo dT primers followed by amplification of the second strand cDNA with LongAmp Taq polymerase (New England Biolabs, Ipswich, Massachusetts). Libraries were prepared using the Nextera XT kit (Illumina, San Diego, California), individually barcoded, pooled to a 2 nM final pooled concentration, and sequenced on a NextSeq500 instrument (Illumina) using either the 75 single‐end or the 75 × 75 paired‐end mode. After sequencing, the base‐called demultiplexed (fastq) read qualities were determined using FastQC (v0.11.2), aligned to the GENCODE v25 human genome (GRCh38.p7), and gene counts were generated using STAR (v2.5.2a). Postalignment qualities were generated with QoRTS (v 1.1.6). An expression matrix of raw gene counts was generated using R and filtered to remove low count genes (less than five reads in at least one sample). The filtered expression matrix was used to generate a list of differentially expressed genes between the sample groups using three statistical methods: DESeq2, EdgeR, and Limma‐voom.

Statistical analysis

Each experiment was repeated independently twice with three biological replicates within each experiment unless stated otherwise in the figure legends. Results were presented as mean ± SEM. Statistical analysis was performed using GraphPad Prism (GraphPad Software, La Jolla, California). One‐way or two‐way analysis of variance was used for multiple comparisons. P values were calculated by one‐tailed Student's t test, and significant differences were defined by P <.05.

RESULTS

Stepwise differentiation of OPs from hPSCs in vitro

To mimic gastrulation during which the primitive streak (PS) forms, hPSCs were differentiated into PS using the GSK inhibitor, CHIR99021, as demonstrated previously , (Figure 1A). hPSC‐derived PS cells showed significantly higher expression of primitive streak markers (BRACHYURY (T), MESP1, MIXL1, and FOXF1) at the mRNA level compared with undifferentiated hPSCs (Figure S1A, B). T, MIXL1, and MESP1 are early, reliable markers of gastrulation , , ; FOXF1 expression signifies the development of mesoderm during the late primitive streak stage.

FIGURE 1

Stepwise differentiation of hPSCs into paraxial mesoderm‐like (PM) cells, lateral plate mesoderm‐like (LP) cells, and neural crest‐like (NC) cells. A, Schematic diagram of stepwise differentiation and culture conditions of hiPSCs and hESCs toward primitive streak‐like (PS) cells, PM cells, LP cells, and NC cells. B, Flow cytometric analysis showing marker expression of PM cells derived from hiPSCs and sorting for KDR−/CD34−/CD271dim/PDGFRα+ cells. C, Flow cytometric analysis showing the marker expression of PM cells derived from hESCs and sorting for KDR−/CD34−/CD271dim/PDGFRα+/RUNX2−. D, Flow cytometric analysis of marker expression of hiPSC‐derived LP cells and sorting for CD34+/KDR+/PDGFRα−/CD271− cells. E, Flow cytometric analysis of marker expression of hESC‐derived LP cells and sorting for CD34+/KDR+/PDGFRα−/CD271−/RUNX2− cells. F, Flow cytometric analysis showing marker expression of hiPSC‐derived NC cells and sorting for CD34−/KDR−/PDGFRα−/CD271dim cells. (G) Flow cytometric analysis showing marker expression of hESC‐derived NC cells and sorting for CD34−/KDR−/PDGFRα−/CD271dim/RUNX2− cells. (B‐G) Postsorted qRT‐PCR analysis for PM, LP, and NC markers (mean ± SEM, n = 3 biological replicates, *P < .05 vs respective undifferentiated hPSCs). Differentiation efficiencies (%) are shown, with sorted population in red text. PM, LP, and NC analyses were done on day 7 of differentiation

Given that the primitive streak gives rise to paraxial mesoderm, lateral plate mesoderm, and definitive endoderm, PS cells were differentiated into PM using the BMP inhibitor, LDN193189, and TGF‐β inhibitor, SB431542 (Figure 1A), a protocol modified from Tan et al (2013). Gene expression analysis of PM cells revealed higher expression of paraxial mesoderm markers (TBX6, PDGFRα) compared with primitive streak (T, MIXL1), lateral plate (CD34, KDR), endoderm (FOXA2) and pluripotent (OCT4) markers (Figure S2A). PDGFRα has been used as a key marker of paraxial mesoderm, , , and in vivo studies have shown that TBX6 is expressed in nascent and maturing paraxial mesoderm. FOXA2 specifies endoderm in the posterior epiblast and OCT4 is expressed in embryonic stem and germ cells. Because CD34 and KDR are established markers for lateral plate mesoderm , , , PM cells were enriched by sorting for KDR−/CD34− phenotype and subsequently divided into CD271high/PDGFRα+ and CD271dim/PDGFRα+ subpopulations. BMSCs sourced from mesoderm‐derived bone, such as the iliac crest, reliably express CD271. , The CD271dim/PDGFRα+ population expressed higher paraxial mesoderm markers, TBX6 and PAX3 (Figure S2B), and was further characterized. Cells of this phenotype had significantly higher levels of paraxial mesoderm markers (CDX2, MSGN1, TBX6, and PAX3) compared with undifferentiated hPSCs (Figure 1B,C). CDX2 is essential for PM‐derived axial bone embryogenesis. PAX3 and MESOGENIN 1 (MSGN1) are master regulators of paraxial mesoderm. Because RUNX2 is the master transcription factor for osteogenic commitment, RUNX2‐YFP expression was measured to detect osteogenic commitment, of which there was none (Figure 1C).

hPSC‐derived PS cells were also differentiated into LP cells with BMP4 and VEGF, using a protocol developed by Tan et al (Figure 1A) and corroborated by others. LP cells had higher expression of known lateral plate markers (CD34, KDR) compared with the primitive streak, paraxial mesoderm, endoderm, and pluripotent markers (Figure S3A). To confirm an enriched LP population, cells were sorted into KDR+/CD34−/CD271−/PDGFRα−, KDR+/CD34+/CD271−/PDGFRα− and KDR−/CD34+/CD271−/PDGFRα− populations. The KDR+/CD34+/CD271−/PDGFRα− subpopulation showed the highest expression of the LP marker, EOMES, compared with the other two populations (Figure S3B). EOMES plays a crucial role in early gastrulation, and its deficiency results in loss of LP formation. Therefore, enriched LP cells (KDR+/CD34+/CD271−/PDGFRα−) were further analyzed for LP lineage markers (EOMES, HHEX, NKX2‐5, and ISL) (Figure 1D,E). Activated by HHEX, NKX2‐5 is reported to be a key modulator of LP maturation. Additionally, ISL1 is reported to be upstream of the sonic hedgehog pathway for LP differentiation. RUNX2‐YFP expression in hESC‐derived LP cells was analyzed and confirmed to be absent (Figure 1E).

NC cells were differentiated directly from hPSCs by using the protocol reported by Fukuta et al, which includes both TGF‐β and GSK inhibitors (Figure 1A). Subsequently, the enrichment of cells was accomplished by sorting for KDR−/CD34−. CD271 was also used as part of our selection criteria, as it is a known marker for neural crest cells. The KDR−/CD34−/CD271dim/PDGFRα− cells expressed higher neural crest markers (ALX4, SOX10) compared with KDR−/CD34−/CD271high/PDGFRα− cells (Figure S4A). ALX4 is upregulated in NC, and its mutation is associated with craniofacial disorders. SOX10 appears during neural crest migration and regulates both neural crest survival and differentiation. Therefore, KDR−/CD34−/CD271dim/PDGFRα− cells were sorted to achieve a further enriched population of NC cells. NC cells were also characterized for cranial neural crest markers (ALX4, GSC), neural crest specifier genes (SOX10, SOX9), and neural plate border genes (DIL, FOXC2, and HAND1) to confirm the neural crest lineage (Figure 1F,G); all genes were expressed significantly higher than undifferentiated hPSCs. GSC is required during embryogenesis and normal formation of craniofacial structures. , , SOX9 precedes markers of migratory neural crest. DIL acts in the dorsal part of the neural tube. , FOXC2 and HEART AND NEURAL CREST DERIVATIVES EXPRESSED 1 (HAND1) expression has been reported in neural crest cells, with the latter involved in RUNX2‐IHH‐regulated endochondral ossification. The hESC‐derived NC‐enriched population did not display signs of osteogenic differentiation (Figure 1G).

Derivation of PM‐OPs, LP‐OPs, and NC‐OPs in a serum‐free microenvironment

hPSC‐derived progenitors were further differentiated into three OP groups in serum‐free medium as previously reported : PM‐OPs, LP‐OPs, and NC‐OPs (Figure 2A). These OPs expressed significantly higher levels of RUNX2 and DLX5 compared with respective undifferentiated hPSCs (Figure 2B). DLX5 is a homeobox protein that drives osteoblast differentiation. Using flow cytometry, a significantly higher percentage of RUNX2‐YFP‐positive hESC‐derived OPs compared with undifferentiated hESCs confirms high osteogenic differentiation efficiency (Figure 2C). PM‐OPs, LP‐OPs, and NC‐OPs derived from both iPSCs and hESCs demonstrated positive staining for key bone matrix proteins OSTEOPONTIN and OSTEOCALCIN at Day 28 (Figure 2D,E). SAOS‐2, an osteosarcoma cell line, served as a positive control.

FIGURE 2

In vitro characterization of hPSC‐derived lineage‐specific osteogenic progenitor (OP) cells. A, Schematic diagram showing stepwise differentiation and culture conditions of PM, LP, and NC cells into their respective OP cells. B, Quantitative mRNA analysis of osteogenic markers, RUNX2 and DLX5 (mean ± SEM, n = 3 biological replicates, *P < .05 vs hiPSCs, #P < .05 vs hESCs). C, Analysis of osteogenic differentiation efficiency by flow cytometry of RUNX2‐YFP expression of hESC‐derived OPs (mean ± SEM, n = 3 biological replicates, *P < .05 vs hESCs). D and E, Immunofluorescence staining for OSTEOPONTIN and OSTEOCALCIN in OPs derived from their respective hPSC. Isotype was negative control. SAOS‐2 was used as a positive control. Analyses were done on 28 days after osteogenic differentiation. Scale bar = 100 μm

The osteogenic fate of OPs was validated by subcutaneous transplantation of OPs attached to ceramic particles into mice. OPs and BMSCs (positive control) formed bone at 16 weeks (Figure 3A), which was confirmed to be of human origin by in situ hybridization probes for human‐specific ALU DNA sequences (Figure 3B). PM‐OPs and LP‐OPs also formed chondrocytes after 8 weeks of transplants (Figure 3C), which may represent the cartilage intermediate in endochondral ossification that occurs in axial and appendicular skeleton. On the other hand, evidence of chondrogenic ossification in NC‐OPs was not found (Figure 3C). In contrast with endochondral ossification, intramembranous ossification does not include a cartilage template. Instead, the flat bones of the calvaria form from neural crest‐derived OPs that proliferate and condense into compact nodules containing osteoblasts, which deposit osteoid matrix that later calcifies. Overall, histological data indicates that OPs may form bone in vivo in a process similar to native bone development. However, further investigation is required for confirmation.

FIGURE 3

In vivo characterization of hPSC‐derived lineage‐specific osteogenic cells. A, H&E staining of subcutaneous transplants formed by hPSC‐derived PM‐OPs, LP‐OPs, and NC‐OPs in immunocompromised mice. Scaffold without cells served as the negative control and transplanted BMSCs as a positive control. Transplants were analyzed at 16 weeks. Scale bar = 100 μm. B, In situ hybridization for human‐specific ALU DNA sequences in calcified tissue formed by hPSC‐derived OPs transplanted with scaffold into immunocompromised mice. Transplants with BMSCs and human bone were positive controls. Negative control was a normal mouse bone. Transplants were harvested and analyzed at 16 weeks. Scale bar = 100 μm. C, Toluidine blue staining of transplants with hPSC‐derived OPs with the scaffold. Red arrows indicate hypertrophic chondrocytes, and blue arrows indicate bone. Analyses were performed at 8 weeks. Scale bar = 100 μm

The unique transcriptomic patterns of PM‐OPs, LP‐OPs, and NC‐OPs

To describe differences in lineage‐specific transcriptomic patterns, principal component analysis (PCA) was performed using the top 500 most variable genes (Figure 4A). Variation was appreciated among the three populations of OPs independent of the cell line from which they were derived, with tight grouping among experimental replicates. Clear differences were also observed in transcriptional profiles of OPs across cell lines. Among hiPSC‐derived OPs, the transcriptional profiles of LP‐ and NC‐OPs were most dissimilar. On the other hand, among hESC‐derived OPs, the transcriptional profile of LP‐ and PM‐OPs were most dissimilar (Figure 4A).

FIGURE 4

Transcriptomic patterns of PM‐OPs, LP‐OPs, and NC‐Ops. A, Principal component analysis plot showing OPs derived from hiPSCs and hESCs. B, Venn diagrams showing unique and shared genes expressed by PM‐OPs, LP‐OPs, and NC‐OPs derived from hiPSCs or hESCs or combined. C, Heat maps showing genes expressed exclusively by PM‐OPs, LP‐OPs, and NC‐OPs derived from both hiPSC and hESC lines (log2 transformed normalized expression). Analyses were performed 28 days after beginning osteogenic differentiation. After PCA, genes selected for further analysis were expressed at least 5.5× log‐fold higher than expression levels in undifferentiated cells

To further characterize these differences, the signature transcriptomic patterns of PM‐OPs, LP‐OPs, and NC‐OPs were analyzed. Any gene with at least a 5.5× log‐fold or higher expression in OPs relative to their respective undifferentiated hPSCs was selected for signature transcriptomic pattern analysis. A total of 613 uniquely expressed genes in hiPSC‐derived PM‐OPs (192 genes), LP‐OPs (277 genes), and NC‐OPs (144 genes) were found (Figure 4B). A total of 573 uniquely expressed genes were found in hESC‐RUNX2‐YFP‐derived PM‐OPs (102 genes), LP‐OPs (328 genes), and NC‐OPs (143 genes) (Figure 4B). A list of genes shared between hiPSC‐ and hESC‐derived OPs was compiled. Nine exclusive genes in PM‐OPs (Figure 4B, Table S3), 68 in LP‐OPs (Figure 4B, Table S4), and 14 in NC‐OPs (Figure 4B, Table S5) were found. Also identified were genes shared between LP‐OPs and PM‐OPs (2 genes) (Table S6), PM‐OPs and NC‐OPs (7 genes) (Table S7), and NC‐OPs and LP‐OPs (6 genes) (Figure 4B, Table S8). Lastly, 104 genes were shared by all populations of OPs derived from both hiPSCs and hESCs (Figure S4B, Table S9). These highly expressed genes shared by OPs are displayed in heat maps (Figure 4C).

Next, we selected the four highest genes expressed in both hiPSC‐ and hESC‐derived PM‐OPs, LP‐OPs, and NC‐OPs compared with undifferentiated cells. Out of nine total common genes in PM‐OPs, cadherin 6 (CDH6), chemokine C‐X‐C motif chemokine ligand 6 (CXCL6), fibrinogen beta chain (FGB), and insulin‐like growth factor binding protein 1 (IGFBP1) were most expressed (Figure 5A). The highest expressed genes in LP‐OPs were fibroblast growth factor 7 (FGF7) and a cluster of homeobox (HOX) family genes (HOXA13, HOXC10, and HOXC11) (Figure 5B). Fibroblast growth factor 1 (FGF1), keratinocyte 17 (KRT17), nerve growth factor (NGF) and pregnancy‐associated plasma protein‐A‐antisense 2 (PAPPA‐AS2) in NC‐OPs were expressed the most (Figure 5C).

FIGURE 5

Genes most expressed in PM‐OPs, LP‐OPs, and NC‐OPs. A, Box plot of the top 4 genes expressed by PM‐OPs. B, Box plot of the top 4 genes expressed by LP‐OPs. C, Box plot of the top 4 genes expressed by NC‐OPs. Analyses were performed 28 days after beginning osteogenic differentiation. Red arrowhead: >5.5× log‐fold higher vs hPSCs. Dotted red line: less than five reads in at least one sample (mean ± SEM, n = 3 biological replicates)

Significantly higher FGF1 in NC‐OPs

The link between craniofacial disorders and mutations in neural crest cells, coupled with the high expression of FGF1 in NC‐OPs, motivated deeper investigation into the role of FGF1 in NC‐OPs. Expression of FGF1 in NC‐OPs was confirmed using immunofluorescent imaging, with a high FGF1 signal in NC‐OPs compared with OPs of other lineages (Figure 6A). The data also showed more noticeable levels of FGF1 in neural crest‐derived BMSCs from the jawbone compared with those derived from the iliac crest and femur, which are both mesodermal in origin. Because FGF1 is readily exported from cells by direct translocation across the cell membrane to serve paracrine and autocrine functions, ELISA was performed on OP cell culture medium to quantify FGF1 release. In support of the fluorescent signal patterns, NC‐OPs released significantly more FGF1 extracellularly compared to PM‐ and LP‐OPs, which showed a more moderate release (Figure 6B). Craniofacial BMSCs had similarly high release and SAOS‐2 served as the positive control. FGF1 expression was then shown to be maintained in vivo by transplanted NC‐OPs (Figure 7A). FGF1 displayed affinity toward FGFR3, one potential receptor of FGF1, with a co‐localizing signal (Figure 7A). Human brain tissue was taken as a positive control for both ligand and receptor. , The reaction of the human‐specific antibody with the transplants compared with mouse calvaria served as a negative control.

FIGURE 6

High FGF1 expression by NC‐OPs in vitro and in vivo. A, Immunofluorescence staining for FGF1 expression in PM‐OPs, LP‐OPs, and NC‐OPs. Isotype control was taken as a negative control. BMSCs derived from iliac crest, femur, and a craniofacial site were osteogenically differentiated and taken as positive controls. 4′,6‐diamidino‐2‐phenylindole (DAPI) was used as a counterstain. Scale = 100 μm. Analyses were don e 28 days after the start of osteogenic differentiation. B, Enzyme‐linked immunosorbent assay for FGF1 in culture medium for PM‐, LP‐, and NC‐OPs, with SAOS‐2 as the positive control (mean ± SD, n = 2 biological replicates, *P < .05 vs PM‐OPs and LP‐OPs)

FIGURE 7

FGF1 inhibition reduces RUNX2 expression and ERK1/2 signaling. A, Immunofluorescence staining for FGF1 and FGFR3 expression in transplants generated by hPSC‐derived PM‐OPs, LP‐OPs, and NC‐OPs. Mouse calvaria and antibody isotype controls were taken as negative controls. Human brain tissue was taken as positive control. DAPI was used as a counterstain. Scale = 100 μm. Analyses were done on transplants harvested after 16 weeks. B, Measured using an enzyme‐linked immunosorbent assay, FGF1 knockdown efficiency for hPSC‐derived NC‐OPs compared to Con‐siRNA (negative control) (mean ± SD, n = 2 technical replicates). C, Flow cytometric analysis for RUNX2 in NC‐OPs treated with Con‐siRNA and FGF1‐siRNA. D, Western blot analysis of ERK1, ERK2, and β‐actin (housekeeping) after NC‐OP treatment with Con‐siRNA and FGF1‐siRNA

Inhibition of FGF1 reduces expression

Studies have demonstrated various FGFs directly stimulate RUNX2 expression and increase its binding to promoters, which in turn upregulate OP proliferation and differentiation into osteoblasts , , . Therefore, to explore the potential role of FGF1 in bone morphogenesis, NC‐OPs were transfected with FGF1‐siRNA, and its downstream effect on RUNX2 was evaluated. Knockdown efficiency, assessed through ELISA of cell extracts, showed a mean decrease in FGF1 of 57% and 49% in hiPSC‐derived and hESC‐derived NC‐OPs, respectively, compared with corresponding scrambled siRNA control groups as negative controls (Figure 7B). Subsequent flow cytometric analysis showed marked decreases in RUNX2 levels in NC‐OPs when compared with controls (Figure 7C). Previously shown to be major molecules in the FGF‐MAPK signaling pathway, MAPK3 and MAPK1 (ERK1/2) promote RUNX2 transcriptional activity through phosphorylation and subsequent activation. , Interestingly, western blot analysis showed a decrease in ERK1/2 levels in NC‐OPs after inhibition of FGF1 (Figure 7D).

DISCUSSION

This study recapitulated the development of lineage‐specific osteogenic subpopulations by stepwise differentiation of hPSCs and identified their signature transcriptomic patterns. From PS to PM, LP, and NC, the success of the differentiation system was supported by the upregulated expression of corresponding lineage markers compared with the parental lines. Although we drew upon established differentiation protocols in this study, optimization was accomplished with the addition of cell sorting coupled with gene expression analysis using reliable markers and treatment modifications that enriched populations of interest. Of note, in PM differentiation, the initial approach was the treatment of PS cells with only SB431542, a TGF‐β‐mediated SMAD 2/4 inhibitor. TGF‐β inhibition blocks endoderm differentiation and induces PM differentiation. However, high contamination with LP cells (CD34+/KDR+) was encountered. Therefore, in addition to TGF‐β, PS cells were treated with the BMP‐mediated SMAD 1/5 inhibitor, LDN193189. BMP is known to induce LP formation as opposed to PM.

The combined treatment successfully reduced LP cell contamination for more efficient derivation of an enriched population of PM cells. Interestingly, in contrast to the protocol set forth by Fukuta et al upon which ours is based, we encountered a higher expression of NC markers in CD271dim subsets as opposed to CD271high. This difference may be attributed to alternative choices in basal differentiation medium and sorting procedures. Nonetheless, higher ALX4 and SOX10 expression in our CD271dim subset compared with that of CD271high confirmed the purity of the NC population. In the final stages of differentiation, OPs clearly demonstrated osteogenic characteristics in vitro and in vivo, confirming their identity as bona fide OPs.

Lineage contamination is quite common in the differentiation of PS cells, preventing the study of lineage‐specific transcriptomic patterns. Our lineage‐specific differentiation allowed the characterization of separate OP subpopulations. To ensure that transcriptomic signatures are representative of OPs independent of the cell type from which they were derived, only highly expressed genes from both hiPSCs and hESCs were analyzed. Concerning the highest expressed PM‐OP genes, it has been shown that CDH6 is a target of TGF‐β and is regulated by RUNX2. FIBRINOGEN induces RUNX2 activity through the SMAD1/5/8 signaling pathway. CXCL6 has been reported to play an important role in bone formation during embryogenesis and in response to hormonal and mechanical stimuli. Similarly, IGFBP1 and its ligands, insulin‐like growth factors, play key roles in bone metabolism. In LP‐OPs, FGF7 is expressed in connective tissues and plays an essential role in regulating long bone development. The cluster of HOX genes is known to provide cells with specific positional identities on the anterior‐posterior axis. They also have important roles in long bone embryogenesis. Finally, KRT17 is expressed by NC‐OPs in bone marrow stromal cells. Osteoblast‐secreted NGF stimulates nerve sprouting, which points to an interesting crosstalk between the skeletal and nervous systems. , ANTISENSE A2 OF PAPPA (PAPPA‐AS2) controls the upstream signaling pathway(s) during adipogenesis. Antisense transcripts cause cis‐repression of transcription of their sense counterparts, which usually results in a negative correlation between the two transcriptional states. This may play a role in the potential of NCs to give rise to tissues of ectodermal and mesodermal (connective) lineages. When contextualized based on a particular OP's germ layer of origin, transcriptomic patterns meet expectations and further support the notion that OP differentiation remained faithful to lineage type.

Particularly interesting genes encountered from the transcriptomic analyses were those belonging to the FGF family. FGF pathways play critical roles in osteogenesis, mineral homeostasis, and ossification. , In fact, several congenital bone diseases have been directly linked to mutations in FGFs and their receptors. Demonstration of high expression of FGF7 appears unique to LP‐OPs, as other osteogenic progenitor types showed more modest expression. This phenotype is consistent with previous studies that have shown that FGF7 appears in the mesodermal embryonic mesenchyme and perichondrium, which helps give rise to the appendicular skeleton. , Building upon the understanding of the tissue‐specific expression of FGFs, we also presented, for the first time, that FGF1 is found in high levels in NC‐OPs. FGF1 has been reported in diverse cell types, including preadipocytes and astrocytes. However, to the best of our knowledge, there has not yet been evidence of endogenous FGF1 activity in neural crest OPs. As such, this finding may have important implications regarding craniofacial bone disorders. FGF2, which possesses remarkable sequence homology to FGF1, acts on neural crest cells in the development of murine frontal bone, and dysfunction of the FGF2 pathway can cause dysmorphic facial features. Mutation in the TGF‐β type II receptor in mice, for example, causes downregulation of FGF2 and compromised osteoblast differentiation of the orbital and calvarial components of the frontal bone primordium. Notably, PM‐OPs and LP‐OPs also demonstrated high expression of FGF1. However, interestingly, immunofluorescent FGF1 signal and extracellular protein levels were markedly higher in NC‐OPs compared to other lineages. This discrepancy between RNA sequencing and immunofluorescent data may be due to post‐transcriptional regulation, activity, and/or induction of FGF1 in different tissue types.

Our study provided evidence that signaling induced by FGF1 likely plays vital bone‐related functions by modulating RUNX2. The significance of RUNX2 is highlighted by mice harboring a homozygous mutation in RUNX2, which died immediately after birth and displayed a complete absence of bone formation and expression of bone markers. RUNX2 is also required for the maturation of prehypertrophic to hypertrophic chondrocytes during endochondral ossification. However, evidence that connects FGF1 with RUNX2 and describes the molecular mechanisms of regulation is limited. RUNX2 is mediated through multiple mechanisms associated with the various cascades initiated by FGFR. One study showed that treatment of osteoblast‐like MC3T3‐E1 cells with FGF2 and FGF4 strongly stimulated RUNX2 expression, with many lines of evidence pointing toward the PLCγ‐protein kinase C pathway as the responsible route. , Other studies have demonstrated that Smad‐induced junB and p38 MAPK pathway also induce RUNX2 expression after TGF‐β1 and BMP‐2 stimulation. Taken together, various cytokines and considerable crosstalk among pathways likely contribute to its finely tuned regulation. In fact, by showing that RUNX2 expression was maintained in OPs without high FGF1 levels, namely PM‐ and LP‐OPs, our study indicates that FGF1 is neither an early specifier of the OP lineage nor the sole regulator of RUNX2. In the context of osteoblastogenesis, most lines of evidence implicate FGF2 and FGF18 as other major players. Disruption of FGF2 (FGF2−/−) in mice causes marked decreases in long bone mass and rate of bone formation, with cultured FGF2−/− BMSCs displaying decreased osteoblast differentiation. , FGF18 knockout mice exhibited normal RUNX2 expression in the perichondrium/periosteum of the humerus, indicating normal OP cell numbers; however, delayed ossification was observed in the cortical bone, which arises from the perichondrium, suggesting FGF18's role in osteoblast maturation and less so in determining OP lineage. By demonstrating the reduction in RUNX2 by inhibiting FGF1, our study adds to the body of evidence that FGF1, like FGF2, 4 and 18, is a positive regulator of RUNX2 but a unique molecule in neural crest‐derived osteoblast differentiation.

Lastly, regulation of RUNX2 also occurs by phosphorylation of specific serine residues. Ser‐to‐Ala mutations at sites Ser‐301 and ‐319 can diminish its transcriptional activity at promoter regions of other osteogenic genes and interaction with other transcription factors, such as OSTERIX. Phosphorylation of RUNX2 appears to mainly involve the MAPK pathway. Transfection of preosteoblast cells with constitutively active MEK1 (MAPK kinase or MAPKK), an upstream activator of ERK1/2, corresponded with an increase in phosphorylation and activity of RUNX2, which another study demonstrated occurred by enhancing DNA‐binding capacity in the context of FGF2. The decrease in ERK1/2 after FGF1 inhibition suggested that in addition to influencing expression levels of RUNX2, FGF1 may also regulate its transcriptional activity through MAPK signaling. Although there are still many questions remaining about how FGF1 influences RUNX2, we put forward evidence that FGF1 is an important molecule in the context of bone formation in NC‐OPs. To better describe FGF1's role in the context of craniofacial bone morphogenesis, it will be necessary to further investigate FGF1's receptor binding profile, signaling cascades, and FGF1‐RUNX2 pathways in vivo.

The strengths of our study include a comprehensive, reliable list of markers against which the success of the differentiation system was measured; derivation of all three OP lineages in multiple cell lines; in vivo confirmation of OP bone‐forming potential; selection of a high threshold for highly expressed genes; and demonstration of FGF1 expression in vivo and quantification at the protein level. Reflecting lineage‐specific developmental processes, our stepwise differentiation system offers a platform to model human bone organogenesis and disease through hPSCs harboring mutations of interest. Additionally, the discovery of FGF1 in NC‐OPs presents an opportunity to explore its role in neural crest‐derived structures, including its receptor binding profile and signaling pathways. Limitations include differences in expression levels and differentiation efficiencies between hiPSCs and hESCs. However, it is not uncommon for cell types to contribute to variation in performance, even when cultured under identical conditions. For example, in a study generating osteoclasts from hESC‐ and hiPSC‐derived hematopoietic progenitors, the percentage of CD45+ cells varied widely across two hESC and one hiPSC lines tested. These differences may be due, in part, to epigenetic memory or accumulated aberrations during the reprogramming process. Importantly, however, despite these differences, we emphasize the identity and purity of final differentiated cell populations. Furthermore, although our study provides strong evidence of FGF1's influence on RUNX2, direct confirmation of its role in bone formation will be needed through in vivo studies.

CONCLUSION

An optimized method of differentiating hPSCs into lineage‐specific osteogenic populations was developed and validated through extensive characterization at every stage. Differential transcriptomic signatures and shared genes were also identified among OPs, which demonstrated osteogenic potential in vivo, the gold standard by which to determine osteogenic differentiation. The new discovery of high FGF1 found in NC‐OPs and its influence on RUNX2 indicates its potentially important role in the development of craniofacial structures. Overall, this study strengthens the utility of hPSCs to model early bone developmental processes.

CONFLICT OF INTEREST

D.S.K. declared intellectual property rights with UC‐San Diego; advisory role, research finding and ownership interest with Fate Therapeutics. The other authors declared no potential conflicts of interest.

AUTHOR CONTRIBUTIONS

F.K.: conception and design, data analysis and interpretation, manuscript writing, collection and/or assembly of data, final approval of the manuscript; B.M.: collection and/or assembly of data, data analysis and interpretation, manuscript writing; D.A.: collection and/or assembly of data, data analysis and interpretation; K.I.: assembly of data, data analysis and interpretation, manuscript writing; M. Z., L.F.D.C.D., N.C., D.M., V.D.M., M.A., K.F., R.M.: collection and/or assembly of data; S.A.: manuscript proofread; D.S.K., J.L.: data analysis and interpretation; P.G.R.: conception and design, data analysis and interpretation, manuscript writing, final approval of the manuscript.

Figure S1. Gene expression analysis of hPSC‐derived primitive streak‐like (PS) cells. (A) Expression levels of hiPSC‐derived PS markers (*P < .05 vs hiPSCs). (B) Expression levels of hESC‐derived PS markers (*P < .05 vs hESCs). All analyses were performed on Day 2 of differentiation and represented as mean ± SEM with n = 3 biological replicates.

Click here for additional data file.

Figure S2. Gene expression analysis of presorted and sorted hiPSC‐derived paraxial mesoderm‐like (PM) cells. (A) Presorted qRT‐PCR analysis for hiPSC PM markers (*P < .05 vs hiPSCs). (B) qRT‐PCR analysis as initial screen for hiPSC PM markers in sorted cells (*P < .05 vs hiPSCs TBX6; #P < .05 vs hiPSCs PAX3). All analyses were performed on Day 7 of differentiation and represented as mean ± SEM with n = 3 biological replicates.

Click here for additional data file.

Figure S3. Gene expression analysis of presorted and sorted hiPSC‐derived lateral plate mesoderm‐like (LP) cells. (A) Presorted qRT‐PCR analysis for hiPSC LP markers (*P < .05 vs hiPSCs). (B) qRT‐PCR analysis as the initial screen for the hiPSC LP marker in sorted cells (*P < .05 vs hiPSCs). All analyses were performed on Day 7 of differentiation and represented as mean ± SEM with n = 3 biological replicates.

Click here for additional data file.

Figure S4. Gene expression analysis of presorted and sorted hiPSC‐derived neural crest‐like (NC) cells. (A) qRT‐PCR analysis as an initial screen for hiPSC NC markers in sorted cells (*P < .05 vs hiPSCs ALX4; #P < .05 vs hiPSCs SOX10). All analyses were performed on Day 7 of differentiation and represented as mean ± SEM with n = 3 biological replicates.

Click here for additional data file.

Table S1. Sequences of primers used for real time RT‐PCR

Click here for additional data file.

Table S2. Antibodies used for Immunofluorescence and Flow Cytometry

Click here for additional data file.

Table S3. Showing transcriptomic patterns of hiPSC PM OPs vs hESC PM OPs

Click here for additional data file.

Table S4. Showing transcriptomic patterns of hiPSC LP OPs vs hESC LP OPs

Click here for additional data file.

Table S5. Showing transcriptomic patterns of hiPSC NC OPs vs hESC NC OPs

Click here for additional data file.

Table S6. Showing transcriptomic patterns of hiPSC PM‐LP OPs vs hESC PM‐LP OPs

Click here for additional data file.

Table S7. Showing transcriptomic patterns of hiPSC PM‐NC OPs vs hESC PM‐NC OPs

Click here for additional data file.

Table S8. Showing transcriptomic patterns of hiPSC NC‐LP OPs vs hESC NC‐LP OPs

Click here for additional data file.

Table S9. Showing transcriptomic patterns of hiPSC OPs vs hESC OPs

Click here for additional data file.