STRA6 as a possible candidate gene for pathogenesis of osteoporosis from RNA‑seq analysis of human mesenchymal stem cells.

To identify novel candidate genes associated with osteoporosis, RNA‑sequence analysis of human mesenchymal stem cells (hMSCs) from patients with osteoporosis (G3) and osteopenia (G2), and healthy controls (G1) was performed. Differentially expressed genes (DEGs) from among the three groups were identified. DEGs were separated into nine groups according to their gene expression patterns: UU (up and up), UF (up and flat), UD (up and down), FU (flat and up), FF (flat and flat), FD (flat and down), DU (down and up), DF (down and flat), and DD (down and down). Among the 42 DEGs between G3 and G1, eight candidate genes, namely stimulated by retinoic acid 6 (STRA6), melanophilin, neurotrophic receptor tyrosine kinase 2, cartilage oligomeric matrix protein, collagen type XI α 1 chain, integrin subunit β 2, monooxygenase DBH‑like 1 and selenoprotein P, were selected, as they demonstrated consistent gene expression patterns of UU, FU, FD, and DD. Among these eight genes, STRA6 was highly expressed in the osteoporosis group and based on additional data from quantitative polymerase chain reaction analysis, it was selected for further study. In order to investigate whether STRA6 served a functional role in osteoblast or adipocyte differentiation, the effects of STRA6 expression changes in pluripotent stem cell C3H10T1/2, preosteoblast MC3T3‑E1 and stromal ST2 cell lines were examined. Bone morphogenetic protein 2 enhanced STRA6 expression only at the early stage of osteoblast differe-ntiation, and overexpression of STRA6 temporally inhibited the expression of osteoblastogenesis markers, including runt related transcription factor 2, bone sialoprotein and osteocalcin. Furthermore, the knockdown of STRA6 slightly enhanced nodule formation at the late stage of osteoblast differentiation, and overexpression of STRA6 in ST2 cells enhanced adipocyte differentiation. Taken together, STRA6 expression could be associated with the pathogenesis of osteoporosis by promoting adipocyte differentiation over osteoblast differentiation in the hMSC population.

Introduction

Osteoporosis is characterized by skeletal fragility due to low bone mineral density (BMD). It leaves the afflicted susceptible to low-trauma fractures, particularly among the elderly and is a major public health concern (1–2). Although, susceptibility to fractures depends on many factors, including non-skeletal factors such as propensity to fall, diminished soft-tissue cushioning, and more broadly on the physical environment, genetic factors indeed play a major role in the determination of BMD and osteoporosis risk (3). For example, a family history of hip fractures carries a two-fold increased risk of fracture among the descendants. Also, approximately 60% of peak bone mass is genetically determined (3,4). However, the majority of causal genes and molecular mechanisms have yet to be identified.

With the advent of whole-transcriptome sequencing using next-generation sequencing (NGS) technologies and RNA sequencing (RNA-Seq) with massively parallel sequencing, the transcriptome analyses of genomes at a far higher resolution is available; this is compared with previously available Sanger sequencing and microarray-based methods (5). NGS has several advantages over microarrays. In addition to accurate quantitation of RNA expression, NGS has allowed for detection of rare genes, analyses of transcript variants, and discovery of novel transcripts (6,7).

Bone strength is an integration of bone density and bone quality, and BMD is the most important factor for bone strength, accounting for about 70% of bone strength (8). Based on the World Health Organization (WHO) diagnostic classification, BMD categories are: i) normal bone mass (T-score of −1.0 and above), ii) osteopenia (T-score between −1.0 and −2.5), and iii) osteoporosis (T-score of −2.5 and below) (9). However, there have only been few studies comparing gene and expression profiles according to BMD classes (normal bone mass, osteopenia and osteoporosis) among various populations.

