Loss of β‐Actin Leads to Accelerated Mineralization and Dysregulation of Osteoblast‐Differentiation Genes during Osteogenic Reprogramming

Abstract Actin plays fundamental roles in both the cytoplasm and the cell nucleus. In the nucleus, β‐actin regulates neuronal reprogramming by consolidating a heterochromatin landscape required for transcription of neuronal gene programs, yet it remains unknown whether it has a role in other differentiation models. To explore the potential roles of β‐actin in osteogenesis, β‐actin wild‐type (WT) and β‐actin knockout (KO) mouse embryonic fibroblasts (MEFs) are reprogrammed to osteoblast‐like cells using small molecules in vitro. It is discovered that loss of β‐actin leads to an accelerated mineralization phenotype (hypermineralization), accompanied with enhanced formation of extracellular hydroxyapatite microcrystals, which originate in the mitochondria in the form of microgranules. This phenotype is a consequence of rapid upregulation of mitochondrial genes including those involved in oxidative phosphorylation (OXPHOS) in reprogrammed KO cells. It is further found that osteogenic gene programs are differentially regulated between WT and KO cells, with clusters of genes exhibiting different temporal expression patterns. A novel function for β‐actin in osteogenic reprogramming through a mitochondria‐based mechanism that controls cell‐mediated mineralization is proposed.

technique has been widely used for visualizing single atoms of high atomic number elements supported on a lower atomic number matrix [1] . STEM imaging of ultrasmall clusters of heavy atoms affords the necessary contrast and Signal to noise ratio for visualization at electron doses that are compatible with biological specimens. In one important application, ultrasmall clusters are used for site-specific labeling of supramolecular assemblies to identify precisely the location of specific subunits within these assemblies [2] . The samples were scanned at spot size 9 and with screen current of 60 pA. The data was analysed using Velox analytical software.

C. Z-score normalization
Z-score normalization was performed in R. First, mean and standard deviation of expression values were computed for every gene. Expression values for each condition and for each gene were normalized using the z-score method to center the data: , where x is a vector of mean gene expression in all conditions (WTM, WT4, WT14, KOM, KO4, KO14) for a single gene.

D. Violin plots
Violin plots in Figure 2L-M were produced using z-score normalized data of significantly DE mitochondrial genes (n=1136) and oxidative phosphorylation genes (n=68) using R and the ggplot package. Significance of pairwise comparisons between two conditions was determined using two-tailed matched-pairs t-test.

E. Hierarchical clustering
Hierarchical clustering was performed on z-score normalized gene expression data. Agglomerate coefficients were calculated for 4 different linkage methods for each dataset to determine the one that captures the maximum structure of the data (Table S5). For each dataset, the euclidean distance matrix was calculated, Ward.D 2 linkage clustering was computed, and dendrograms were constructed. The number of clusters was determined using Silhouette plots; for most analyses where the optimal value was too low (k < 3), the next best optimal value was used in order to obtain a more granular clustering of the data. Heatmaps were then generated based on hierarchical Ward.D2 clustering and using the sub-optimal number of gene clusters. Statistical analysis was carried out using R and associated packages.

F. Gene Ontology (GO) enrichment analysis
Differentially expressed genes were selected for enrichment analysis by comparing two different conditions. Gene ontology (GO) analysis was performed using the Web-based DAVID bioinformatics resources [6] . For the GO terms to be considered over-represented or enriched in each gene list, the following criteria was applied: 1-the test p-value is less than 0.05 for the enrichment, 2-fold of enrichment (observed number of genes in the term/expected number of genes in the term) is ≥ 2 (Table S1D) or ≥ 1.5 (Table S2B-C). genes that were differentially expressed (FDR-adjusted p-value < 0.05, |fold change|>2) were subjected to GO enrichment analysis. The statistically significant GO terms (p-value < 0.05, fold of enrichment ≥ 2) are presented for biological processes (A) and molecular function (B). The GO terms shared by all 3 pairwise comparisons are highlighted. (B,D) Differentially expressed genes across the three pairwise comparisons for specific GO terms were pooled together. The list of unique genes was subjected to heatmap clustering for osteoblast differentiation (B) and calcium-ion binding (D). Clustering is based on the covariance of gene expression. Scale bar: log 2 CPM.

Figure S2: Direct reprogramming WT MEF (WTM) to osteoblast-like cells (WTO) with small molecules. (A)
Diagram of the experimental setup for direct reprogramming of WTM to WTO using three different mediums composed of OsteoMax medium and cocktails of small molecules: OsteoMAX-XF™ medium only, OsteoMAX-XF™ medium + Dex/F, OsteoMAX-XF™ medium + Dex/F/Chir. Morphology (B), calcium-based mineralization of the ECM by Alizarin Red staining (C) and the expression of osteoblast biomarkers by qPCR (D) was monitored over 14 days. Addition of Chir produced more elongated cells and reduced mineralization at day 14 ( Fig S2B-C), and the OsteoMAX-XF™ + Dex/F medium resulted in strongest mineralization at day 14 ( Fig S2C). All media induced the RNA expression of bone matrix proteins such as osteocalcin (Bglap) and osteopontin (Spp1), but none induced Runx2 expression ( Fig S2D) suggesting possible differences between the small molecules-driven reprogramming and Runx2-dependent osteoblast differentiation in vivo. Notwithstanding, the OsteoMAX-XF™ + Dex/F medium showed effective induction of a functional osteoblastic phenotype with mineralized extracellular matrix at day 14 in WT cells, which was used for all subsequent experiments.    Optimal k for mitochondriarelated genes. K=4 was selected as a sub-optimal cluster to compromise between the determined k=2 optimal clusters and enhancing the granularity of clustering analysis. (B) Optimal k for OXPHOS-related genes. K=6 was selected as a sub-optimal cluster to compromise between the determined k=2 optimal clusters and enhancing the granularity of clustering analysis. Figure S6. Exploratory RNA-seq DE analysis on osteogenesis-related genes. (A) Principal Component Analysis (PCA) shows clear separation of MEFs and osteoblast-like cells (PC1, 63% variance), as well as of WT and KO (PC2,19% variance), resulting in 6 distinct clusters for the 6 conditions. (B) Analytical pipeline for RNA-seq data normalization, differential expression and gene ontology (GO) enrichment analysis. Summary of all resulting data is provided in Table S2B-C). (C-E) Summary of GO analysis focusing on osteoblast-differentiation (OBD)-related genes (C), Wnt signaling-related genes (D) and BMP signaling-related genes (E), including the parent GO term (black), positive regulation GO terms (red), negative regulation GO terms (blue) and other relevant GO terms (grey). Numbers in parentheses reflect number of genes, and significance of GO term: *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001, (ns) P>0.05.

Figure S7. OBD, Wnt and Bmp signaling pathways in reprogrammed WT and KO cells. (A)
Analytical pipeline for RNA-seq normalization, DESeq2 DE, Z-score normalization and hierarchical clustering. (B-D) Silhouette plots for determining optimal and sub-optimal number of clusters (k) for OBD (B), Wnt (C) and BMP-related genes (D). The sub-optimal cluster number was selected for each dataset (k=5, k=4 and k=4, for OBD, Wnt and BMP, respectively  Genes that are involved in positive regulation of OBD are marked with "1", and genes that are involved with negative regulation of OBD are marked with "-1".