Single-Cell Atlas Reveals Fatty Acid Metabolites Regulate the Functional Heterogeneity of Mesenchymal Stem Cells

Bone marrow mesenchymal stem cells (MSCs) are widely used clinically due to their versatile roles in multipotency, immunomodulation, and hematopoietic stem cell (HSC) niche function. However, cellular heterogeneity limits MSCs in the consistency and efficacy of their clinical applications. Metabolism regulates stem cell function and fate decision; however, how metabolites regulate the functional heterogeneity of MSCs remains elusive. Here, using single-cell RNA sequencing, we discovered that fatty acid pathways are involved in the regulation of lineage commitment and functional heterogeneity of MSCs. Functional assays showed that a fatty acid metabolite, butyrate, suppressed the self-renewal, adipogenesis, and osteogenesis differentiation potential of MSCs with increased apoptosis. Conversely, butyrate supplement significantly promoted HSC niche factor expression in MSCs, which suggests that butyrate supplement may provide a therapeutic approach to enhance their HSC niche function. Overall, our work demonstrates that metabolites are essential to regulate the functional heterogeneity of MSCs.

The versatile functions of MSCs enable increasing clinical applications of MSCs to treat diseases such as graft versus host disease (GvHD), Crohn's disease, heart failure, bone marrow failure syndrome, and osteogenesis imperfecta and bone fractures (Galipeau and Sensebe, 2018;Andrzejewska et al., 2019;Martin et al., 2019;Yin et al., 2019). MSC infusion is also applied for facilitating implantation after HSC transplantation (Mendez-Ferrer et al., 2010;Zhou et al., 2017). However, the diverse therapeutic targets challenge the clinical trials of MSCs in which the cell functional heterogenicity, culture methods, and expansion levels could potentially influence the therapeutic consistency and limit clinic efficacy of MSCs (Pattappa et al., 2013;Moll et al., 2014;Yuan et al., 2019).
Recent works show that butyrate, a natural short-chain fatty acid (SCFA) produced by mammalian intestinal microbiota, can inhibit histone deacetylase (HDAC) activity and impairs intestinal epithelial stem/progenitors in wound repair in vivo (Kaiko et al., 2016) but promotes iPSC reprogramming efficiency (Liang et al., 2010). However, the role of butyrate in regulating MSCs remains elusive. Here, using single-cell RNA sequencing (scRNA-seq), we identified FA pathways that are involved in lineage commitment of MSCs, and further functional assays prove that metabolite butyrate alters MSC cell fate in self-renewal, apoptosis, and HSC niche factor expression.

Mice
C57BL/6 mice were obtained from the Jackson Laboratory and were maintained in the C57BL/6 background. All animal experiments were performed according to protocols approved by the institutional animal care and use committee.

Bone Marrow Digestion
Bone marrow digestion was performed as described with small changes (Zhou et al., 2014). In brief, intact bone marrow from mice at the age of 6-8 weeks were flushed from mouse femora and tibiae and subjected to two rounds of enzymatic digestion at 37 • C for 20 min each. The digestion buffer contained 0.2 mg/ml liberase (Roche) and 200 µg/ml DNAse I (Roche) in 1 × HBSS with calcium and magnesium. The digested marrow cells then underwent red blood cell lysis using 0.16 M ammonium chloride solution.

scRNA-seq
Sorted CD45 − Ter-119 − CD31 − PDGFα + CD51 + single cells were processed through the Chromium Single Cell Platform using the Chromium Single Cell 3 Library and Gel Bead Kit v3 (10X Genomics, PN-1000075) and the Chromium Single Cell B Chip Kit (10X Genomics, PN-1000074) as the manufacturer's protocol. In brief, 15,000 cells were loaded onto the chromium instrument to generate single-cell barcoded droplets. Cells were lysed and barcoded reverse transcription of RNA was occurred. Libraries were prepared by following the amplification, fragmentation, adaptor, and index attachment and then sequenced on an Illumina NovaSeq platform. The scRNA-seq data generated in this study are deposited in GEO (GSE167035 1 ).