Mesenchymal stem cells (MSCs) from bone marrow are pluripotent cells that can differentiate into multiple cell lineages, including osteoblasts, adipocytes, fibroblasts and chondrocytes (10). Although MSCs are the origin of osteoblasts, it is currently not well understood if the individual intrinsic genetic differences in these cells contribute to the pathogenesis of osteoporosis (11). Therefore, RNA-seq information from MSCs of patients with osteoporosis may provide novel insights into the mechanisms underlying the pathogenesis of osteoporosis.

In the present study, we isolated human MSCs from the bone marrow of typical severe osteoporosis patients, osteopenia patients and normal BMD controls and identified DEGs among the groups using NGS technology. The screened DEGs were further analyzed and STRA6 gene was selected as a primary candidate for a replication study. We demonstrated experimentally the functional role of the STRA6 in osteoblast and adipocyte differentiation.

Materials and methods

Preparation of hMSC

hMSCs were obtained from the bone marrow of femurs during total hip or knee arthroplasty due to osteoarthritis or hip fracture; the study protocol was approved by the Ajou University Hospital Institutional Review Board (IRB; AJIRB-GEN-GEN-11-332). All hMSCs were obtained from postmenopausal women with an average age of 68.5 years (60–77 years old). Densitometric examinations were performed using a Lunar iDXA apparatus (GE Lunar, Madison, WI, USA). Diagnosis of osteopenia or osteoporosis were made using the WHO T-score criteria (−2.5 osteoporosis group had vertebral fractures.

Separation of hMSCs

Human MSCs were isolated using an RosetteSep Isolation kit (STEMCELL Technologies Inc., Vancouver, Canada) and according to the manufacturer's instructions. The cells were cultured at 37°C in a humidified atmosphere containing 5% CO2 with Complete MesenCult® Medium (human) and MesenCult® MSC Basal Medium (human) containing Mesenchymal Stem Cell Stimulatory Supplements (human) (STEMCELL Technologies) and antibiotics (Welgene, Gyeongsan, Korea). Half of the culture medium was replaced once a week with fresh medium. At the end of the second passage, as the cells reached confluence, they were trypsinized and were either cryopreserved or used immediately.

RNA isolation

Total RNA was isolated using an RNeasy Micro kit (Qiagen, Valencia, CA, USA) and treated with a Turbo DNA-Free™ kit (Life Technologies, Grand Island, NY, USA) to remove potential contaminating DNA that could lead to false-positive amplification. RNA quality and purity was confirmed with a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA, USA). The resulting intact RNA was used for either RNA-seq analysis or QPCR analyses.

RNA sequencing (RNA-seq)

RNA-seq was carried out at TheragenEtex (Suwon, Korea). Prior to cDNA library construction, quality was evaluated electrophoretically using an Agilent High Sensitivity DNA kit (Agilent Technologies, Santa Clara, CA, USA). The fragment sizes were 350–450 bps. Subsequently, the library was sequenced using a HiSeq 2500 sequencer (Illumina, San Diego, CA, USA) in rapid-run mode. Cluster generation was then performed, followed by 2×100 cycle sequencing reads, separated by a paired-end turnaround. Image analysis was performed using the HiSeq control software, version 1.8.4.

Reagents and constructs

STRA6 was prepared by RT-PCR using RNA from C3H10T1/2 cells. The primer sequences were as follows: STRA6 sense, 5′-CGGGAATTCATGGAGTCCCAGGCATCTCGG0-3′; STRA6 antisense, 5′-ATACTCGAGTATGGGCTGGGTGCCATTGGCC-3′. The amplified PCR fragments were cloned into the HA-pcDNA3.1 vector. ShRNA-STRA6 for knockdown was obtained from Sigma-Aldrich (St. Louis, MO, USA).

Semi-quantitative reverse transcription polymerase chain reaction (RT-PCR) analysis

