Gene Network Analysis for Osteoporosis, Sarcopenia, Diabetes, and Obesity in Human Mesenchymal Stromal Cells

The systemic gene interactions that occur during osteoporosis and their underlying mechanisms remain to be determined. To this end, mesenchymal stromal cells (MSCs) were analyzed from bone marrow samples collected from healthy individuals (n = 5) and patients with osteoporosis (n = 5). A total of 120 osteoporosis-related genes were identified using RNA-sequencing (RNA-seq) and Ingenuity Pathway Analysis (IPA) software. In order to analyze these genes, we constructed a heatmap of one-way hierarchical clustering and grouped the gene expression patterns of the samples. The MSCs from one control participant showed a similar expression pattern to that observed in the MSCs of three patients with osteoporosis, suggesting that the differentiating genes might be important genetic determinants of osteoporosis. Then, we selected the top 38 genes based on fold change and expression, excluding osteoporosis-related genes from the control participant. We identified a network among the top 38 genes related to osteoblast and osteoclast differentiation, bone remodeling, osteoporosis, and sarcopenia using the Molecule Activity Predictor program. Among them, 25 genes were essential systemic genes involved in osteoporosis. Furthermore, we identified 24 genes also associated with diabetes and obesity, among which 10 genes were involved in a network related to bone and energy metabolism. The study findings may have implications for the treatment and prevention of osteoporosis.


Introduction
Osteoporosis is a common, but severe, systemic skeletal disease caused by low bone mass, which results in bone fragility and results in increased risk of hip, spine, wrist, and other bone fractures [1,2]. Risk factors such as aging, low body weight, smoking, alcoholism, and a history of fracture contribute to osteoporosis development [3]. Osteoporosis is characterized by reduced bone mineral density (BMD), deteriorated bone tissue microarchitecture, and a low amount and lack of variety in altered bone proteins [4]. However, systemic gene interactions and molecular mechanisms underlying osteoporosis remain unelucidated.
Obesity is a global health concern. In 2016, over 650 million adults worldwide were classified as obese, with high-energy food intake and physical inactivity considered as the main causes of obesity and overweight [5]. Numerous studies have reported that certain genes affect the progression of obesity, including peroxisome proliferator-activated receptor γ (PPARG), CCAAT enhancer binding protein α (C/EBPα), and interleukin 6 (IL6) [6][7][8]. Moreover, adiponectin, a protein secreted by adipocytes that is involved in regulating glucose levels and lipid metabolism, can be regulated by osteoblasts [9]. Mesenchymal stromal cells (MSCs) can be isolated from various tissues and differentiated into multiple mesenchymal families, such as osteogenic, chondrogenic, adipogenic, and neurogenic cells [10,11]. In addition, the transcriptional profile of osteoporotic stem cells contains numerous genes that predispose individuals to osteoporosis [12]. Osteogenesis and adipogenesis are controlled by the Wnt signaling pathway [13], and nuclear receptor subfamily 4A (NR4A) is commonly associated with osteoporosis and obesity [14].
Osteoporosis is a common complication of diabetes mellitus, and patients with diabetes usually have skeletal disorders [15][16][17]. Diabetic osteopathy, an important comorbidity of both type-1 and type-2 diabetes mellitus, increases the risk of bone fracture [18]. Mutation in the gene encoding lipoprotein receptor-related protein 6 (LRP6), a co-receptor in the Wnt signaling pathway, is associated with osteoporosis, diabetes mellitus, and coronary artery disease [19].
Although previous studies have investigated the individual mechanisms responsible for the development of osteoporosis, obesity, and diabetes, the common molecular mechanisms underlying the development of these conditions have not yet been determined. In addition, systematic genetic analyses of endocrine diseases need to be performed. Therefore, the aim of the current study was to identify the genetic network related to osteoporosis using next generation sequencing to analyze MSCs isolated from skeletal bone marrow samples collected from healthy individuals with normal bone mineral density (BMD) and patients with osteoporosis. Subsequently, we aimed to discover the network comprised of differentially expressed genes (DEGs) related to osteoporosis, diabetes, and obesity using Ingenuity Pathway Analysis, thus elucidating the relationship of these genes with osteogenesis and adipogenesis.