scRNA-seq Data Processing
The scRNA-seq reads were aligned to the mm10 reference genomes, and unique molecular identifier (UMI) counts were obtained by CellRanger 3.0.2. Normalization, dimensionality reduction, and clustering were performed with the Seurat 3.0 R package (Butler et al., 2018) on RStudio, and 3005 of 11,888 cells with Ptprc (CD45) expression were excluded to remove potential contamination. Cells were filtered to have > 200 and < 7000 detected genes and less than 5% of total UMIs mapping to the mitochondrial genome. Data set normalization was performed by dividing the UMI counts per gene by the total UMI counts in the corresponding cell and log-transforming, followed by the results scaling and centering. Cells underwent dimensionality reduction with the t-stochastic neighboring embedding (tSNE) method and partition-based graph abstraction (PAGA) using scanpy (Wolf et al., 2018). Integration of published MSC and immune cell scRNA-seq with our data was performed by Seurat with the function SCTransform() and Harmony algorithm. Signature genes of each cluster were obtained using the Seurat function FindMarkers with the Wilcox test. All correlations were calculated based on values with the function cor() and the parameter "method = 'spearman'." The pseudotime trajectory was analyzed by monocle2 based on the Seurat clustering (Qiu et al., 2017). A signature gene heat map was generated by pheatmap R packages. GOChord plots and GOClust plots were generated by GOplot R packages. Gene set enrichment analysis (GSEA) was performed using the gsea R package (Subramanian et al., 2005). Gene lists were preranked by the fold change values of the differential expression analysis using Seurat. GSEA plots were generated by enrichplot R packages. Gene sets were obtained from the gene ontology database as indicated. Signature genes feature plots and violin plots were generated with Seurat R packages.

Cell Counting Kit-8 (CCK-8) and Lactate Dehydrogenase (LDH) Activity Assay
Freshly isolated marrow cells were plated at a density of 2 × 10 4 cells/well in a 96-well plate and cultured with or without butyrate, ZVAD, or Nec-1 as indicated for 24 h in 37 • C, 5% O 2 , and 5% CO 2 chambers. For the CCK-8 assay, 10 µl CCK-8 reagent (Solarbio) was added in wells and continued incubated for 4 h in chambers. The optical density was measured at 450 nm using a microplate reader. LDH activities were performed following the manufacturer's instructions (Njjcbio). Cell death frequency was calculated as LDH release into cell culture medium dividing LDH in the total cell lysate.

Statistical Analyses
Data are presented as mean ± s.e.m. For multiple comparisons analysis, data were analyzed by repeated-measures one-way ANOVA followed by Dunnett's test. The difference was considered statistically significant if p < 0.05, and P < 0.05, P < 0.01, P < 0.001. For other experiments, except for scRNAseq analysis, data were analyzed by two-tailed Student's t test. The difference was considered statistically significant if p < 0.05, and * p < 0.05, * * p < 0.01, * * * p < 0.001.

Single-Cell Atlas Identifies the Heterogeneity of MSCs
To explore the heterogeneity of bone marrow MSCs, we applied droplet-based scRNA-seq with non-hematopoietic (CD45 − Ter-119 − ), non-endothelial (CD31 − ) bone marrow cells that express MSC markers (PDGFα + CD51 + ) (Morikawa et al., 2009;Yue et al., 2016;Boulais et al., 2018;Lee et al., 2018) ( Figures 1A,B). We successfully detected a total of 8883 cells with an average of ∼2000 genes per cell in bone marrow MSCs (Supplementary Figure 1A). We then employed t-SNE and identified 11 clusters in MSCs ( Figure 1C). Results of PAGA ( Figure 1D) and enriched gene ontology (GO) (Figure 1E) annotated the 11 clusters into six populations, including pre-MSCs, adipogenic MSCs, chondrogenic MSCs, osteogenic cells, angiogenic MSCs, and immunomodulating MSCs ( Figure 1F). Pre-MSCs, constituted of clusters 6 and 8, enriched pathways in protein transport, nuclear transport, and ribosome biogenesis, which are critical for MSC lineage commitment (Herencia et al., 2012;Chen et al., 2015). Adipogenic MSCs, chondrogenic MSCs, and osteogenic cells, constituted of clusters 1 and 5, cluster 3, and clusters 4 and 9, enriched pathways in fat cell differentiation, chondrocyte development, and bone development, respectively. Angiogenic MSCs were constituted of cluster 2, enriching pathways in endothelial and epithelial cell migration. Immunomodulating MSCs were constituted by clusters 7, 10, and 11, which enriched pathways associated with leukocyte proliferation and myeloid cell homeostasis ( Figure 1E). Furthermore, these immunomodulating MSCs overlapped with other MSC subpopulations and published bone marrow MSCs (Baryawno et al., 2019;Leimkuhler et al., 2020) but not with immune cells (Choi et al., 2019) ( Figure 1G). In line with this, immunomodulating MSCs did not express any immune cell markers, such as T cells, B cells, NK cells, and macrophages (Supplementary Figure 1B). This ruled out the potential contamination of immune cells in this immunomodulating MSC population. Furthermore, adipogenic MSCs, chondrogenic MSCs, and osteogenic cells were referred to as lineage-committed MSCs. Consistently, our unsupervised trajectory analysis by Monocle 2 showed that the pre-MSCs clustered in the root of the trajectory, and adipogenic, chondrogenic, osteogenic, and angiogenic MSCs clustered in different branches ( Figure 1H).
Overall, using scRNA-seq, we identified pre-MSCs and lineage-committed MSC clusters and revealed their potential regulatory mechanisms in MSC lineage commitment.
Overall, these findings indicate the distinguished regulation mechanisms for multipotent pre-MSCs and lineage-committed MSCs.

