Osteoblasts are key players in bone remodeling. The accessibility of human primary osteoblast-like cells (HObs) from bone explants makes them a lucrative model for studying molecular physiology of bone turnover, for discovering novel anabolic therapeutics, and for mesenchymal cell biology in general. Relatively little is known about resting and dynamic expression profiles of HObs, and to date no studies have been conducted to systematically assess the osteoblast transcriptome. The aim of this study was to characterize HObs and investigate signaling cascades and gene networks with genomewide expression profiling in resting and bone morphogenic protein (BMP)-2- and dexamethasone-induced cells. In addition, we compared HOb gene expression with publicly available samples from the Gene Expression Omnibus. Our data show a vast number of genes and networks expressed predominantly in HObs compared with closely related cells such as fibroblasts or chondrocytes. For instance, genes in the insulin-like growth factor (IGF) signaling pathway were enriched in HObs (P = 0.003) and included the binding proteins (IGFBP-1, -2, -5) and IGF-II and its receptor. Another HOb-specific expression pattern included leptin and its receptor (P < 10−8). Furthermore, after stimulation of HObs with BMP-2 or dexamethasone, the expression of several interesting genes and pathways was observed. For instance, our data support the role of peripheral leptin signaling in bone cell function. In conclusion, we provide the landscape of tissue-specific and dynamic gene expression in HObs. This resource will allow utilization of osteoblasts as a model to study specific gene networks and gene families related to human bone physiology and diseases.
- gene expression
- bone morphogenic protein-2
osteoporosis is a skeletal disorder characterized by compromised bone strength and increased risk of fracture in which an imbalance in the regulation of bone remodeling, e.g., bone resorption and bone formation, occurs and the amount of resorption exceeds the formation (1). The disease can be idiopathic or a consequence of identifiable predisposing factors, such as gonadal hormone deficiency or glucocorticoid (GC) treatment. Osteoporosis can be treated by antiresorptive and anabolic agents. The antiresorptive agents (such as bisphosphonates) are known to decrease bone resorption but only modestly increase bone mineral density (BMD) by 5–10% (36, 45). In contrast, the anabolic agent parathyroid hormone 1-34 can increase BMD by up to 50% by a process that results in improved bone formation (38). Consequently, the search for novel therapeutic approaches focuses on discovering anabolic agents that will increase bone formation.
A key player in the bone remodeling process is the osteoblast, which, in addition to initiating the resorption process, is responsible for the production of the bone matrix (10). The osteoblasts are thus the focal point for development and assessment of novel therapeutic agents for osteoporosis. Studying mature human osteoblasts in vitro most often involves collection of bone fragments from orthopedic surgery and subsequently cultivation of human trabecular bone cells (35, 42, 46). The relative ease of human osteoblast-like cell (HOb) isolation from an abundant clinical source, which would otherwise be disposed, makes them an accessible primary cell model for studying cell physiology in human cell panels.
The osteoblasts are of mesenchymal origin and are in cell culture morphologically similar to fibroblasts. The only known osteoblast-specific feature is the ability to produce mineralized extracellular matrix. There is no evidence to date that bone matrix mineralization is orchestrated by genes selectively expressed in osteoblasts. Only two osteoblast-specific transcripts have been identified: one encoding Cbfa1 (11) and the other encoding osteocalcin (9). In addition, all fibroblast-expressed genes are suggested to be expressed in osteoblasts as well.
Two factors with important but different roles in bone remodeling are GCs and bone morphogenic proteins (BMPs). GCs are known to induce osteoporosis, most likely via a process involving accelerated osteoblast apoptosis. In contrast, other studies have shown that GCs stimulate bone formation in vitro, indicating the complex role of GCs in bone metabolism. In addition, a recent paper reported that GCs mediate bone-suppressive effects indirectly via the GC receptor on osteoclasts (29). BMPs, and especially BMP-2, are potent stimulators of bone formation (6, 57). Moreover, BMP-2 is a necessary component of the signaling cascade that governs fracture repair (53).
In the past years, extensive research has been devoted to studying gene expression involved in bone metabolism with cultured cell lines of humans and mice (13). It has been shown that serial passage has an effect on the gene expression in MC3T3-E1 preosteoblastic cells (25), which indicates the risk of introducing significant bias when using immortalized and passaged cell lines. The aim of the present study was to characterize cultured primary HObs and identify signaling cascades and gene regulatory networks involved in bone remodeling by using genomewide gene expression profiling in resting as well as BMP-2- and dexamethasone-induced bone cells. We obtained HObs from bone fragments collected from the proximal femoral shaft of patients undergoing total hip replacement. We measured expression profiles in HObs in a resting state and after stimulation with BMP-2 and dexamethasone. We also compared the osteoblast transcriptome to human primary fibroblasts and other commonly utilized human cell systems used in molecular biology and human genomics. Our results highlight previously unappreciated features of the osteoblast transcriptome.
MATERIALS AND METHODS
Bone cell culture for whole genome expression profiling.
Human trabecular bone from the proximal femoral shaft was collected from two male donors: a 55-year-old and a 25-year-old man, both undergoing total hip replacement. The bone chips obtained from each donor were thoroughly minced and washed with PBS (National Veterinary Institute of Sweden, Uppsala, Sweden). The minced bone chips were divided into three portions, and each portion was cultured separately, i.e., three biological replicates from each donor. The cells were grown in medium containing α-MEM (Sigma-Aldrich, Haverhill, UK) supplemented with 2 mmol/l l-glutamine, 100 U/ml penicillin, 100 mg/ml streptomycin (National Veterinary Institute of Sweden), and 10% fetal bovine serum (Sigma-Aldrich) at 37°C with 5% CO2 until confluence was reached. The culture medium was changed twice weekly. At 70–80% confluence, each cell culture was trypsinized and subcultured in six-well plates (100 000 cells/well) for 12 days. Before treatment, the cells were starved for 20 h by addition of complete cell medium containing 0.5% fetal bovine serum. The cells were then incubated for 2 h and 24 h with 10−7 M dexamethasone and 10−4 mg/ml of recombinant human (rh)BMP-2 with the same concentration of vehicle (control), respectively.
At the two time points, the cell medium was removed and the cells were harvested by adding 500 μl of lysis buffer (Qiagen). The cell lysates were homogenized with QIAshredder (Qiagen, Germany) homogenizers and stored at −70°C until RNA extraction.
The study was approved by the local ethics committees (Dnr Ups 03-561\ and McGill University Institutional Review Board).
RNA was extracted from the cell lysates with the RNeasy Mini Kit (Qiagen). High RNA quality was confirmed for all samples with an Agilent 2100 BioAnalyzer (Agilent Technologies, Palo Alto, CA), and the concentrations were determined with Nanodrop ND-1000 (NanoDrop Technologies, Wilmington, DE).
For the comparison study, expression profiling of control samples at t = 24 h were performed for each biological replicate with the Affymetrix Human Genome U133 plus 2.0 Array (Affymetrix, Santa Clara, CA), i.e., the three biological replicates from the two donors were analyzed separately. For the treatment study, expression profiling of all samples at both time points from one individual were performed for each biological replicate.
One microgram of RNA was reverse transcribed into cDNA, and in vitro transcription was performed to generate biotin-labeled cRNA for subsequent hybridization. Hybridized target cRNA was then stained with streptavidin-phycoerythrin, and the arrays were scanned with a GeneArray Scanner at an excitation wavelength of 488 nm.
Data extraction from Gene Expression Omnibus.
Microarray expression data of 75 publicly available samples from 25 data series that were hybridized to Affymetrix Human Genome U133 Plus 2.0 Arrays (Affymetrix) were downloaded from GEO (14) (Supplemental Table S1).1 The samples comprised gene expression data from a wide variety of different homogeneous tissues and primary cells (e.g., adipose, liver, cartilage, skeletal muscle, cerebellum, smooth muscle cells, fibroblasts, mesenchymal stem cells, keratinocytes, embryonic stem cells, dendritic cells, monocytes, neutrophils, lymphocytes, and macrophages) and were selected on the basis of the following criteria: 1) raw data (CEL) files were available for downloading; 2) at least three biological replicates were analyzed and available; and 3) primary tissues/cells were collected from healthy donors. No cancer cell lines were included in this initial comparison. In addition, we included expression data from frequently used immortalized cell lines, e.g., Jurkat (T cell lymphocytic cell line), Caco-2 (colonic adenocarcinoma cell line), LCL (lymphoblastoid cell line), MCF-10F (human breast epithelial cell line), LNCaP (prostate cancer cell line), H292 (human pulmonary mucoepidermoid carcinoma cell line), THP-1 (human acute monocytic leukemia cell line), and NB4 (human promyelocytic leukemia cell line), to compare with the human osteoblast transcriptome (Supplemental Table S2).
Robust multichip average adjustments.
All expression values from the GEO samples as well as the human primary bone samples were generated by the robust multichip average (RMA) algorithm (3, 26) in Bioconductor with the R 2.5.0 package (18). The RMA normalization was applied to each GEO data set individually and across all data sets as described in results.
Cluster and statistical analysis.
To cluster the samples on the basis of their similarities in gene expression we initially used unsupervised clustering methods, e.g., principal component analysis (PCA) (56) and hierarchical clustering (HC). PCA condenses the variables in the original data into a new set of variables (principal components) to summarize the features of the data. These new variables together account for as much of the variance in the original variables as possible while remaining mutually uncorrelated and orthogonal. Unsupervised HC has been described for gene expression and has clearly proven valuable (15, 49) in identifying sets of correlated genes and/or samples with similar behavior across the experiments. We used MeV (MultiExperiment Viewer) version 4 software (48) and average linkage for determining cluster-to cluster distances. PCA and HC are frequently used for clustering samples that show similar expression pattern but have limitations in their ability to capture the full structure of the data. Since the HOb samples included in the present study were compared with samples deposited in the GEO, biological and technical differences may have had an impact on the array comparison. Nonnegative matrix factorization (NMF) (31) was recently introduced as a new clustering method to extract relevant biological correlations from gene expression data (4). The method is designed to capture alternative structures inherent in the data and provides biological insight by organizing both sample and gene information. NMF analysis was performed in MATLAB with codes for NMF divergence-reducing equations, as well as for model selection and reordering of consensus matrixes (http://www.broad.mit.edu/cancer; Ref. 4). The numbers of classes to be determined were k = 2–5.
The Cyber-t test was used to determine the significance between the observed differences in gene expression. The Cyber-t test is based on simple t-tests and uses the observed variance of replicate gene measurements across replicate experiments, thereby accommodating noise, variability, and low replication often typical of microarray data (2). Since the number of tests (genes) is large, the P values are adjusted for multiple comparisons with the false discovery rate (FDR) algorithm (50) in Bioconductor with the R 2.5.0 package (18).
Gene network and pathway analysis.
To visualize whole genome expression data in the context of biological networks, the gene expression data of treated and untreated samples were analyzed with the Ingenuity Pathway Analysis (IPA) system (Ingenuity Systems, Mountain View, CA; www.ingenuity.com). Data sets containing gene identifiers and corresponding expression values were uploaded in the application. Each gene identifier was mapped to its corresponding gene object in Ingenuity Pathways Knowledge. A fold change cutoff was set to identify genes whose expression was significantly differentially up- or downregulated as a result of the treatment. These genes, called Focus Genes, were overlaid onto a global molecular network developed from information contained in the Ingenuity Pathways Knowledge Base. Networks of these Focus Genes were then algorithmically generated on the basis of their connectivity.
The functional analysis identified the biological functions that were most significant to the data set. Genes from the data set that met the cutoff and were associated with biological functions in Ingenuity Pathways Knowledge were considered for the analysis. Fisher's exact test was used to calculate a P value determining the probability that each biological function assigned to the data set is due to chance only.
Canonical pathways analysis identified the pathways from the IPA library of canonical pathways that were most significant to the data set. Genes from the data set that met the cutoff in the functional analyses and were associated with a canonical pathway in the Ingenuity Pathways Knowledge Base were considered for the analysis. The significance of the association between the data set and the canonical pathway was measured in two ways. 1) A ratio of the number of genes from the data set that map to the pathway divided by the total number of genes that map to the canonical pathway was displayed. 2) Fischer's exact test was used to calculate a P value determining the probability that the association between the genes in the data set and the canonical pathway is explained by chance alone.
cDNA synthesis and real-time PCR.
The effect of dexamethasone and BMP-2 on gene expression in HObs was verified by quantitative real-time (RT) PCR. The same sources of total RNA used in the microarray experiments were used for the data validation by RT-PCR. RNA from the three biological replicates were each annealed to 250 ng of random hexamers (Invitrogen, Carlsbad, CA) at 70°C for 10 min. First-strand cDNA synthesis was performed with SuperScript II reverse transcriptase (Invitrogen) according to manufacturer's instructions. Each pair of primers was tested by agarose gel electrophoresis to verify the presence of a single band of the predicted size. Primers used are presented in Supplemental Table S3, and all primers were designed with Primer3 software (http://frodo.wi.mit.edu/).
Each gene in the quantitative RT-PCR assay was analyzed in duplicate as well as a calibration curve from a twofold dilution series of nonstimulated trabecular bone cell cDNA and nontemplate control samples.
The RT-PCR assays were performed on the Rotor-Gene 6000 real-time rotary analyzer (Corbett Life Sciences, Sydney, Australia) with Platinum SYBR Green qPCR SuperMix-UDG (Invitrogen) according to the manufacturer's recommendation. The cycling conditions on the Rotor-Gene 6000 real-time rotary analyzer were 4 min at 95°C and 40 cycles of 20 s at 95°C, 30 s at 58°C, and 30 s at 72°C, followed by the dissociation protocol at 72°C.
Results of the experimental samples were analyzed with the relative standard curve method. The mean and SD values of each technical replicate sample were calculated, and the mean value was then normalized to the 18S mean value to obtain the target-to-reference ratio. This ratio was then averaged from the biological replicate (n = 3) samples. Finally, the fold differences between control sample normalized target and dexamethasone or BMP-2 sample normalized target were calculated.
Cell culture for protein analysis.
Human trabecular bone from the proximal femoral shaft was collected from a female patient undergoing total hip replacement. The bone chips were minced thoroughly and washed with PBS (National Veterinary Institute of Sweden), cultured in cell medium containing α-MEM (Sigma-Aldrich) supplemented with 2 mmol/l l-glutamine, 100 U/ml penicillin, 100 mg/ml streptomycin (National Veterinary Institute of Sweden), and 10% fetal bovine serum (Sigma-Aldrich), and grown at 37°C with 5% CO2 until confluence was reached. The culture medium was changed twice every week. At 70–80% confluence, the cells were trypsinized and subcultured in a 24-well plate (100,000 cells/well). Before treatment, the cells were starved for 20 h by addition of complete cell medium containing 0.5% fetal bovine serum. The cells were then incubated for 24 h and 48 h with 10−7 M of dexamethasone and with the same concentration of vehicle, respectively. Each treatment was performed in three replicates. At the two time points, the cell medium was removed and stored at −20°C until protein analysis.
Assays of protein production.
The concentration of IL-6 and leptin in the cell culture media of HObs treated with dexamethasone and rhBMP-2 was measured with the Human IL-6 Immunoassay (QuantiGlo, R&D Systems, Minneapolis, MN) and the Human Leptin Immunoassay (Quantikine, R&D Systems) according to the manufacturer's instructions.
Gene expression profiling of human tissues and cells.
Whole genome expression profiling of untreated HObs from two unrelated individuals, each with three biological replicates, was conducted with Affymetrix GeneChip U133 +2 microarrays. In addition, 25 data sets in triplicate generated on the same GeneChip were downloaded from the GEO, and the normalized expression values from the total data set of 81 samples were used in the clustering analyses. In the GEO data set selection we focused on cell types with independently (in different labs) generated data to allow a more robust inference of differences between cell types. In a further attempt to detect true biological differences between cell types as opposed to methodological differences (in culture, sample preparation, or chip hybridization) we used several normalization (data set-by-data set or global normalization, respectively) and clustering (PCA, HC, NMF) approaches.
In the PCA, a total of 23 components explaining >0.1% of the overall variability were detected. Determining the true dimensionality of the data and eliminating noisy components is often ad hoc, and many exclusion criteria exist. Eliminating low-variance components, while reducing noise, may also result in discarding information. We used a criterion that discarded all components that account for <3% of the overall variability. This resulted in three remaining components that together accounted for 87% of the overall variability. Four major groups of samples were identified by the PCA analyses and included cells from the hematopoietic cell lineage (monocytes, macrophages, lymphocytes, neutrophils, dendritic cells); samples derived from the central nervous system (CNS; 2 independent sources of cerebellum); cells of mesenchymal origin (fibroblasts, cartilage, mesenchymal stem cells, trabecular bone-derived cells); bronchial smooth muscle cells; keratinocytes and embryonic stem cells; and finally a cluster containing adipose, skeletal muscle, and liver samples (Fig. 1A).
The HC algorithm confirmed the clustered groups of samples with the HObs closely related to fibroblast, mesenchymal stem cells, and fetal cartilage (Fig. 1B). With average linkage and Euclidean distance measure, four clusters were identified with a linkage threshold of 0.62 (distance = 1.0); adipose tissue and skeletal muscles were closely related to the cluster including the mesenchymal stem cells, embryonic stem cells, bronchial smooth muscle cells, keratinocytes, cartilage, fibroblasts, and HObs.
With NMF the same robust clustering pattern demonstrated in the PCA and HC was seen (Fig. 1C). With rank k = 2 NMF separated the hematopoietic cells from the rest of the samples, and a higher k revealed further partitioning of the samples. For k = 3, the partition reflects the HOb, fibroblast, cartilage, embryonic stem cell, mesenchymal stem cell, bronchial smooth muscle cell, and keratinocyte distinction. For k = 4, a fourth class appeared with embryonic stem cells. Finally, for k = 5, CNS-derived samples were subclustered in the fifth class. The NMF approach for data set-by-data set normalization yielded nearly identical clusters, suggesting that our analysis across all studies was not significantly influenced by the RMA normalization (Supplemental Fig. S1, Supplemental Table S4).
Common cell models compared with trabecular bone-derived osteoblasts.
In a separate set of analyses we compared the HOb transcriptome to the transcriptome of common immortalized cell models such as LCLs and MCF-10F. These comparisons allowed the evaluation of gene networks observed in experiments involving these cell systems in relation to bone cells. Our results showed three distinct clusters with PCA (data not shown), which was confirmed by HC analysis (Supplemental Fig. S2). With an average linkage threshold of 0.50 (distance = 1.3) the different cancer cell lines were clustered in two separate groups: cells of hematopoietic origin (Jurkat, LCLs, THP-1, and NB4 cells) and cells of epithelial origin (Caco-2, H292, MCF-10F, LNCaP). The HObs were grouped in a third cluster.
Pathway/network-specific genes in cultured trabecular bone cells.
Both PCA and unsupervised HC revealed a close link between normal human fibroblasts and HObs with an average linkage threshold of 0.76 (distance = 0.75, Fig. 1, A and B), which is high compared with most other cell types (e.g., adipocytes and skeletal muscle show linkage threshold of 0.62, leukocytes and lymphocytes 0.60, and neuronal tissue an average of 0.57). However, when comparing genes expressed in cultured HObs with those expressed in fibroblasts, we detected a total of 571 probe sets with greater than fourfold differences in gene expression. This corresponds to 486 unique genes that were expressed at a higher level in HObs compared with fibroblasts. From these genes we listed the top highest differentially expressed genes in HObs that are known to be involved in bone metabolism (Table 1). A full gene list of differentially expressed genes in HObs compared with fibroblasts is shown in Supplemental Table S5. These highly expressed genes were associated with formation of bone (P = 0.006), volume of bone (P = 0.01), osteoclastogenesis of bone (P = 0.008), and quantity of trabecular bone (P = 0.02) based on literature mining by the IPA software (see materials and methods). Moreover, the insulin-like growth factor (IGF)-I signaling pathway was one of the top three canonical signaling pathways that was specifically linked to the HOb transcriptome and included the IGF binding protein (IGFBP) genes as well as the protein kinases PDPK1, PIK3CB, and PRKCH (ratio of number of genes identified per number of genes in the pathway = 8/92, P value 0.003). Genes that were expressed lower in human osteoblasts than in fibroblasts included MMP3, GREM2, MMP1, and SFRP2.
Fetal cartilage showed an average linkage to the bone cells of 0.68 (distance = 1.0), and comparative analyses by fold change and IPA software identified key differences that were related to leptin (LEP)-leptin receptor (LEPR) signaling (bone specific) and IGF-I signaling (bone specific). The most relevant functions associated with cartilage tissue were development and function of connective tissue including condensation (P = 0.005) and development (P = 0.01) of cartilage tissue and development of endochondral bone (P = 0.02).
Other pairwise comparisons, including those of HObs to more distantly related cell types and selected immortalized cells, are shown as supplemental data (Supplemental Results). The top canonical pathways expressed in LCLs are presented in Supplemental Table S6. A full gene list of differentially expressed genes in HObs compared with LCLs is shown in Supplemental Table S7.
Gene expression profiling of human bone cells after treatment with BMP-2.
Expression profiling of three biological replicates of HObs exposed to BMP-2 was conducted to identify immediate-early (2-h exposure) and late (24-h exposure) genes.
Overall, at the two time points significantly more genes were upregulated in response to BMP-2 than were downregulated. At the 2 h time point, a total of 61 probe sets exhibited fold changes above the twofold change threshold with RMA intensities. This corresponds to 47 unique genes that were upregulated and no genes that were downregulated after BMP-2 stimulation (Table 2). All genes were included in the network analysis as well as in the function and pathway analyses using the IPA software. The top molecular and cellular functions significantly associated with the regulated genes after 2 h of BMP-2 exposure by IPA are presented in Table 3.
At the 24 h time point, a total of 427 probe sets showed a change larger than the twofold threshold with the RMA intensities, corresponding to 306 unique genes.
Genes that were altered during 2 h and 24 h of treatment with BMP-2 are presented in Table 4. Full gene lists at 2 h and 24 h of BMP-2 treatment are shown in Supplemental Tables S8 and S9.
Validation of BMP-2-responsive genes by RT-PCR.
A total of seven genes (CXCL2, DLX2, FGF18, HES1, ID2, SMAD6, and SMAD7), upregulated after 2 h of BMP-2 exposure, were validated by RT-PCR (Fig. 2). These genes were selected on the basis of significant results from gene expression profiling as well as their potential effect on bone formation. In both cases, the direction of change was similar between microarray and RT-PCR. Moreover, microarray and RT-PCR fold change findings were highly correlated (r2 = 0.89), suggesting that microarray findings for other genes in the study should be qualitatively correct, even though absolute change levels are different.
Gene expression profiling of human bone cells after treatment with dexamethasone.
Whole genome expression profiling of three biological replicates of HObs was conducted on one individual exposed to dexamethasone for 2 h and 24 h. At t = 2 h, a total of 60 genes exhibited fold changes above the twofold threshold with RMA intensities. This is a response quantitatively similar to that observed for BMP-2 at the same data point. In contrast, striking alterations in gene expression were observed at the 24 h data point of dexamethasone treatment. In total, 346 genes were upregulated and 369 genes were downregulated. The altered top genes during 24 h dexamethasone treatment are presented in Table 5. Of the 715 genes, 520 genes were eligible for function and pathway analyses with IPA whereas 554 genes were eligible for network analysis.
A well-known effect of GCs on bone-forming cells is accelerated apoptosis, which was the top molecular and cellular function significantly associated with the 24 h expression data set. Of the 520 genes that were included in the function analysis, 183 genes were classified as apoptosis associated (P value <10−15–10−3). Other functional classes that were significantly represented in the dexamethasone-induced transcriptome included cellular growth and proliferation (200 genes, P value <10−11–10−3) and cell morphology (106 genes, P value <10−8–10−3).
The IGF-I signaling pathway was one of the top canonical signaling pathways that was downregulated after 24 h of dexamethasone exposure (ratio 14/90, P value 0.0003). The pathway analysis indicated that the phosphatidylinositol 3-kinase (PI3K) signaling and the ERK signaling downstream of IGF-I action might also be downregulated (Fig. 3).
However, the most significant canonical pathway to the dexamethasone-induced data set, with IPA, was the circadian rhythm signaling (ratio 5/31, P value 0.004; Table 6).
Full gene lists at 2 h and 24 h of dexamethasone treatment are shown in Supplemental Tables S10 and S11.
Validation of dexamethasone-responsive genes by RT-PCR and protein analyses.
The expression of seven genes (IL-6, LEP, LEPR, sFRP1, SOX9, SPP1, and TNFAIP6) that were up- or downregulated after 24 h of dexamethasone exposure was validated by RT-PCR (Fig. 4A). Two of those seven genes (IL6 and LEP) were further validated by protein analyses (Fig. 4, B and C). These two genes were selected based on significant results from gene expression profiling and because of their potential effect in bone remodeling.
The effects of dexamethasone on IL-6 and leptin protein levels in human trabecular bone cells were found to increase with time, and the maximum effect was seen after 48 h of stimulation (Fig. 4, B and C). No difference in IL-6 protein levels was seen after 24 h. However, after 48 h of dexamethasone stimulation, IL-6 protein production was significantly decreased (P value <0.01, Student's t-test). After 24 h and 48 h, the leptin levels in control samples were below the detection level of the immunoassay. In all cases, the direction of change was similar between microarray, RT-PCR, and protein analyses. Moreover, microarray and RT-PCR fold change findings were highly correlated (r2 = 0.92).
Our detailed and comprehensive comparative transcriptome study demonstrates the relatively close relationship of HObs to human fibroblasts and to fetal cartilage. These results, as well as the comparison to other human cells, provide a point of reference for assessing the relevance of functional data observed in human osteoblasts. Numerous other variables (such as cell culture conditions and sample handling or variability in GeneChip protocols) could account for differences observed between our data and publicly available expression data. Arguments favoring true biological differences between transcriptomes in different cell types detected in our analysis are that 1) similar cell and tissue types from independent GEO data sets showed coclustering, 2) the clustering yielded similar results with three different methods, and 3) different data normalization strategies gave similar results.
Despite the close relation between the expression profiles of osteoblasts compared with that of fibroblasts, we found interesting differences. The osteoblast-enriched transcriptome in relation to fibroblasts and chondrocytes includes upregulated genes known to be involved in bone metabolism (collagen genes, osteocalcin, ANKH, ALPL, etc.). The top upregulated pathways are the IGF signaling and the LEP-LEPR pathway, respectively. Furthermore, the transcriptome of osteoblasts is clearly distinct compared with cells of lymphoid origin. The latter cells are commonly used in biobanking and increasingly used as a general model system in human genomics (reviewed in Refs. 7, 44). Given these differences as well as the accessibility of abundant clinical material, the simplicity of osteoblast isolation and culture makes them a lucrative primary cell model, and the use of osteoblasts will complement studies carried out on more commonly employed models.
The IGF system plays an important role in the autocrine and paracrine regulation of bone formation and remodeling and are the most abundant growth factors produced by osteoblasts. The most pronounced upregulated genes in the IGF signaling pathways include the binding proteins (IGFBP-1, -2, and -5) and IGF-II and its receptor. This is in agreement with previous data on bone tissues and osteoblastic cells (5, 8, 17, 23, 24, 34) that IGF-II and binding proteins IGFBP-2 and IGFBP-5 are most represented in bone cells and tissues (43). IGF-II and IGFBP-2 are known to be closely interrelated. IGF-I stimulates IGFBP-2 synthesis, and IGFBP-2 shows greater affinity for IGF-II than for IGF-I (43). IGFBP-5 stimulates markers of bone formation both in vivo and in vitro, probably by binding directly to sites that are independent of the IGF receptor (37). In addition, IGFBP-5 is abundant in serum and bone during normal skeletal development, but levels decrease in osteoporosis (27).
A striking HOb-specific expression pattern includes leptin (LEP) and its receptor (LEPR), a single-transmembrane-domain receptor of the cytokine receptor family (52). Leptin exerts its effect by binding to the LEPR and activates the JAK/STAT pathway (19, 54). Recently it was shown that this process appears to be mediated by the PI3K/Akt pathway through activated JAK1 (41), a gene that was highly expressed in HObs (data not shown). Also, genes in the JAK/STAT and the PI3K/Akt pathways were highly expressed in the HObs (data not shown), indicating that osteoblasts possess the receptor apparatus necessary to convey the leptin signaling. It is also known that leptin enhances the expression of transforming growth factor (TGF)-β (12) and alkaline phosphates (ALP) (55), which are all genes that were upregulated in the HObs compared with fibroblasts. The LEPR was expressed in most of the samples included in our transcriptome comparison. On the other hand, LEP was only expressed in adipose tissue and in HObs. LEP has been advocated as a centrally acting factor responsible for inhibiting the accumulation of bone mass (reviewed in Ref. 51). However, recent investigations and data from our study unequivocally indicate that leptin is expressed in HObs and might act as a local (autocrine) factor in bone remodeling. It is known that exogenously added leptin causes osteoblastic cell proliferation and differentiation, while also rendering osteoblasts more efficacious in terms of mineralization (20).
On stimulation of the HObs with BMP-2, a potent stimulator of bone formation, a large number of negative regulators were upregulated. The genes included the BMP inhibitor Noggin, the inhibitory Smads (SMAD6 and SMAD7), the Notch pathways genes (HEY1 and HES1), and other genes inhibiting proliferation pathways (DUSP2, TNFAIP3, GADD45b, and NFKBIA). In general, the BMP signaling is complex, and there are multiple potential points of regulation. Expression studies as well as gain- and loss-of-function experiments suggest that SMAD6 participates in an important negative feedback loop in which BMP-2-mediated effects on differentiation are reduced by induction of SMAD6 (32). Smad signaling has negative feedback regulations at the cytoplasmic or nuclear level, which are important to restrict or terminate the biological effect of BMPs. There is also evidence that a negative feedback loop involving Noggin might account for the precise localization of BMP signaling (28).
In addition, a number of BMP-2-responsive genes, known to promote bone formation, were identified. The most upregulated genes, after both 2 h and 24 h of BMP-2 exposure, were the inhibitory helix-loop-helix transcription factors Id1 and Id3 and the distal-less homeobox 2 and 5 transcription factors (DLX2 and DLX5). All four genes are known to code for important transcription factors promoting bone formation (21, 33). Of note, the transcription factor Cbfβ, which is known to be expressed in developing bone and forms a functional interaction with Runx2, was significantly upregulated after 2 h of BMP-2 stimulation. Similarly, the gene Heme oxygenase 1 HMOX1 (HO-1), which has been shown to suppress osteoclastogenesis (58), was significantly upregulated after BMP-2 stimulation.
GCs are known to induce osteoporosis, most likely via a process involving accelerated osteoblast apoptosis (21, 29). This was the top molecular and cellular function significantly associated with the dexamethasone expression data set. Of the 520 genes included in the function analysis, 183 genes are classified as apoptosis associated.
The IGF-I signaling pathway is one of the top canonical signaling pathways that were downregulated after 24 h of dexamethasone exposure. The dexamethasone-regulated genes in the IGF-I signaling pathway included Cyr61, c-fos, Foxo1a, IGF-1, IGF1R, IGFBP-1, IGFBP-5, IRS2, MAPk1, PIK3R3, PRKAG2, PRKAR2B, RASA1, and RRAS. The pathway analysis indicated that the PI3K signaling and the ERK signaling downstream of IGF-I action were downregulated. A more detailed analysis of the top downregulated genes revealed a large number of genes known to be involved in bone remodeling. The HMOX1 (HO-1) gene that was significantly upregulated after BMP-2 stimulation was highly downregulated after dexamethasone stimulation. The Sequestosome 1 gene, SQSTM1, recently associated with predisposition to Paget disease (30), is the most exaggerated example of abnormal bone remodeling. We found SQSTM1 significantly downregulated after dexamethasone stimulation but significantly upregulated after 24 h of BMP-2 exposure. Two other downregulated genes in dexamethasone stimulation, known to be involved in the bone remodeling, were the bone sialoprotein gene, the osteopontin gene, and the matrix metalloproteinase 1 (MMP1) gene. Moreover, genes in the Wnt/β-catenin pathway were shown to be regulated by dexamethasone. The Wnt antagonists frizzled-related protein (FRZB), secreted frizzled-related protein 1 (SFRP1), and Dickkopf (DKK1) were all upregulated as well as the Wnt receptor Frizzled 5. In addition, the Wnt genes Wnt5b, -2B, and -3 were all downregulated by dexamethasone. Finally, leptin and its receptor were highly upregulated after dexamethasone stimulation, which further indicates a role of peripheral leptin signaling in bone cell function.
Additional findings related to dexamethasone exposure included downregulation of the circadian rhythm signaling. In bone, clock genes are involved in the circadian oscillation of bone formation and extracellular matrix expression. The expression of type I collagen, type II collagen, and Runx2 have been shown to have rhythmic expression in bone (16, 22, 39, 40, 47). In the present study, a total of 12 clock genes were shown to be significantly downregulated after dexamethasone stimulation, suggesting the negative role of GCs in circadian rhythm signaling in human bone cells.
In conclusion, our comprehensive comparative study of transcriptomes demonstrates a vast number of genes and gene networks that are predominantly expressed in HObs compared with closely related cells of mesenchymal origin such as fibroblasts. Furthermore, induced gene expression by factors with important but different roles in bone remodeling showed upregulation of known and novel genes and pathways. For example, our data support the role of peripheral leptin signaling in bone cell function. Together, our data show that osteoblasts are an excellent model and a useful model to study specific gene networks and gene families related to human bone physiology and disease.
The authors thank research technicians and nurses at the Departments of Medical and Surgical Sciences at Uppsala University for skillful work and handling of the samples. The authors also thank Haig Djambazian at the Genome Quebec Innovation Centre for valuable help with statistical analyses. The work was supported by the Swedish Research Council, Genome Canada, and Genome Quebec.
↵1 The online version of this article contains supplemental material.
Address for reprint requests and other correspondence: T. Pastinen, McGill Univ. and Genome Quebec Innovation Centre, 740 Dr Penfield Ave., Rm 6202, Montreal, QC, H3A1A4, Canada (e-mail:).
Article published online before print. See web site for date of publication (http://physiolgenomics.physiology.org).
- Copyright © 2008 the American Physiological Society