Study Subjects
Ten subjects aged 40 to 85 years, regardless of gender and with and without osteoporosis, were enrolled in the present study. MSCs were obtained from samples of femur bone marrow collected from participants during total hip or knee arthroplasty due to osteoarthritis, hip fracture, or traffic-related accident. Subjects who presented with other metabolic disorders and secondary causes of osteoporosis were excluded. BMD was measured using dual-energy X-ray absorptiometry (iDXA, GE Lunar, Madison, WI, USA). Subjects with osteoporosis had vertebral fractures or hip fractures and low BMD (T-scores of the lumbar spine, femoral neck, and total hip were ≤−2.5). The five subjects in the control group had no prior fractures and normal BMD (T-scores of the lumbar spine, femoral neck, total hip, and whole BMD were ≥−1.0). The obesity criterion was BMI ≥ 25 kg/m 2 based on the guidelines of the Korean (Asian) Society for the Study of Obesity [20]. Subjects with osteoporosis and type-2 diabetes mellitus were also enrolled in the study if the onset of diabetes was >30 years of age and supported by medical history/records.

Human MSC Isolation and Culture
Human MSCs were isolated by negative immunoselection (RosetteSep Isolation Kit, STEMCELL Technologies Inc., Vancouver, BC, Canada) according to manufacturer's instructions. Bone marrow samples were incubated at room temperature (20-22 • C) for 25 min with a depletion cocktail of tetrameric antibodies. The samples were then diluted in 1× phosphate-buffered saline and isolated by density-gradient centrifugation (CR3-22, Jouan, Saint Herblain, France). Isolated human MSCs were seeded into vented 25 cm 2 tissue culture flasks (Thermo Fisher Scientific, Waltham, MA, USA) containing Complete MesenCult ® Medium (human) and MesenCult ® MSC Basal Medium (human) mixed with Mesenchymal Stem Cell Stimulatory Supplements (human) (STEMCELL Technologies Inc.) and streptomycin (Welgene, Gyeongsan, Korea). The flasks were incubated in a humidified atmosphere containing 5% CO 2 at 37 • C, and half of the medium was replaced with fresh medium every 3 days. When cells reached 90% confluence, the adherent cells were subcultured in a 100 mm 2 culture dish (first passage) and re-seeded into three 100 mm 2 culture dishes (second passage). Once cells reached confluence at the end of the second or third passage, they were used immediately or cryo-preserved [21].

RNA Isolation
Total RNA content was extracted from second or third passaged human MSCs using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions and quantified using a Nanadrop1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The resulting complete RNA was used for RNA-sequencing.

