Generation of healthy human iPSCs and osteogenic differentiation transcriptomic data.

(A) Flowchart of iPSC establishment, iPSC-derived MSC generation, MSC to OB differentiation (preOBs, preosteoblasts; OBs, osteoblasts), and RNA-seq data generation and analyses. (B) In vitro osteogenic differentiation of MSCs (Day 0), preOBs (Day 7), and OBs (Day 21) stained with alkaline phosphatase (ALP), alizarin red, and von Kossa. (C) Expression level of known osteogenic genes at the three osteogenic stages. Data are shown as mean ± SEM. * MSC vs preOB, # preOB vs OB, or + MSC vs OB, adjusted p value < 0.0001. (D) Principal component analysis (PCA) of RNA-seq data from our iPSC-derived MSCs and differentiated OBs as well as previously published human primary MSCs and OBs, iPSCs, and other tissues in GTEx.

Gene expression profile during osteogenic differentiation.

(A) Gene expression during osteogenic differentiation with the total number of expressed genes and percentages of expressed coding and noncoding genes (middle). Number and percentages of differentially expressed (DE) and non-differentially expressed (non-DE) coding genes (left; TFs and non-TFs (white pie chart)) and of noncoding genes (right). (B) DEGs during osteogenic differentiation with total and percentages of coding and noncoding genes (left), and TFs and non-TFs (right). (C) Heatmap showing hierarchical clustering of 60 RNA-seq datasets from 20 iPSC-derived MSC, preOB, and OB lines (columns) and significant differentially expressed genes (DEGs) (rows), fold change ≥ 1.2 and adjusted p value < 0.05. Up-regulated and down-regulated gene expression is colored in red and blue, respectively. (D) and (E) Volcano plots illustrate the distribution of down- and up-regulated genes (blue and red, respectively) with adjusted p values and fold changes when comparing gene differential expression from MSC to preOB in (D) and preOB to OB stages in (E). Cutoffs of fold change ≥ 1.2 and adjusted p value < 0.05 were applied to define DEGs. The total number of up-regulated and down-regulated genes are noted at the top (red and blue, respectively). The genes are labeled for the top five up-regulated (red, bottom right), downregulated (blue, bottom left), and statistically significant down-regulated (blue, top left) and up-regulated (red, top right).

TF regulatory network in osteogenic differentiation.

(A) Principal component analysis (PCA) of all 20 healthy cell lines at three stages of osteogenic differentiation using all differentially expressed genes (DEGs). (B) PCA using only differentially expressed TF genes. (C) TF regulatory network during osteogenic differentiation. Each node represents a TF, with known bone formation associated regulators underlined in red. Two nodes are connected by a line where ReMap data suggest regulation and our RNA-seq data suggest the association between them. Nodes labeled with the gene name represent the top 100 strongest TFs based on betweenness centrality. The size of the nodes reflects the regulation strength of the TF, with the top 5 strongest circled in pink. (D) and (E) Top significantly enriched GO BP terms in (D) and Reactome pathways in (E) of TFs in each network module.

Gene co-expression network in osteogenic differentiation.

(A) Sunburst plots represent the hierarchy structure of the MEGENA co-expression network constructed on gene expression during the osteogenic differentiation. The structure is shown as concentric rings where the center ring represents the parent modules, and outer rings represent smaller child modules. Subnetwork modules are colored according to the enrichment of differential gene expression between stages (FDR < 0.05; left, MSC to preOB; right, preOB to OB; blue, down-regulated; red, up-regulated; cyan, both up-and down-regulated). The subnetwork branch for Module M204 is outlined and labeled. (B) Co-expression network Module M204. Diamonds indicate KNR genes, and circles indicate non-KNR genes. Blue indicates DEGs from MSC to preOB stages, red indicates DEGs from preOB to OB stages, and cyan indicates shared DEGs for both comparisons. Genes known to be related to bone are in red.

Inhibitory role of Klf16 in osteogenic differentiation in vitro and in vivo.