C3H10T1/2 and MC3T3-E1 cells were induced by 50 µg/ml ascorbic acid, 10 mM beta-glycerophosphate, and 100 ng/ml BMP-2 (C3H10T1/2 only) to evaluate the expression of osteogenic marker genes during osteoblast differentiation. The PCR primers were as follows: STRA6 sense, 5′-GAACCTCACTGCAGCATCAA-3′; STRA6 antisense, 5′-GCTTTCTGCCAGACCTGAAC-3′; Runx2 sense, 5′-CCCAGCCACCTTTACCTACA-3′; Runx2 antisense, 5′-CAGCGTCAACACCATCATTC-3′; Bsp sense, 5′-AAAGTGAAGGAAAGCGACGA-3′; Bsp antisense, 5′-GTTCCTTCTGCACCTGCTTC-3′; OCN sense, 5′-GCGCTCTGTCTCTCTGACCT-3′; OCN antisense, 5′-ACCTTATTGCCCTCCTGCTT-3′; Sox9 sense, 5′-AAGGAAGCTGGCAGACCAGTA-3′; Sox9 antisense, 5′-CCCTCTCGCTTCAGATCAAC-3′; β-actin sense, 5′-AAGAGCTATGAGCTGCCT-3′; and β-actin antisense, 5′-CACAGGATTCCATACCCA-3′.

Cell culture and assay of ALP and nodule formation and Oil Red O staining

To overexpress or knockdown STRA6, MC3T3-E1 and ST2 cells were transfected with HA-pcDNA3.1, HA-STRA6-pcDNA3.1 or ShRNA-STRA6 on 48-well plates at a density of 3×104, one day before transfection. Plasmid DNA was mixed with TransIT-2020 (Mirus, Madison, WI, USA) and transfected into the cells, as per manufacturer's protocol. After 48 h of transfection, the cells were induced with osteogenic media containing 50 µg/ml ascorbic acid and 10 mM β-glycerophosphate for 5 days ~2 weeks and culture medium was replaced every 3 days. Alkaline phosphatase (ALP) and nodule formation assay were performed as previously described (12). For the assay of adipocyte differentiation, transfected cells induced with adipogenic media including 0.5 mM IBMX, 5 µg/ml insulin and 1 µM dexamethasone for 3, 5 and 7 days and the culture medium was replaced every 3 days. To determine lipid accumulation of Stra6 induced cell, cells fixed with 4% formaldehyde were stained with 0.5% Oil Red O (Sigma) for 10 min and washed with PBS.

Statistical analysis

To compare the subject basal characteristics, statistical comparisons were performed using IBM SPSS Statistics 23 (SPSS Inc., Chicago, IL, USA) and the data are presented as mean ± SEM.

Results

Clinical characteristics of the study subjects for RNA-seq analysis

Clinical characteristics of each four sample of normal, osteopenia and osteoporosis patients group (60–77 years old) are shown in Table I. There were no statistically significant differences in basal characteristics such as age, height and body mass index (BMI) except for weight and BMD.

Table I.

Clinical characteristics of the study subjects for RNA-seq.

	Patient no.	Age (years)	Height (cm)	Body weight (kg)	L-spine T-score	Total hip T-Score	Femur neck T-score
	1	60	147.8	77.0	−0.8	0.0	−0.4
	2	70	153.1	70.8	0.9	−0.3	−0.7
Normal	3	71	56.8	83.7	0.5	0.5	−0.2
	4	68	155.2	74.0	−0.9	0.0	−0.3
	AVG	67.3	153.7	76.4	−0.1	0.1	−0.4
	1	64	156.0	66.7	−1.7	−1.6	−2.0
	2	70	155.0	63.4	−1.4	−0.7	−0.7
Osteopenia	3	66	155.0	64.3	0.6	0.0	−1.2
	4	77	149.5	56.9	−1.5	−1.0	−1.6
	AVG	69.3	153.9	62.8	−1.0	−0.8	−1.4
	1	65	153.1	59.7	L1-L4 Fx	−2.7	−2.4
	2	73	154.0	73.8	−2.8 (L3-L4 Fx)	−1.8	−2.4
Osteoporosis	3	74	150.0	58.3	−4.1 (L1 Fx)	−2.1	−2.7
	4	62	151.3	61.4	−1.1 (L3-L4 Fx)	−2.7	−2.4
	AVG	68.5	152.1	63.3	−2.7	−2.3	−2.5
	P-value	0.912	0.535	0.039	0.071	0.012	0.008