RNA Library Preparation and Sequencing
To construct cDNA libraries, 1 µg of total RNA was used with the TruSeq strandspecific mRNA library kit (Illumina Inc., San Diego, CA, USA) to achieve polyA-selected RNA extraction, RNA fragmentation, random hexamer primed reverse transcription, and 100 nt paired-end sequencing using the HiSeq 4000 sequencing system (Illumina Inc.). The cDNA libraries were quantified by using quantitative polymerase chain reaction (qPCR) according to the qPCR Quantification Protocol Guide, and library quality was assessed using the 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA).
To estimate expression levels, RNA-seq reads were mapped to the human genome using TopHat version 1.3.3 software (accessed on 2 August 2018) [22]. The human reference genome sequence (hg19) and annotation data were downloaded from the University of California Santa Cruz Genome Browser website (http://genome.uscs.edu, accessed on 2 August 2018). The transcript counts were calculated at the gene level, the relative transcript abundance was reported in terms of fragments per kilobase of exon per million (FPKM), and fragments were mapped using Cufflinks version 1.2.1 software [23].

Analysis of Gene Expression
Relative gene abundance was calculated in Read Count using StringTie. Statistical analysis was performed to identify DEGs using the estimates of abundance for each gene in the sample. Sample genes that had Read Count values of one more than zero were excluded. The statistical significance of the DEG set was assessed using the edgeR package in Bioconductor and fold change with a null hypothesis showed that no difference existed between groups.

Hierarchical Clustering
Hierarchical clustering analysis was performed using complete linkage and Euclidean distance as a measure of similarity to show the expression patterns of DEGs that satisfied the conditions of |fold change| ≥ 2 and raw p < 0.05.

Network Analysis
To construct pathways between genes of interest and biofunctions, we used the Path Explorer tool in Ingenuity Pathway Analysis (IPA) software (QIAGEN Inc., Hilden, Germany) (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis, accessed on 2 June 2021). Pathways were connected by adding intermediate molecules from the QIAGEN Knowledge Base. The Molecule Activity Predictor tool was used to identify downstream effects. Red and pink nodes indicated increased expression of DEGs, while green nodes indicated reduced expression. Orange nodes indicated predicted gene activation and inhibition, whereas blue nodes indicated biofunctions and diseases. The degree of expression or prediction was indicated by node color intensity. Solid lines indicated direct interaction, whereas dashed lines suggested possible gene interaction.

Statistical Analysis
We performed statistical analyses with R software, version 4.0.1 (R Foundation for Statistical Computing, Vienna, Austria). Participant data of clinical characteristics were converted to standardized score (Z-score) and then calculated by the univariate linear regression model with age. To translate the data converted to Z-score to the original scale, an age-adjusted value was multiplied by the sample standard deviation of raw data, and the average value of raw data was added.

Clinical Characteristics of Study Participants
We compared MSCs extracted from the bone marrow of five patients with osteoporosis (osteoporosis group) with those extracted from the bone marrow of five healthy individuals with normal BMD (control group) (Table 1 (a). The mean weights and BMI values of the osteoporosis group were also lower than those of the control group. Additionally, most participants in the osteoporosis group were diagnosed with type 2 diabetes. There was a significant difference between control and osteoporosis patients with lumbar spine, femur neck, total hip, and whole-body BMD values in age-adjusted analysis data (Table 1 (b)). Other variables including sex, diabetes, obesity, and menopause did not show any significance.

Comparison of Gene Expression Patterns between Control and Osteoporosis Groups
To determine the essential systemic genes related to osteoporosis, DEG analysis was performed on MSCs isolated from each participant. A total of 120 DEGs were identified, and a heatmap of one-way hierarchical clustering was constructed ( Figure 1A). Subsequently, 120 DEGs were regrouped by similarity of expression pattern via hierarchical clustering analysis of significant genes. According to the heatmap, three participants in the osteoporosis group (no. 6, 8, and 9) demonstrated similar gene expression patterns, which were similar to that displayed by participant no. 3 in the control group ( Figure 1B). This result suggested that the differentiating genes between control subject no. 3 and those commonly upregulated and downregulated (i.e., similar expression patterns) among the other four participants in the control group could be used as genetic factors to classify participant no. 3 as the control subject in further analysis.

RNA Analysis of Bone Remodeling Epigenome
To compare the gene expression pattern observed for control subject no. 3 with that of the other participants in the control (excluding the chemical molecules registered in IPA), nine genes related to bone remodeling were selected among the top 120 genes (marked in red in Figure 1A). To connect these nine genes with the bone remodeling epigenome, intermediate molecules expressed in the MSCs of all participants in the control group were added to produce 38 genes. The expression levels of these 38 genes for the participants in the control group are listed in Tables 2 and S1. The expression of some genes was upregulated for control subject no. 3 compared to that of the other control group participants, while that of other genes was downregulated (Table S2). Specifically, the expressions of BCL2 apoptosis regulator (BCL2), bradykinin receptor B2 (BDKRB2), insulin-like growth factor 1 (IGF1), IL6, PPARG coactivator 1 α (PPARGC1A), transferrin (TF), and 19 other genes were upregulated for control subject no. 3 compared to that of other control group participants, while the expressions of insulin-like growth factor 2 (IGF2), interleukin 1 receptor type 1 (IL1R1), and PPARG were downregulated. In contrast, the expression of myocyte enhancer factor 2C (MEF2C) was upregulated for control subject no. 3 but downregulated for the other participants in the control group (Table 3). Direct and indirect relationship networks of the 38 genes related to osteoblasts, osteoclasts, bone remodeling, osteoporosis, and sarcopenia were subsequently generated (Figures 2 and S1). Based on these data, we identified 25 commonly upregulated genes that might be essential systemic genes involved in normal bone remodeling and osteoporosis prevention.

Genetic Relationship among Osteoporosis, Obesity and Diabetes
To clarify the related pathways among osteoporosis, obesity, and diabetes, the DEGs in obesity-related and diabetes-related networks were identified among the MSCs of all participants in the control group genes using IPA software. A total of 24 osteoporosisrelated genes, including IGF1, IL6, pyruvate dehydrogenase kinase 4 (PDK4), PPARGC1A, and TF, were also associated with obesity and diabetes. All 24 genes were highly expressed with direct or indirect relationships established via intermediate molecules such as calcitriol, colony stimulating factor 2 (CSF-2), interferon-γ (IFN-γ), IL1, IL4, transforming growth factor-α (TGF-α), transforming growth factor-β (TGF-β), and tumor necrosis factor (TNF) (Tables 4 and S3). The relationship networks are illustrated in Figures 3 and S2, and they demonstrate that 24 candidate genes reduced diabetes and increased obesity. These results suggested that the 24 candidate genes might be involved in the development of osteoporosis as well as obesity and diabetes.

Genetic Relationship among Osteoporosis, Obesity and Diabetes
To clarify the related pathways among osteoporosis, obesity, and diabetes, the DEGs in obesity-related and diabetes-related networks were identified among the MSCs of all participants in the control group genes using IPA software. A total of 24 osteoporosisrelated genes, including IGF1, IL6, pyruvate dehydrogenase kinase 4 (PDK4), PPARGC1A, and TF, were also associated with obesity and diabetes. All 24 genes were highly expressed with direct or indirect relationships established via intermediate molecules such as calcitriol, colony stimulating factor 2 (CSF-2), interferon-γ (IFN-γ), IL1, IL4, transforming growth factor-α (TGF-α), transforming growth factor-β (TGF-β), and tumor necrosis factor (TNF) (Tables 4 and S3). The relationship networks are illustrated in Figures 3 and S2, and they demonstrate that 24 candidate genes reduced diabetes and increased obesity. These results suggested that the 24 candidate genes might be involved in the development of osteoporosis as well as obesity and diabetes.  Entrez gene names and p-values are listed in Table S3.

Common Related Genes in Osteoporosis, Sarcopenia, Diabetes, and Obesity
The IPA results revealed 10 common candidate genes (CFB, CXCL2, HSD11B1, IGF1, IL6, NTN1, PCSK1, PPARGC1A, SOD2, and TF) that were inhibited in osteoporosis, sarcopenia, and diabetes but activated in obesity for all participants in the control group ( Figure 4). Furthermore, IPA data suggested that activation of HSD11B1 and PPARGA1A resulted in obesity, but the other eight candidate genes were not significantly associated with obesity.

Common Related Genes in Osteoporosis, Sarcopenia, Diabetes, and Obesity
The IPA results revealed 10 common candidate genes (CFB, CXCL2, HSD11B1, IGF1, IL6, NTN1, PCSK1, PPARGC1A, SOD2, and TF) that were inhibited in osteoporosis, sarcopenia, and diabetes but activated in obesity for all participants in the control group ( Figure 4). Furthermore, IPA data suggested that activation of HSD11B1 and PPARGA1A resulted in obesity, but the other eight candidate genes were not significantly associated with obesity.

Discussion
MSCs are multipotent cells with the potential to differentiate into multiple lineages, including osteoblasts, adipocytes, and chondroblasts [24,25]. Therefore, MSCs were employed in the present study to investigate the relationship between bone and fat cells. RNA-seq was performed to identify DEGs in MSCs collected from healthy individuals and subjects with osteoporosis [26][27][28][29]. Previous studies have pooled numerous samples together to obtain DEG data during investigations of the overall trends of disease [14,30,31]. In the traditional pooling system, individual specific genes may not be selected or lost by other control samples. In the current study, each of the five samples from the control and osteoporosis groups was analyzed separately to determine the minimum individual character gene groups involved in osteoporosis.
Using IPA software, we identified 25 important systemic genes related to bone remodeling and 24 candidate osteoporosis-related genes associated with obesity and diabetes. Moreover, 10 common genes that influence metabolic disorders such as osteoporosis,

Discussion
MSCs are multipotent cells with the potential to differentiate into multiple lineages, including osteoblasts, adipocytes, and chondroblasts [24,25]. Therefore, MSCs were employed in the present study to investigate the relationship between bone and fat cells. RNA-seq was performed to identify DEGs in MSCs collected from healthy individuals and subjects with osteoporosis [26][27][28][29]. Previous studies have pooled numerous samples together to obtain DEG data during investigations of the overall trends of disease [14,30,31]. In the traditional pooling system, individual specific genes may not be selected or lost by other control samples. In the current study, each of the five samples from the control and osteoporosis groups was analyzed separately to determine the minimum individual character gene groups involved in osteoporosis.
Using IPA software, we identified 25 important systemic genes related to bone remodeling and 24 candidate osteoporosis-related genes associated with obesity and diabetes. Moreover, 10 common genes that influence metabolic disorders such as osteoporosis, sarcopenia, diabetes, and obesity were identified. Most of the candidate genes were related to bone remodeling and/or sarcopenia, diabetes, and obesity, as supported by the results of previous studies. Complement factor B (CFB) plays a significant role in the pathogenesis of diabetic retinopathy [32], and there might be a relationship between CFB and diabetes onset. Furthermore, CFB plays an important role in the differentiation of 3T3-L1 preadipocytes in CFB transgenic mice [33], and CFB and IL6 are associated with the rapid metabolism of glucagon in metabolically compromised subjects [34]. Additionally, IL6 individually controls osteoblast differentiation and induces bone resorption [35]. Since hypertrophied white adipose tissue produces chemokines, CXCL2 expression is increased in osteoclastogenesis [36] compared with osteoblastogenesis [37], and white adipose tissue is activated by CXCL2 to prevent the onset of obesity [38]. HSD11B1 might be related to BMD [39] and is associated with the onset of type 2 diabetes [40,41]. IGF-1 is known to enhance muscle mass and alleviate sarcopenia [42], and low serum IGF-1 levels could increase the risk of idiopathic osteoporosis [43,44]. Low IGF-1 levels not only affect BMD but also results in the onset of type 1 or 2 diabetes mellitus [45,46]. NTN1 is an essential factor involved in the differentiation of osteoclasts [47] and is also involved in osteoblast differentiation [48]. Low PCSK1 expression might be related to obesity and diabetes [49], and PPARGC1A expression can increase obesity risk [50]. SOD2 increases osteoblast differentiation and bone formation [51,52] and protects neurons in diabetic neuropathy [53]. The expression of SOD2 is higher in patients with obesity than in healthy individuals [54]. TF can repress alkaline phosphatase activities in osteoblasts due to high oxidative stress [55], and increased TF expression increases the risk of type 1 and/or type 2 diabetes [56]. Taken together, these findings support that the 10 candidate genes identified in the current study are associated with metabolic disorders.
Osteoporosis commonly occurs with other endocrine system-related diseases, such as obesity and diabetes. Bone and fat cells are known to have a common origin [57], and osteoporosis and diabetes can both be induced by vitamin D deficiency [58]. Moreover, insulin signaling pathway-associated proteins can induce both osteoporosis and diabetes [59]. However, studies that include osteoporosis, sarcopenia, obesity, and diabetes are extremely rare.
In the current study, we analyzed individual differences in MSC gene expression among participants in the control and osteoporosis groups to overcome the disadvantage of previous studies that investigated overall expression trends. Another merit of the present study is that using MSCs to identify osteoporosis-related DEGs enabled us to evaluate the genetic relationship among osteoporosis, sarcopenia, obesity, and diabetes using cells from which osteoblasts and adipocytes commonly originate. Nevertheless, this study has some limitations. First, the mean age difference between the control and osteoporosis groups is statistically significant. The subjects in the control group underwent surgery due to unexpected accidents; therefore, no medical history was available and the statistical difference was unavoidable. Second, we analyzed transcriptional differences associated with osteoporosis-related genes (mRNA). The process of mRNA translation is not always successful; thus, a mechanistic study of these genes is needed. The molecular mechanisms of the identified genes will be explored in future studies.
In summary, we identified the top 25 essential systemic genes involved in osteoporosis using IPA software and the Molecule Activity Predictor tool. Furthermore, we constructed a network of the top 24 genes also associated with diabetes and obesity to verify their association with osteoporosis. The study findings provide insight into the relationship among osteoporosis, obesity, and diabetes, suggesting potential targets for osteoporosis treatment and prevention.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes13030459/s1, Figure S1: Functional relationship networks of genes associated with the canonical pathways related to osteoporosis development in control group participants; Figure S2: Functional relationship networks of genes associated with the canonical pathways related to obesity and diabetes onset in control group participants; Table S1: RNA-sequencing data of 38 genes related to osteoblasts, osteoclasts, bone remodeling, osteoporosis, and sarcopenia in mesenchymal stromal cells among participants in the control group; Table S2: Genes commonly upregulated or downregulated for control subject no. 3 according to RNA-sequencing data, compared to control participants no. (A) 1, (B) 2, (C) 4, and (D) 5; Table S3: Expression of 24 genes related to obesity and diabetes in mesenchymal stromal cells among participants in the control group based on RNA-sequencing data.