(A) The expression of KLF16 at three human osteogenic differentiation stages. Data are shown as mean + SEM. * MSC vs. preOB, # preOB vs. OB, or + MSC vs. OB, adjusted p value < 0.001. (B) Analysis of Klf16 expression by RT-qPCR in MC3T3-E1 cells transduced vectors containing either stuffer sequence or Klf16 cDNA. Data are presented as the mean ± SEM (n=3, unpaired t-test p value < 0.01) (Supplementary Table 7). (C) Osteogenic differentiation of MC3T3-E1 cells without or with overexpression of Klf16, stained for ALP at Day 7 and alizarin red and von Kossa at Day 14 and Day 21. (D) Length, fat mass, and lean mass of wild type (WT) and Klf16+/- mice. (E) DEXA analysis of whole body (head excluded) bone mineral content (BMC), bone area (B-area), and bone mineral density (BMD) of WT and Klf16+/- mice. (F) and (G) Representative microCT images of distal femur trabecular bone in (F) (left, top view; right, side view) and cortical bone in (G) from WT and Klf16+/- mice. Scale bar: 1 mm. (H) and (I) Graphs show trabecular bone volume/tissue volume (BV/TV), trabecular thickness (Tb.Th), trabecular number (Tb.N), and trabecular separation (Tb.Sp) in (H); cortical bone area (Ct.Ar), cortical periosteal perimeter (Ct.Pe.Pm), and cortical endosteal perimeter (Ct.En.Pm) in (I). For (D), (E), (H), and (I), data are presented using BoxPlotR (Spitzer et al., 2014) as the mean ± SEM (WT n = 6, Klf16+/- n = 6, 3 males and 3 females for each group, aged 17 weeks, paired t-test, n.s, not significant, * p value < 0.05, ** p value < 0.01).

Characterization of healthy human iPSCs and iPSC-derived MSCs.

(A) Identification of stem cell markers: ALP staining, fluorescent immunocytochemistry (ICC) of pluripotency markers TRA-1-60, SOX2, OCT4, SSEA4, and NANOG. Scale bar: 100 μm. (B) Embryoid body (EB) formation in suspension from aggregates of iPSCs. ICC of three germ layer markers: AFP, endoderm; α-SMA, mesoderm; and TUBB3, ectoderm. Scale bar: 50 μm. (C) Fluorescence-activated cell sorting of CD105+/CD45- MSCs. FITC, fluorescein isothiocyanate; PE, phycoerythrin. (D) The percentage of mesenchymal surface marker CD105+/CD45- cells in total cells differentiated from iPSCs originated from PBMCs and fibroblasts. Data are shown as mean + SEM. (E) and (F) Analyses of MSCs after sorting and expansion with flow cytometry for MSC negative markers (CD31, CD34, and CD45) in (E) and MSC positive markers (CD29, CD73, CD90, and CD105) in (F) labeled with different fluorochromes (V450, APC, PE, and FITC). Marker expression is presented as histograms (green). Unstained cells were used as controls (gray). Data are shown as mean + SEM. (G) Gene expression of MSC markers in transcripts per kilobase million (TPM). (H) Spindle-like morphology of iPSC-derived MSCs. Scale bar: 50 μm.

Gene Ontology (GO) enrichment of TF regulatory network modules.

(A) Top enriched GO cellular component (CC) terms of each module. (B) Top enriched GO molecular function (MF) terms of each module.

Gene expression pattern across osteogenic differentiation stages using pseudobulk single-cell RNA-seq data in Housman et al., 2022.

KNR_TF_up: top five up-regulated KNR transcription factors; KNR_TF_down: top five down-regulated KNR transcription factors; TFs: top five transcription factors (TFs) based on TF regulatory network; osteogenic_markers: five known osteogenic markers.

Expression of osteoclastogenesis markers in the femur bone of WT and Klf16+/- mice.

(A) IHC for CALCR in the distal femur bone. The CALCR-positive areas in epiphysis and metaphysis of two Klf16+/- mice and two WT mice were quantified and shown as the percentage of CALCR-positive area (fluorescence-positive area/total bone area) in (B). (C) IHC for RANKL in the distal femur bone. The percentage of RANKL-positive area were quantified and shown in (D). For (B) and (D), three replicate sections from each animal were analyzed. Data are presented as the mean ± SEM and visualized with GraphPad Prism 10. n.s: not significant. * p value < 0.05 by unpaired t-test. Scale bar: 500 µm.