The FA Metabolic Process Regulates MSC Lineage Commitment
To explore the underlying mechanism in regulating different MSC clusters, we performed GSEA on our scRNA-seq data between pre-MSCs and lineage-committed MSCs using gene sets in the GO database (Ashburner et al., 2000; The Gene Ontology Consortium, 2019). As expected, pre-MSCs had much less osteoblast cell differentiation (GO: 0001649) and fat cell differentiation (GO: 0045444) genes compared with lineage-committed MSCs (NES = -1.60 and -1.56, respectively, Figures 3A,B). Furthermore, we discovered that pre-MSCs are relatively quiescent as they had much less activated cell cycle (GO: 0045787) genes (NES = -1.45, Figure 3C). We next investigated how metabolites regulate pre-MSCs and lineage-committed MSCs in their stem cell fate decision. The GSEAs showed that pre-MSCs had less activated genes under the carbohydrate metabolic process (GO: 0005975) compared with lineage-committed MSCs (NES = -1.28) but were not significantly different under glycogen (GO: 0005977) or amino acid metabolic process (GO: 0006520) (Figures 3D,E). However, pre-MSCs significantly enriched genes under the FA metabolic process (GO: 0006631) (NES = 1.61, Figure 3F). We also confirmed that FA metabolic genes are much enriched in pre-MSCs compared with lineage-committed MSCs, such as Pam, Cyp1b1, and Lpl (Satani et al., 2003;Rozovski et al., 2015;Bushkofsky et al., 2016) (Figure 3G). These indicate that the FA metabolic process might contribute to the regulation of pre-MSCs and lineage-committed MSCs through their distinguished metabolic profile.
Overall, our scRNA-seq data indicate a novel potential role of MSCs responding to FA treatment, which could be applied to in vitro culture.

Supplement of Butyrate Suppresses Self-Renewal and Differentiation Potential of Bone Marrow MSCs
To investigate whether supplement butyrate influences the maintenance of MSCs during in vitro culture, we performed a bone marrow MSC CFU-F assay and found that supplement of butyrate at biological serum concentration (500 nM) (Tyagi et al., 2018) or higher concentration (5 µM and 500 µM) effectively reduced the CFU-F clone numbers compared with vehicle control (33.5%, 31.7%, and 60.8% decreased, respectively; Figures 4A,B). These show that butyrate suppressed MSC proliferation and self-renewal potential.   We next asked whether butyrate impairs the self-renewal potential of MSCs due to enhanced differentiation. To this aim, we performed an ex vivo differentiation assay to examine the osteogenesis and adipogenesis capacities of MSCs. Interestingly, butyrate-treated MSCs have reduced osteogenic and adipogenic differentiation ability (Figures 4C,D). Consistently, the block of osteogenic differentiation ability was confirmed by multiple osteogenic-specific marker genes, such as Runx2, Ocn, and Col1a1, which were markedly increased after osteogenic differentiation (7.75-fold, 83.7-fold, and 3.19fold increased, respectively, compared with undifferentiated MSCs) in control MSCs but were inhibited in butyrate-treated MSCs (76.1%, 96.2%, and 92.3% decreased, respectively, compared with differentiated control MSCs, Figures 4C,E). Furthermore, the adipogenic-specific marker genes, such as Fabp4, Adiponectin, and Cebpa, were increased in control MSCs (219-fold, 2071-fold, and 27.9-fold increased, respectively, compared with undifferentiated MSCs) but were dramatically compromised in butyrate-treated MSCs (70.3%, 77.2%, and 64.0% decreased, respectively, compared with differentiated control MSCs; Figures 4D,F).
Overall, our data show that supplement of butyrate during in vitro culture impaired MSCs in their self-renewal potential and differentiation abilities.