AVG, average; Fx, fracture.

Differentially expressed genes (DEGs) from the RNA-seq analysis among G1, G2 and G3

To understand the different gene expression patterns among three groups; normal (G1), osteopenia (G2) and osteoporosis (G3), gene expression patterns of G3 vs. G2, G3 vs. G1, and G2 vs. G1 were compared separately. Significantly differential genes for MSCs between the groups were defined as q-values <0.05 and log2 FC <-2. When comparing G3 with G2, 47 DEGs were up-regulated and 10 DEGs were down-regulated (Table II). Comparing G3 with G1, 20 DEGs were up-regulated and 22 DEGs were down-regulated (Table II). For the G2 with G1 comparison, 16 DEGs were up-regulated and 49 DEGs were down-regulated (Table II). According to the gene expression patterns within three groups, we classified the gene data into 9 groups (P-value <0.001): UU (up and up), UF (up and flat), UD (up and down), FU (flat and up), FF (flat and flat), FD (flat and down), DU (down and up), DF (down and flat), DD (down and down) (Fig. 1).

Table II.

DEGs and Gene ontology (GO) classification of MSCs between osteoporosis (G3) vs. osteopenia (G2) vs. normal (G1).

No.	Name	DEG	GO	Molecular function (%)	Biological process (%)	Cellular component (%)
DEG-1	G3 vs. G2	57	3791	14/484 (2.9)	204/3005 (6.8)	18/302 (6.0)
		(Up 47: Down 10)
DEG-2	G3 vs. G1	42	4018	35//488 (7.2)	334/3192 (10.5)	55/338 (16.3)
		(Up 20: Down 22)
DEG-3	G2 vs. G1	65	4243	27/506 (5.3)	492/3424 (14.4)	28/313 (8.9)
		(Up 16: Down 49)

DEG, differentially expressed gene; GO, gene ontology.

Figure 1.

Expression patterns and variations for each group (cutoff: P-value <0.001).

Classification of gene ontology (GO) and analysis of gene function category

For gene and gene product grouping, we used gene ontology classification for categorizing gene function. The cluster for ‘molecular function’ (2.9%), ‘biological process’ (6.8%) and ‘cellular component’ (6.0%) show in 3791 GO cases for G3 vs. G2. Secondly, the cluster of 4018 GO cases of G3 vs. G1 shows ‘molecular function’ (7.29%), ‘biological process’ (10.5%) and ‘cellular component’ (16.3%). Thirdly, the cluster for ‘molecular function’ (5.3%), ‘biological process’ (14.4%) and ‘cellular component’ (8.9%) are in 4243 GO cases for G2 vs. G1 (cutoff: P-value <0.05) (Table II).

Selection of the candidate genes among the DEGs and confirmation of the gene expressions by real-time PCR

To find promising candidate genes, among 42 DEGs between G1 vs. G3, we selected 8 DEGs such as STRA6, Mlph, Ntrk2, Comp, Col11a1, ITGb2, Moxd1 and Sepp1, which had shown meaningful consistent gene expression patterns (UU, FU, FD and DD) (Fig. 1). Real-time PCR was performed to confirm the gene expression patterns for the selected 8 DEGs. Among them, STRA6 and Mlph were up-regulated in osteoporosis patients compare to normal and osteopenia cases. The others were down-regulated in osteoporosis patients. In addition, eight genes expression pattern were approximately consistent with the results of the RNA-seq data (Fig. 2). Finally, we selected STRA6 gene for further analysis for association of osteoblast differentiation and adipocyte differentiation, as it showed the highest expression in the osteoporosis group (Fig. 2).

Figure 2.

