Identification of key genes of human bone marrow stromal cells adipogenesis at an early stage

Expression profile data

The raw gene profile dataset GSE107789, based on the platform of Agilent-028004 SurePrint G3 Human GE 8 × 60 K Microarray chip (GPL14550), was downloaded from GEO database (http://www.ncbi.nlm.nih.gov/geo/). The data were from human BMSCs line (hMSC-TERT) undergoing adipogenic differentiation on day 7. Total RNA of three replicates from control and adipocytic differentiation group was extracted and hybridized to the one-color Agilent Human GE 8 × 60 K Microarray chip.

Analysis of DEGs

Based on the data from the GEO dataset GSE107789, we analyzed DEGs in BMSCs undergoing adipogenic differentiation vs. BMSCs treated with control agent on day 7. Significant values for DEGs were obtained using web analysis tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) (Barrett et al., 2013). The adjusted p values were calculated using Benjamini and Hochberg false discovery rate method by default. Hence, fold changes were transformed to log2 values for volcano plot. |Logarithm of fold change (LogFC)| cutoff > 1.25 and adjusted p < 0.01 were used to determine significantly changed transcripts.

Functional and pathway enrichment analysis

Gene Ontology terms, including biological process (BP), cellular component (CC), and molecular function (MF), have been widely used for gene annotation and characterization of high-throughput genome or transcriptome data (Ashburner et al., 2000). KEGG is a website database for pathways enrichment analysis of target genes sets, to obtain integration and interpretation of these data sets (Kanehisa et al., 2012). The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) was used to identify GO terms and KEGG pathway enrichment analysis of DEGs (Da Huang, Sherman & Lempicki, 2009). Moreover, GSEA software (GSEA 4.0.3, https://www.gsea-msigdb.org/gsea/index.jsp) was used to analyze the dataset as well. KEGG subset of canonical pathway (c2.cp.kegg.v7.0) from Molecular Signatures Database was used for analysis and gene sets with normalized enrichment score >1 and normalized p < 0.01 were considered as statistical significance. Leading edge analysis was carried out to screen out the key genes.

PPI network analysis

Protein–protein interactions network analysis was performed in order to identify crucial proteins and protein modules involved in BMSC adipogenic differentiation. In the present study, STRING database (http://www.string-db.org/) was used to converse DEGs to corresponding protein and PPI network (Szklarczyk et al., 2015). In consideration of the 2,000 genes count limit of STRING, DEGs set containing 1,668 genes with |LogFC| > 2 was put into analyzing. Cytoscape (https://cytoscape.org/), an open source software platform for visualizing molecular interaction networks and biological pathways and integrating these networks with annotations, gene expression profiles and other state data, was used to visualize and analyze topologic parameters, including degree of distribution, closeness centrality and betweenness centrality of PPI data obtained from STRING database. Nodes (standing for corresponding proteins) with a higher degree which form more edges (standing for PPIs) with other proteins were defined as hub proteins. Moreover, module analysis of the present PPI network was performed, using MCODE, a plugin of Cytoscape, with the following settings: degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max. depth = 100 (Bader & Hogue, 2003). Interested modules were screened out, manually by MCODE grade and proteins involved in the corresponding module, and further investigated by GO terms and KEGG pathway enrichment analysis.

Experimental animals

C57BL/6J mice were provided by the Experimental Animal Center at Southern Medical University, Guangzhou, China. All animal protocols were approved by the Animal Care and Use Committee at Nanfang Hospital, Southern Medical University (NFYY-2018-34) and all animal experiments were performed in accordance with the guidelines of Animal Care and Use Committee. Mice were raised under specific pathogen-free conditions at 18–28 °C with humidity of 40–70%, free access to food and water with a 12 h light/dark cycle. 12-week-old mice were randomly divided into two groups with bilateral ovariectomy (OVX) and sham surgery. At week 8 after the surgery, mice were euthanized and bone samples were harvested for further analysis. Bone marrow of tibias and femurs from the left side was flushed out and stored at −80 °C for mRNA expression analysis. Tibias and femurs from the right side were dissected free of soft tissue and processed for histological analysis.

Isolation and culture of BMSCs

Bone marrow stromal cells were isolated from 4–6 week-old mice following the protocol previously described (Maridas et al., 2018). Briefly, after both ends of tibia and femur were removed, bone marrow was flushed out with α-MEM under sterile condition. Flushed cells were passed through 70 μm cell strainer and seeded at a density of 1 × 10⁶ cells/cm² in 100 mm dish for expanding. Cells were cultured with growth medium (α-MEM with 15% fetal bovine serum, 55 μM β-mercaptoethanol, 2 mM glutamine, 100 IU/ml penicillin and 100 μg/ml streptomycin). For adipogenic differentiation, cells were seeded at a density of 2 ×10⁵ cells/ cm² in six-well plates with growth medium. When cells reach 100% confluence, medium was replaced with adipogenic differentiation media (MUBMX-03031, Cyagen Biosciences Inc., China) following manufacturer’s instructions. Cells were harvested for mRNA expression analysis or Oil red staining at day 7 after adipogenic differentiation.

Oil red O staining

Lipid droplets were detected by Oil Red O (ORO) staining. BMSCs were washed with PBS and fixed with 4% paraformaldehyde for 30 min at room temperature, washed with distilled water and 60% isopropanol, and then incubated with ORO staining solution (0.3% ORO powder in 60% isopropanol, Sinopharm Chemical Reagent Co. Ltd., China) for 30 min at room temperature. For ORO staining in bone tissue, tibias were fixed, decalcified followed by frozen embedding. Sections at 8-μm were used for ORO staining. Images were acquired under a BX53 microscope (OLYMPUS Co., Japan). Quantification of ORO was obtained by ImageJ software. For cell culture experiments, five different fields were randomly imaged and counted for ORO-positive cells and the total cells. For frozen sections, the ratio of ORO-positive staining area to the total area of each photograph was calculated.

RNA isolation and real-time quantitative PCR analysis

Total RNA was isolated using RNAiso Plus (Takara, Kyoto, Japan) and reverse-transcribed using the PrimeScript RT reagent Kit (Takara, Kyoto, Japan) according to the manufacturer’s instructions. The cDNA samples were then used for real-time quantitative PCR (qRT-PCR) analysis with TB Green Premix Ex Taq II (Takara, Kyoto, Japan) following the manufacturer’s instructions. Relative expression was calculated using 2^−ΔΔCt method and normalized with GAPDH. The primer sequences of corresponding genes used in qRT-PCR analysis are listed in Table 1.

Statistical analysis

All statistical analyses were performed using GraphPad Prism 7.0 (GraphPad software, Inc., La Jolla, CA, USA) with Student’s t-test. All experiments in BMSCs were repeated independently at least three times. All data were represented as mean ± SEM. Values of p < 0.05 were considered statistically significant.

Identification of DEGs

We screened out 2,869 DEGs using GEO2R, including 1,357 up-regulated genes and 1,512 down-regulated genes, based on the dataset GSE107789. Corresponding volcano plot is shown in Fig. 1 and heat map of top 100 DEGs in Fig. 2.

GO terms analysis

To get a better insight into how DEGs are orchestrated in the process of BMSC adipogenic differentiation, DAVID database was used to perform GO terms enrichment analysis. For Biological Functions (BP), we identified 38 GO terms. The enriched terms were mostly related to cell cycle and mitosis, such as cell division, DNA replication, mitotic nuclear division, G1/S transition of mitotic cell cycle and DNA replication-dependent nucleosome assembly (Fig. 3A). For Cell Content (CC), DEGs were enriched in 21 terms while the top terms were mostly related to chromosome and extracellular content, such as chromosome space, nuclear chromosome and extracellular region. For MF, we identified seven GO terms, including protein binding, protein heterodimerization activity, heparin binding, single-stranded DNA-dependent ATPase activity, growth factor activity, cytoskeletal protein binding, histone binding.

Gene Ontology terms were further analyze for the top 200 DEGs. As shown in Fig. 3B, the top 10 BP terms enriched were skeletal system development, regulation of cell proliferation, response to lipopolysaccharide, aging, female pregnancy, brown fat cell differentiation, inflammatory response, positive regulation of cell migration, cell adhesion and extracellular matrix organization.

KEGG pathway analysis

DAVID database was used to obtain KEGG pathway enrichment. DEGs count and enriched KEGG pathway of all DEGs are shown in Fig. 4A. The top five pathways were DNA replication, cell cycle and systemic lupus erythematosus, mismatch repair and Alcoholism. In addition, KEGG pathway analysis of top 200 DEGs was performed (Fig. 4B). The top five pathways enriched were pathways in cancer, tyrosine metabolism, drug metabolism-cytochrome P450, PPAR signaling pathway and PI3K-Akt signaling pathway.

GSEA enrichment analysis

A total of 19 significant enriched pathways were identified by GSEA (as shown in Table 2), three of which were closely associated with adipogenesis or lipogenesis (Figs. 5A–5C). Furthermore, leading edge analysis was conducted of these three pathways (Figs. 5D and 5E) in which 3 of 51 genes, CPT1A, ACSL1 and ACSL3, appeared in all three pathways (Fig. 5D). Most importantly, ACSL1 had the highest transcriptional level in adipogenic cells compared to CPT1A and ACSL3 (Fig. 5E).

PPI network analysis: hub proteins and PPI module

Protein–protein interactions network was obtained using STRING database. In consideration of a huge number of DEGs, we excluded out of the PPI network analysis the DEGs whose |logFC| was less than 2, leaving 1,668 DEGs for further analysis. As shown in Fig. 6A, 1,481 nodes and 21,274 edges were included in the network. Among proteins constructing the present network, those with high node degrees were considered as hub proteins. In this case, Cyclin-dependent kinase 1 (CDK1), Cyclin-A2 (CCNA2), G2/mitotic-specific cyclin-B1, aurora kinase A (AURKA), aurora kinase B (AURKB), CDC20, PLK1, BUB1, TOP2A, interleukin 6 (IL6) were the top 10 hub proteins. In pursuit of a better understanding of the results of PPI network analysis, MCODE plugin of Cytoscape was introduced to screen out top PPI modules. The top three modules are shown in Figs. 6B–6D, respectively. Module 1 of the present PPI network was consisted of 142 nodes and 8,251 edges. Most of the proteins clustered in module 1 were down-regulated during adipogenic differentiation. The top five proteins in top three modules are shown in Table 3. The top five proteins with higher degrees in module 1 were CDK1, CCNA2, CCNB1, AURKA and AURKB. Module 2 was consisted of 59 nodes and 678 edges, including interleukin 8 (IL8), Stromal cell-derived factor 1 and sphingosine 1-phosphate receptors 3 (S1PR3). Moreover, module 3, which comprised 50 nodes and 358 edges, contained IL6, bone morphogenetic protein 4 (BMP4), apolipoprotein B-100 (APOB), osteopontin (SPP1), zinc finger and BTB domain containing 16 (ZBTB16), glypican 3 (GPC3), etc.

Verification of candidate gene expression

Based on the above data, four candidate genes, S1PR3, ACSL1, ZBTB16 and GPC3, were further verified in adipogenic cells and OVX osteoporosis mice model. For adipogenic differentiation of BMSCs, ORO staining showed more lipid drops in cells undergoing adipogenic differentiation on day 7 (Figs. 7A and 7B). qRT-PCR analysis showed up-regulated expression of adipogenic marker genes (ADIPOQ, CEBPA, PPARG) (Fig. 7C), confirming the adipogenic differentiation of BMSCs. Next, we verified the expression of the candidate genes. Results showed significantly increased mRNA expression of ACSL1, S1PR3 and ZBTB16 (Fig. 7D).

Further, OVX mice model was used to verify the alterations of the above four candidate genes in vivo. ORO staining showed more lipid drops in tibias from OVX mice (Figs. 7E and 7F), confirming the increased adipogenesis. qRT-PCR confirmed up-regulated expression of ACSL1, S1PR3 and ZBTB16 in bone (Fig. 7G). Together, the above data suggest intrinsic connections of ACSL1, S1PR3 and ZBTB16 with adipogenesis.