Supplement of Butyrate Triggers Apoptosis but Promotes HSC Niche Factor Expression in Bone Marrow MSCs
To investigate the potential mechanisms of butyrate treatment on MSCs, we performed 7AAD cell death staining to examine whether butyrate treatment influences their survival. We (F) Quantitative RT-PCR of transcript levels of adipogenic genes (normalized to day 0 control group). Repeated-measures one-way ANOVA followed by Dunnett's test for multiple comparisons in (B). p < 0.01, p < 0.001. Two-tailed Student's t tests were used to assess statistical significance in (E,F). *p < 0.05, **p < 0.01, ***p < 0.001.
Overall, our data show that butyrate supplement dichotomously regulates MSCs in their self-renewal and HSC niche function potentially through altering the MSC metabolic status and inhibiting HDAC activity.

DISCUSSION
MSCs are present in all organs and tissues, and MSCs are a highly heterogenous subset (Uccelli et al., 2008). The heterogeneity of MSCs is closely related to their clinical utilities and also determines the barriers in transferring MSC capacities into the clinic (Costa et al., 2020). The functional heterogeneity of MSCs has been indicated in previous studies. For example, IL-17 + MSCs have enhanced antibacterial effects but with reduced immunosuppressive function compared with bulk MSCs due to altered NFκB-TGF-β signaling (Yang et al., 2013). Studies using RNA fluorescence in situ hybridization (FISH) (Cote et al., 2016) and fluorescent probes  indicate that canonical markers are tenuously linked to the differentiated phenotypes, and it is difficult to use single markers to predict functional potential. Recent advanced single-cell studies have further explored the heterogeneity of bone marrow MSCs with distinct differentiation potential (Baryawno et al., 2019;Leimkuhler et al., 2020). Furthermore, researchers also identified an IL-10 regulated metabolically active mature adipocyte subtype from subcutaneous adipose tissue (Rajbhandari et al., 2019) and Runx2 + /Gli1 + cells in the adult mouse incisor, which maintains Gli1 + MSCs .
The metabolic profile determines the functional heterogeneity of MSCs (Costa et al., 2020). Genetic inhibition of mitochondrial complex III in human MSCs and murine adipocyte precursor cells impacts adipocyte differentiation (Tormos et al., 2011;Joffin et al., 2021). Glutamine metabolism regulates proliferation and osteoblast-adipocyte lineage determination (Yu et al., 2019). Among the diverse tissues from which MSCs could be isolated, the most common source tissue is the bone marrow (Yin et al., 2019). To meet the clinical requirement for MSCs in cell therapy, MSCs have to undergo a rapid cultural expansion and longterm cryopreservation, which are largely different compared with their biological microenvironment (Yuan et al., 2019). The metabolic profile determines MSC cell fate and heterogeneity (Pattappa et al., 2013;Moll et al., 2014;Costa et al., 2020). Interestingly, we discovered that pre-MSCs, despite the reduced expression of proliferation and differentiation genes, enriched genes in FA metabolic process compared with lineage-committed MSCs. This finding indicates that supplement of metabolites, such as FA in MSC during in vitro culture, may impact MSC functional heterogeneity.
Recent studies indicated that butyrate, one of the metabolites produced in the healthy intestinal lumen (Jacobi and Odle, 2012), enhances the effect of parathyroid hormone (PTH) to support bone formation (Fan et al., 2017;Li et al., 2020;Pacifici, 2020). Furthermore, butyrate promotes Treg cell regeneration and differentiation, which stimulates bone formation by activating Wnt signaling in osteoblasts (Arpaia et al., 2013;Furusawa et al., 2013;Tyagi et al., 2018). In our study, we discovered that supplement of butyrate suppressed the selfrenewal capacity and differentiation potential toward osteoblasts and adipocytes in MSCs.
HDACs, a series of critical transcriptional cofactors modulating gene expression by deacetylating histones and transcription factors, participate in stemness maintenance, lineage commitment determination, cell differentiation, and proliferation as well as other activities in normal hematopoiesis (Akimova et al., 2012;Cortiguera et al., 2019). Homozygous deletion of HDAC3 in Prrx1-expressing cells reduced chondrocyte and osteoblast differentiation in vitro (Feigenson et al., 2017). HDAC inhibitors exhibit antitumor activity for multiple myeloma (Song et al., 2011;Santo et al., 2012) and B cell lymphoma (Cortiguera et al., 2019) by suppressing cells survival and differentiation and inducing apoptosis. Butyrate is one of the extensively studied HDAC inhibitors (Steliou et al., 2012). Butyrate blocks the activity of class I and II HDACs and increases histone acetylation globally in multiple types of cells, including CD8 + T cells, B cells, hepatocytes, and some tumor cell lines, such as MCF-7 (human breast cancer cells) and HCT116 (human colon carcinoma cells) (Candido et al., 1978;Davie, 2003;Fellows et al., 2018;Luu et al., 2018;Ji et al., 2019;Sanchez et al., 2020). In our work, we discovered that butyrate supplement reduced MSC proliferation and differentiation abilities at biological serum butyrate concentration without cell death. Butyrate treatment reduced the expression of proliferation gene Ccnb1 but upregulated apoptotic genes, including Bax, Bak1, and Bid, which might explain the butyrate-induced growth inhibition. Furthermore, we found that lineage-committed MSCs enriched more genes associated with HDAC activity than pre-MSCs did. HDACs support cell growth for multiple tumor cells, including diffuse large B cell lymphoma, lung adenocarcinoma, and breast cancer (Gupta et al., 2012;Lapierre et al., 2016;Wang et al., 2016); therefore, HDAC inhibitors are clinically used for cancer treatment (Khan and La Thangue, 2012;Li and Seto, 2016). Our findings show that, unlike quiescent pre-MSCs, lineage-committed MSCs are more sensitive to HDAC inhibition, potentially due to their high proliferation potential.
Our data show that a high dose (500 µM) of butyrate treatment increased apoptosis in MSCs and also increased their HSC niche function. MSCs engulfed apoptotic bodies to enhance their differentiation ability and ameliorate the ovariectomyinduced osteopenia through activation of the Wnt/β-catenin pathway (Liu D. et al., 2018). Moreover, MSC apoptosis is related to enhanced osteoblast differentiation (Schaub et al., 2019), and apoptotic MSCs have enhanced immunosuppression activity when infused in patients with GvHD (Galleu et al., 2017;Burnham et al., 2020). These findings inspired us that the apoptosis of MSC might benefit their tissue regeneration or immunomodulation functions, which is consistent with our observation that a high dose of butyrate supplement increased apoptosis and enhanced their HSC niche function. Inhibition of HDAC2 and 3 promotes the proliferation of hematopoietic stem and progenitor cells (HSPCs) (Dhoke et al., 2016;Wang et al., 2020), but whether HDAC inhibition alters the HSC niche remains unclear. In our work, we demonstrate that in vitro supplement of butyrate-promotes HSC niche factor expression in MSCs, including Kitl, Angpt1, Cxcl12, and Jag2. As bone marrow MSCs produce SCF (Kitl), CXCL12 (Cxcl12) or Angiopoietin-1 (Angpt1) and endothelial cells secrete Jagged-2 (Jag2), which are critical for HSC maintenance (Greenbaum et al., 2013;Zhou et al., 2015;Guo et al., 2017;Comazzetto et al., 2019), our observation indicates that butyrate supplement may enhance HSC niche function potentially through inhibiting HDAC in MSCs. Overall, our findings indicate the possibility that the application of butyrate in MSC culture can amplify their HSC niche function and shed light on MSC treatments for patients with ineffective hematopoiesis and patients who underwent HSC transplantation.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GEO (GSE167035, https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE167035).

ETHICS STATEMENT
The animal study was reviewed and approved by the Sun Yatsen University.

AUTHOR CONTRIBUTIONS
JX performed the scRNA-seq data analysis and the functional assays, generated figures, and wrote the manuscript. QL, SX, and QX contributed to the scRNA-seq library construction and functional assays. YZ, YL, and LY contributed to the functional data analysis. JW performed and analyzed the scRNA-seq. LJ, LM, DL, and MZ supervised the project. All authors contributed to the article and approved the submitted version.