mRNA expression of candidate genes, STRA6, Mlph, Ntrk2, Comp, Col11a1, ITGb2, Moxd1 and Sepp1, in normal, osteopenia and osteoporosis from RNA seq. Data and real-time PCR. Results are representative of at least 3 independent sets of similar experiments. Data are presented as mean ± SD of triplicate samples. GAPDH was used as a reference. Group 1, normal; Group 2, osteopenia; Group 3, osteoporosis. *P<0.05 and **P<0.01 vs. normal.

Effects of Stra6 on osteoblastogenesis

BMP2 temporally enhanced STRA6 gene expression at an early time point in the osteogenic culture of C3H10T1/2 cells (Fig. 3A). Interestingly, the up-regulated STRA6 was accompanied by a slight down-regulation of osteogenic marker genes such as Runx2, Bsp, and OCN compared to Sox9 and β-actin in MC3T3-E1 cells. Overexpressed STRA6 gene, however, returned to near basal levels at 48 h (Fig. 3B). ALP and nodule formation were not significantly altered by overexpressed STRA6, while down-regulated STRA6 showed enhanced nodule formation at the late stage of osteoblast differentiation (P=0.06) (Fig. 3C, D).

Figure 3.

Effect of STRA6 on osteoblast differentiation and adipocyte differentiation. (A) BMP2 treatment of C3H10T1/2 for 5 days. (B) Overexpression STRA6 at MC3T3-E1 for 6 and 48 h. (C and D) Effect of overexpression (HA-pcDNA3.1-STRA6) and knockdown (shRNA-STRA6) of STRA6 for alkaline phosphatase (ALP) (day 5) and nodule formation (day 14) in MC3T3-E1. (E and F) Adipocyte differentiation by STRA6 in ST2 day 7. **P=0.01.

Effects of STRA6 on adipogenesis

ST2 cells were cultured for 7 days with adipogenic media and transfected by STRA6 to investigate the effects of STRA6 on adipocyte differentiation. (Fig. 3E, F). At the early stage (3rd day) of adipocyte differentiation, transfected cells temporally had enhanced adipocyte differentiation (P=0.01); however, this returned to near basal levels within 7 days.

Discussion

In this present study, we isolated hMSCs from bone marrow of severe osteoporosis patients, osteopenia patients and normal BMD controls and identified DEGs among the groups using the NGS technology. We identified 42 DEGs from the G3 (osteoporosis) vs. G1 (normal control) analysis and selected 8 candidate genes, which had consistent gene expression patterns within the osteoporosis, osteopenia and normal group with UU, FU, FD and DD characteristics. Among these 8 genes, STRA6 was highly expressed in the osteoporosis group and was selected for further study based on real-time PCR supporting data. From the functional study, STRA6 was shown as a possible candidate gene being involved in the pathogenesis of osteoporosis.

Initially, we investigated the DEGs among normal, osteopenia and osteoporosis through RNA-seq analysis. From this analysis, 57 DEGs between G3 and G2, 42 DEGs between G3 and G1 and 65 DEGs between G2 and G1 were found. From these GO's, various gene functions for these DEGs could be assigned. Bone diseases are from signaling pathways coupled to genetic as well as environmental factors (4,13). From our analysis, osteoporosis might involve various unidentified osteoporosis-related genes, which should be studied in detail in a further study.

For promising candidate genes from among the 42 DEGs from G1 vs. G3 analysis, we selected 8 DEGs, namely STRA6, Mlph, Ntrk2, Comp, Col11a1, ITGb2, Moxd1 and Sepp1, as they had consistent with RNA-seq data (Fig. 2). We selected STRA6 gene for further analysis, showing the highest expression in the osteoporosis group from real time PCR analysis.

STRA6, stimulated by retinoic acid 6 treatment, is well known as a membrane protein involved in the metabolism of retinol as it acts as a receptor for retinol/retinol binding protein complexes (14,15). In a study on craniosynostosis, suture mesenchymal cells expressed the RBP4 receptor, STRA6, on cell surfaces. In a cell culture model of cranial osteogenesis, RBP4 was significantly downregulated during mineralization (16). Moreover, genetic data show that single nucleotide polymorphisms in STRA6 were associated with type 2 diabetes (17), raising the possibility that STRA6 played an important role in altered metabolic states. STRA6 has also been proposed to play a role in insulin resistance, based on the data from a total-body STRA6 knockout mouse (18). Therefore, STRA6 might have a role in pathogenesis of osteoporosis, especially by involvement of differentiation of MSCs to osteoblasts or adipocytes. Therefore, we investigated whether STRA6 was involved in osteoblast and adipocyte differentiation.

We found that BMP2 enhanced STRA6 expression in osteogenic culture of C3H10T1/2 cells at an early time point. However, despite continued presence of BMP2 treatment, STRA6 expression was gradually down-regulated in a time-dependent manner. These results suggest that STRA6 could have an early role in osteoblastogenesis. Subsequently, we examined whether STRA6 could regulate osteogenic marker genes and morphological alteration in osteoblast culture of MC3T3-E1 cells. Interestingly, overexpressed STRA6 temporarily down-regulated appearance of osteogenic marker genes at 6 h; However overexpression STRA6 itself was decreased at 48 h post induction (Fig. 3B). Knockdown of STRA6 slightly enhanced nodule formation at the late stage of osteoblastogenesis. However exogenous overexpression of STRA6 showed no significant differences compared to control cells. These data indicate that STRA6 might act as a minor negative regulator, suppressing osteoblast differentiation at an early stage, and that it could be reversely regulated by other osteogenetic marker genes.

For effects of STRA6 on adipogenesis, STRA6 overexpression led to slight promotion of differentiation for adipocyte precursors. In summary, STRA6 might play a negative role in osteoblast differentiation and positive role in adipocyte differentiation at an early stage. Further studies for detailed mechanisms of STRA6 are required for elucidation of the pathogenesis of osteoporosis and additional therapeutic approaches for bone diseases.

The present study had some limitations. First, the sample size for expression profiling for each group may have appeared to be small. However, we tried to recruit subject representative of each group with overt differences in clinical characteristics in terms of osteoporosis to identify differential DEGs. We carefully recruited normal controls and osteopenia subjects without a medical history that would have affected bone metabolism. Other clinical characteristics of the study subjects such as age, height and BMI were not significantly different among our study groups. Moreover, the osteoporosis subjects were also cautiously selected. Subjects with a definite medical history of osteoporotic fracture were recruited for the osteoporosis group, besides having low BMD values (T-score ≤-2.5). Despite the small sample number, we identified distinct transcriptional differences among the groups with a reliable medical history. Second, our study was conducted in only a specific ethnic group. Therefore, our results cannot be generalized to other races or ethnicities. Actually, among the three study groups, the fragments per kilobase of exon per million fragments mapped (FPKM) of genes of the study subjects of RNA-seq appeared to have a strong correlation (r>0.966) (date was not shown). However, this genetic homogeneity could be rather a strong point of this study. With the careful selection of subjects for the study, the chance to find definite genetic differences for a specific disease could be higher than other study groups.

The present study also has some strengths. MSCs are the common origin of osteoblasts and adipocytes (10), and MSCs may be the best resource for research on the relationship between bone and fat. In this study, NGS technologies and RNA-Seq analysis were used to identify the DEGs, as these techniques have been known to have remarkable sensitivity and accuracy for expression profiling (19,20). Moreover, we performed the RNA-seq analysis for normal, osteopenia and severe osteoporosis subjects and investigated the gradual genetic differences according to the severity of bone loss. To our knowledge, there have only been a few studies performed with this study design to investigate a mechanism related to the pathogenesis of osteoporosis.

In conclusion, we identified 8 candidate DEGs, namely STRA6, Mlph, Ntrk2, Comp, Col11a1, ITGb2, Moxd1 and Sepp1, playing a possible role in pathogenesis of osteoporosis using the RNA-seq analysis from the hMSCs of subjects from osteoporosis and osteopenia and normal controls. Among them, STRA6 could be related to osteoporosis, favoring adipocyte differentiation over osteoblast differentiation in hMSCs.