|Home | About | Journals | Submit | Contact Us | Français|
Rosiglitazone (Rosi), a member of the thiazolidinedione class of drugs used to treat type 2 diabetes, activates the adipocyte-specific transcription factor peroxisome proliferator-activated receptor gamma (PPARγ). This activation causes bone loss in animals and humans, at least in part due to suppression of osteoblast differentiation from marrow mesenchymal stem cells (MSC). In order to identify mechanisms by which PPARγ2 suppresses osteoblastogenesis and promotes adipogenesis in MSC, we have analyzed the PPARγ2 transcriptome in response to Rosi. A total of 4,252 transcriptional changes resulted when Rosi (1 μM) was applied to the U-33 marrow stromal cell line stably transfected with PPARγ2 (U-33/γ2) as compared to non-induced U-33/γ2 cells. Differences between U-33/γ2 and U-33 cells stably transfected with empty vector (U-33/c) comprised 7,928 transcriptional changes, independent of Rosi. Cell type-, time- and treatment-specific gene clustering uncovered distinct patterns of PPARγ2 transcriptional control of MSC lineage commitment. The earliest changes accompanying Rosi activation of PPARγ2 included effects on Wnt, TGFβ/BMP and G-protein signaling activities, as well as sustained induction of adipocyte-specific gene expression and lipid metabolism. While suppression of osteoblast phenotype is initiated by a diminished expression of osteoblast-specific signaling pathways, induction of the adipocyte phenotype is initiated by adipocyte-specific transcriptional regulators. This indicates that distinct mechanisms govern the repression of osteogenesis and the stimulation of adipogenesis. The co-expression patterns found here indicate that PPARγ2 has a dominant role in controlling osteoblast differentiation and suggests numerous gene-gene interactions that could lead to the identification of a “master” regulatory scheme directing this process.
Peroxisome proliferator-activated receptor gamma (PPARγ) is a member of the nuclear receptor family of transcription factors, a large and diverse group of proteins that mediate ligand-dependent transcriptional activation and repression. PPARγ plays an important role in the control of adipocyte development as well as glucose and fatty acid metabolism [Rosen and Spiegelman, 2001]. Alterations in these pathways can lead to the development of obesity, insulin resistance, and Type 2 diabetes mellitus (T2D) [Sharma and Staels, 2007]. Activation of the PPARγ protein is essential for vital processes and affects cell proliferation, differentiation and neoplastic transformation [Lehrke and Lazar, 2005]. The discovery of significant cross-talk between PPARγ and other nuclear receptors, such as retinoid, estrogen, and the vitamin D receptors, strongly supports a pleiotropic role for this protein in the control of cell homeostasis [Gimble et al., 2006].
An essential role of PPARγ in the maintenance of bone homeostasis is implicated in animals and humans. In animal models, a decreased PPARγ activity leads to increased bone mass and osteoblast number [Akune et al., 2004; Cock et al., 2004], whereas increase in PPARγ activity due to treatment with the anti-diabetic thiazolidinedione (TZD) drug rosiglitazone (Rosi) results in significant decreases in BMD, bone volume and changes in bone microarchitecture [Ali et al., 2005; Rzonca et al., 2004; Soroceanu et al., 2004; Sottile et al., 2004]. This bone loss is associated with a decreased number of osteoblasts and an increased number of adipocytes [Ali et al., 2005; Rzonca et al., 2004].
In humans, an administration of TZDs results in progressive bone loss and diminished levels of circulating bone formation markers in older women [Grey et al., 2007; Schwartz et al., 2006]. The results of ADOPT studies, which evaluated the effectiveness of known anti-diabetic drugs on the maintenance of blood glucose levels in prediabetic patients, suggested that Rosi increased the frequency of fractures in women [Kahn et al., 2006; Kahn et al., 2008]. Clinical studies of PPARγ gene polymorphisms in different human populations indicate a role for this transcription factor in the regulation of bone mass and predisposition to the bone loss in a high fat diet conditions [Ackert-Bicknell et al., 2008; Rhee et al., 2005].
The mechanism of TZDs-induced bone loss involves their effects on differentiation of bone cells. In cells of mesenchymal lineage, the TZD-activated PPARγ2 isoform induces adipocyte and suppresses osteoblast differentiation [Lecka-Czernik et al., 1999; Lecka-Czernik et al., 2002]. In cells of hematopoietic lineage, the TZD-activated PPARγ1 isoform positively regulates differentiation of osteoclasts [Wan et al., 2007]. Here, we analyzed the mechanistic consequences of Rosi activation of PPARγ2 on the global transcriptional response in a cellular model of marrow mesenchymal stem cell (MSC) differentiation. Using an ANOVA-based approach we describe time-dependent relationships in MSC-based gene expression and key aspects of cell conversion from osteoblast progenitor cells to mature adipocytes, including connections between genes and gene ontology groupings.
Our cellular model of PPARγ2 control of marrow MSC differentiation was originally developed to study the mechanisms by which PPARγ2 suppresses osteogenesis and promotes adipogenesis. However, we found that activation of this nuclear receptor by the exogenous ligand Rosi modulated gene expression for multiple pathways essential for regulation of cell homeostasis. Moreover, our results suggest that PPARγ2 has a profound effect on gene expression, even in the absence of exogenous ligand. These effects suggest a role in the modulation of MSC phenotype in conditions of increased PPARγ2 expression, including aging.
Murine marrow-derived U-33 cells (previously referred as UAMS-33) represent a clonal cell line spontaneously immortalized in long term bone marrow cultures [Lecka-Czernik et al., 1999]. To study the effect of PPARγ2 on marrow mesenchymal stem cell differentiation, U-33 cells were stably transfected with either a PPARγ2 expression construct (U-33/γ2 cells) or an empty vector control (U-33/c cells), as described previously [Lecka-Czernik et al., 1999]. Several independent clones were retrieved after transfection and carefully analyzed for their phenotype. Clone 28.6 (representing U-33/γ2 cells) and clone γc2 (representing U-33/c cells) were used here. Experimental design, cell maintenance, cell treatment, and harvesting of RNA samples were described previously [Lecka-Czernik et al., 2007]. Primary bone marrow cultures were established from the bone marrow isolated from femora of 6 mo old C57BL/6 mice, which were obtained from the colony maintained by the NIA under contractual agreement with Harlan Sprague Dawley, Inc. (Indianapolis, IN). Cultures of adherent bone marrow cells were maintained as previously described [Lazarenko et al., 2007]. After 10 days of growth, cells were treated with Rosi (1 μM) or vehicle (DMSO) for 3 days followed by RNA isolation [Lazarenko et al., 2007].
Microarray experiments were performed as described previously [Lecka-Czernik et al., 2007]. Briefly, RNA quality was assessed using the Agilent Model 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). For each microarray five micrograms of total RNA was processed using the Affymetrix GeneChip® one-cycle target labeling kit (Affymetrix, Inc., Santa Clara, CA) according to the protocol recommended by the manufacturer. The resulting biotinylated cRNA was fragmented and hybridized to the GeneChip® Mouse Genome 430 2.0 Array (Affymetrix, Inc.). The arrays were washed, stained, and scanned using the Affymetrix Model 450 Fluidics Station and Affymetrix Model 3000 scanner by the University of Iowa DNA Core Facility according to the protocols recommended by the manufacturer. Probe intensity data were generated using the microarray suite (MAS) v5.0 software (Affymetrix, Inc.). The data is deposited in NCBI Gene Expression Omnibus (GEO) with accession number GSE10192.
The 2×2×3 factorial design included the treatment (+/- Rosi), cell type (+/- PPARγ2) and time (2, 24, and 72h) dimensions (see Fig. 1). Probe intensity data from all 24 Affymetrix GeneChip® arrays was read into the R software environment [Ihaka and Gentleman, 1996] directly from .CEL files using the R/affy package [Gautier et al., 2004]. Raw data quality was assessed using image reconstruction, histograms of raw signal intensities and MvA plots. Normalization was carried out using the robust multi-array average (RMA) method using all probe intensity data sets together [Irizarry et al., 2003] to form one expression measure per gene per array. Briefly, the RMA method was used to adjust the background of perfect match (PM) probes, apply a quantile normalization of the corrected PM values and calculate final expression measures using the Tukey median polish algorithm.
Log-transformed expression measures were expressed in fixed effects ANOVA models as the sum of different components contributing to the overall intensity value of each gene on the array [Churchill, 2004; Kerr et al., 2000]. First, the model:
was fit to the log-transformed gene expression measures Yi across each of the twelve treatment conditions, where μ is the mean for each array, CONDITION is the effect for treatment condition (e.g., +Rosi in U33/γ2 cells at 72h) and εi captures random error. A modified F statistic (Fs) that assesses differential expression between treatment conditions was used as a statistical filter to select one probe set for each mapped Entrez gene on the array [Cui et al., 2005].
Other models were used to test focused research hypotheses. In these cases, the data was subset into the arrays corresponding to the experimental state (i.e., a specific condition/time combination) to be tested. For each time point (2, 24, 72h) a model was fit:
where STATE refers to a combination of cell type (i.e., U-33/c or U-33/γ2) and activation status (i.e., with or without Rosi). All statistical tests were conducted with Fs, a modified F-statistic incorporating shrinkage variance components that allows variance estimates to include information from all the probe sets on the array [Cui et al., 2005]. Critical p-values were calculated by permuting model residuals 1000 times and pooling F-statistics [Yang and Churchill, 2007]. False discovery rate (FDR) values were estimated for each test result by implementing the “q-value” method of Storey [Storey, 2002]. A probe set was considered to be differentially expressed for false discovery rates less than 0.01, unless otherwise noted.
k-means clustering was used to identify groups of differentially expressed transcripts that respond in a concerted manner. The k-means approach minimizes variation within a cluster and maximizes variation between clusters. A gene was assigned to a cluster if it was 70% stable over 1000 bootstrap samples [Kerr and Churchill, 2001]. The proportion of variance explained was used to choose an appropriate number of clusters. Here, η represents the ratio of the between groups sum of squares to the total sum of squares, where η2 reduces to the correlation of multiple correlation (R2) for linear relationships. η2 was plotted versus k (see Fig. S1) in order to locate a value of “k” in an “elbow” structure with approximately 90% variance explained. Heat maps of clustering results were visualized in TreeView [Saldanha, 2004].
Over-represented classifications of genes were determined from statistical outcomes by testing for association with ‘biological process’ gene ontology terms [Ashburner et al., 2000]. Mappings between Affymetrix probe sets, Entrez gene identifiers and GO terms were based on mouse Build 36 retrieved from the R/mouse4302 package built on March 15, 2007 (www.bioconductor.org). One Entrez identifier per probe set was used in the analysis, based on the maximum Fs statistic calculated using equation . Hypergeometric tests for over-representation of GO terms were performed using the R/GOstats package [Falcon and Gentleman, 2007]. Unless specifically noted in the text, all enrichment analyses were based on a hypergeometric p < 0.01 significance criterion for groupings enriched with at least 5 genes.
We used a novel algorithm to quantify the hierarchical status of biological processes in the GO tree that we call a “modified GO slim” analysis. Our approach assigns a root to the GO tree based on either the presence of a gene in a term that is a common ancestor of two selected “GO slim” terms or a child of two “GO slim” terms (see Fig. S2). Using this approach, we found that the majority of the differentially expressed transcripts in our study were linked to one of ten biological process GO terms: “GO:0007165: signal transduction”, “GO:0007049: cell cycle”, “GO: 0002376: immune system process”, “GO:0006629: lipid metabolism”, “GO:0001501: skeletal development”, “GO:0006091: metabolic energy generation”, “GO:0008152: metabolism | GO:0006629: lipid metabolism” [i.e., GO:0008152 without (or conditional on) GO:0006629], “GO:0016043: cell organization and biogenesis”, “GO:0051301: cell division” and “GO:0008283: cell proliferation” (see Fig. S3, Table S1). To increase interpretability, we excluded the following “unknown” GO terms (“GO: 0050874”, “GO: 0051242”, “GO: 0050791”, “GO: 0051244”, “GO: 0007582”, “GO: 0043119”, “GO: 0051243”, and “GO: 0050875”) and nonspecific terms (“GO:0048522: positive regulation of cellular process, “GO:0009605: response to external stimulus”, “GO:0009719: response to internal stimulus”, “GO:0048519: negative regulation of biological process”) from consideration.
The most significant biological categories of genes pertaining to changes due to activation of PPARγ2 with Rosi (false discovery rate < 0.01) were determined by testing for association with high-level functions available in an expert-curated database of biological networks (Ingenuity Pathways Analysis™). The Ingenuity Pathways Knowledge Base (www.ingenuity.com) consists of expert-curated information from over 400 journals with known biological relationships between genes and gene products. A right-tailed Fisher's exact test for 2×2 contingency tables was used to determine the significance of overrepresentation of pathway members.
To validate gene expression results obtained from the microarray chips, an analysis of the expression levels of several chosen gene markers was performed using quantitative real-time PCR. This analysis was carried out with the same RNA samples that were analyzed on microarrays, as well as samples from two additional independent experiments performed in U-33/γ2 and U-33/c cells. In addition, the expression of the same gene markers was analyzed in primary bone marrow cells treated with Rosi or vehicle, as described above. Quantitative real time PCR was performed as described previously [Lazarenko et al., 2007]. Table S2 lists pairs of primers used for real time PCR validation of gene expression, whereas Table S3 shows the result of gene expression analysis. The expression of analyzed transcripts closely followed the expression detected on the microarrays and never exhibited change in an opposite direction (versus the microarray) or gave rise to a change detected on the microarray that was not confirmed by PCR.
A search for peroxisome proliferator response elements (PPREs) was carried out using the Target Explorer program [Sosinsky et al., 2003]. DNA sequences spanning 8.0 kb regions upstream of the transcription start sites of the selected genes were obtained from GenBank database (National Center for Biotechnology Information, Bethesda, MD). A set of 15 characterized PPRE sequences was used to create a mononucleotide position weight matrix [Zandbergen et al., 2005]. The following PPREs were used: acyl-CoA oxidase (TGACCTTTGTCCT and TGACCTTCTACCT), fatty acid binding protein 4 (TGAACTCTGATCC), apolipoprotein I (TGACCCCTGCCCC), apolipoprotein II (CAACCTTTACCCT), liver fatty acid binding protein (TGACCTATGGCCT), enoyl-CoA hydratase-3hydroxyacyl-CoA dehydrogenase (TGAACTATTACCT), HMG-CoA synthase (AGACCTTTGGCCC), lipoprotein lipase (TGCCCTTTCCCCC), cytochrome P450IVA6 (TCACTTTTGCCCT, TGGCCTTTGTCCT, and TGACCTTTGCCCA), acyl-CoA synthetase (TGACTGATGCCCT), malic enzyme (TCAACTTTGACCC), and G0/G1 switch gene 2 (TGACCTTTGCAAT). Briefly, the algorithm aligns input sequences, looks for conservation of bases at certain positions and translates the alignment into a matrix. A positive weight implies that the frequency of a nucleotide at a given position is higher than the a priori probability of this base at that position. Motifs with a score of 6.0 or higher were recorded. The threshold score was chosen arbitrarily to accommodate less specific PPREs. Random sequences were generated using FaBox random sequence generator [Villesen, 2007]. Each of these sequences was 8 kb long and had a 45% GC content, which was similar to the length and GC content of the regulatory regions of the analyzed genes.
A system of U-33/γ2 and U-33/c cells represents a cellular model of MSC differentiation under the control of PPARγ2 transcription factor [Lecka-Czernik et al., 1999]. In basal conditions, both types of cells retain the parental osteoblastic phenotype of U-33 cells. However, treatment with the PPARγ agonist Rosi commits U-33/γ2 cells, but not U-33/c cells, towards adipocytes and irreversibly suppresses their osteoblastic phenotype [Lecka-Czernik et al., 1999]. Hence, comparing the transcriptomes of U-33/γ2 and U-33/c cells in the presence or absence of Rosi allowed us to study the exclusive effects of PPARγ2 on marrow MSC lineage allocation.
Our purpose was to find time-dependent relationships in gene expression due to PPARγ2-induced changes in cell phenotype. We addressed three principal questions. First, we sought to characterize the sequence of events which led from osteoblastic cells to cells of adipocytic lineage by Rosi-activated PPARγ2. Next, based on our previous work suggesting that PPARγ2 suppresses osteoblastogenesis and induces adipogenesis by distinct mechanisms [Lazarenko et al., 2006; Lecka-Czernik et al., 2002], we sought to identify gene candidates mediating PPARγ2 anti-osteoblastic activity versus genes mediating pro-adipocytic activity of this transcription factor. Finally, we asked whether PPARγ2 possesses either ligand-independent or endogenous ligand-dependent activities that influence the expression signature of marrow MSC.
Here, U-33/γ2 and U-33/c cells were cultured in the presence or absence of Rosi and gene expression was monitored at three different time points (2, 24, 72h) after exposure to the agonist (Fig. 1). Each time point corresponds to a separate stage of Rosi-treated U-33/γ2 cell conversion from the osteoblast-like phenotype to the adipocyte-like phenotype and includes induction (2h), intermediate alterations in phenotype progression (24h), and a terminally differentiated adipocytic phenotype with completely suppressed osteoblastic characteristics (72h).
To begin, we defined a gene universe from the 45,037 probe sets on the microarray by mapping the putative transcripts to a total of 20,831 unique Entrez gene identifiers (see Experimental Procedures). Next, pairwise contrasts between cell states were conducted at each time point to form three distinct contrast sets (Fig. 2). The expression patterns in Set 1 (U-33/γ2 + Rosi vs. U-33/γ2) and Set 3 (U-33/γ2 vs. U-33/c) provide a basis for addressing our principal research questions, while Set 2 (U-33/c + Rosi vs. U-33/c) is used as a control comparison.
A total of 4,252 and 7,928 genes were found to be differentially expressed at one or more time points in Set 1 and Set 3, respectively. The time-specific overlaps within both contrast sets show an abundance of differentially expressed transcripts after 72h. In Set 1, a total of 147 (51), 822 (1,071), and 1,634 (2,096) genes are up (down-) regulated after 2, 24 and 72h, respectively. In Set 3, a total of 626 (510), 1,907 (2,267), 2,686 (3,355) genes were up (down-) regulated after 2, 24 and 72h, respectively. By contrast, only the expression of adiponectin (adipoq,+38.7 at 72h) was changed in Set 2. An increased expression of adiponectin, which is positively regulated by PPARγ [Sharma and Staels, 2007], in the absence of other transcriptional responses suggests either Rosi nonspecific effect or restricted activation of PPARγ1 isoform, which is naturally expressed in U-33/c cells. Nevertheless, the relative lack of differential expression in Set 2 supports our experimental methodology, since the ligand under investigation should not affect gene expression unless the studied transcription factor is present.
Seven clusters explained approximately 88% of the variance in transcriptional patterns due to time-dependent treatment of U-33/γ2 cells with Rosi (Fig. 3A and B, and Table S4). Genes in two clusters (Clusters 1 and 6) were increasingly upregulated over time with Rosi but not without Rosi. These clusters include genes associated with processes related to lipid and energy metabolism. Two clusters (Clusters 5 and 3) showed genes that were up- (or down-) regulated after 24h and returned to their early expression level after 72h treatment with Rosi while steadily increasing (or decreasing) without Rosi. Cluster 5 includes genes associated with immune processes, skeletal development, cell organization and biogenesis, and cell proliferation. Genes in Cluster 3 are associated with cellular biogenesis, lipid metabolism and energy generation. Finally, genes in the remaining three clusters were mostly downregulated over time with Rosi, but showed increasing (Cluster 7) or decreasing (Clusters 2 and 4) expression levels without Rosi. These genes are associated with metabolism and cell proliferation (Cluster 7), cell division, cell cycle and biogenesis (Cluster 2), and metabolism, skeletal development, and signal transduction (including TGFβ and G-protein signaling, and regulation of MAPK activity) (Cluster 4).
At the beginning of treatment, cells in both cultures exhibited 80% confluency and were exposed to fresh medium. Therefore, the pattern of genes differentially expressed at 2h should reflect basic differences between the two cell lines during the late exponential phase of growth. Although both cell types experienced slower growth after 24h, this intermediate time point reflects basal differences between cell lines and changes due to growth inhibition. However, after 72h both types of cells were in the confluent stage, possessing a fibroblast-like morphology and a termination of cell division. Thus, in contrast to 2h and 24h, the 72h time point reflects phenotypic differences between cell lines, without the confounding effect of growth rate.
Seven clusters explain about 89% of the variance in differential transcription corresponding to Set 3 (Fig. 4A and B, and Table S5). Clusters 4 and 7 contain genes with expression levels that are unchanged without PPARγ2 but increased dramatically with PPARγ2. Significant biological processes associated with these clusters include metabolism and inflammatory response (Cluster 4) and cell biogenesis and signal transduction including G-protein coupled receptor protein signaling (Cluster 7). Clusters 1, 3, 5, and 6 are characterized by significant decrease in gene expression in the presence of PPARγ2 and modest changes in gene expression in the absence of PPARγ2. Genes in these clusters are associated with lipid metabolism, biogenesis, cell cycle and intracellular signaling, including regulation of MAPK activity and apoptosis. Finally, genes in Cluster 2 show increasing expression with and without PPARγ2, but a higher increase over time without PPARγ2. These genes are involved in regulation of transcription regulation of metabolism and bone remodeling.
The Ingenuity Pathways Knowledge Base (IPKB) was used to classify over-represented genes into high-level functions reflecting functional relationships between genes and gene products (Table 1). Genes found in Lipid Metabolism and Carbohydrate Metabolism categories are induced early and are consistently upregulated, reflecting PPARγ2 proadipocytic activity and its role in the control of energy metabolism. The expression of genes in the Skeletal and Muscular System Development and Function category is affected in a manner that is consistent with the known function of these genes and anti-osteoblastic activity of PPARγ2. Positive regulators of osteoblast development are downregulated (e.g. Wisp1, Dmp1 and BMP4), whereas genes with a known or potential role in the negative regulation of this process are upregulated (e.g. Tob1 and Tle3) (see Table 4). In addition, Rosi affects the expression of a significant number of genes controlling cell growth and proliferation; both processes are required for osteoblast differentiation. Interestingly, besides gene categories closely related to fat and bone cell differentiation, PPARγ2 also regulates the expression of a large number of genes controlling cell death and cancer, as well as connective tissue, immune, nervous and cardiovascular systems development.
A modified GO slim analysis of the directed acyclic subgraph (see Experimental Procedures) reveals six important categories of biological processes that are affected at 2h, including morphogenesis, signal transduction, lipid metabolism, apoptosis, immune response, and regulation of cell differentiation (Fig. 5 and Table 2). Most genes involved in cellular lipid metabolism (node 34) and negative regulation of apoptosis (node 53) were induced at 2h and remained induced over 24h and 72h (Fig. 5 and Table 2). The expression of genes involved in the regulation of cell differentiation (node 30) was predominantly increased at all time points, however this node included adipocyte-specific, but not osteoblast-specific, genes. Interestingly, genes corresponding to the Wnt-receptor signaling pathway (node 40) show initial upregulation followed by significant downregulation at 72 h. Components of the immune system process were induced early, however myeloid cell differentiation (node 52) was upregulated at all time points, whereas T cell activation (node 35) and lymphocyte differentiation (node 51) showed diminished expression after 24h and 72h, respectively (Fig. 5 and Table 2).
Considerations of the expression of all 198 early responder genes exhibited three distinctive patterns of expression over time (FDR < 0.01). Table S7 provides an entire list of genes, whereas Table 3 lists only selected examples of genes, with either known or potential role in the regulation of adipocytic and osteoblastic phenotype.
The first group of 66 genes was characterized by consistently increased expression indicating their requirements for differentiation and the maintenance of adipocyte phenotype and function. The majority of these genes represent transcriptional regulators of adipocyte differentiation (e.g. C/EBPα and STAT5a) and genes involved in adipocyte function including adipokine production, lipid accumulation and lipid metabolism (e.g. Acox1, Adipoq, CD36, Cidec, Lipe, Lpin1, FABP4, Pnpla2, and Pex). A number of genes in this group with less characterized function: Angptl4, Aqp7, G0s2, and Ubd, achieved a very high level of expression in mature adipocytes (72h treatment) suggesting their role in the maintenance of adipocyte phenotype. This group also includes Tle3 and Tob1, negative regulators of the pro-osteoblastic Wnt and TGFβ/BMP signaling pathways, respectively [Javed et al., 2000; Yoshida et al., 2000].
A second group of 19 genes show consistently decreased expression in response to Rosi-activated PPARγ2 in comparison with non-treated cells. This group includes known and potential positive regulators of osteoblastogenesis such as Wisp1, BMP4, Plk2, and Fosl2 (Fra-2) [Abe et al., 2000; French et al., 2004; Ma et al., 2003; McCauley et al., 2001], as well as Socs5, Ccl2 and RANKL, cytokines involved in osteoblast communication with hematopoietic environment [Boyce and Xing, 2008; O'Shea and Murray, 2008].
The third group of 113 genes is characterized by rapid upregulation or downregulation at the 2h time point followed by a fluctuating pattern of expression. Two important regulators of phosphate homeostasis and mineralization, Phex1 and Dmp1 [Strom and Juppner, 2008], responded early to PPARγ2 with a pattern of expression over the time of treatment suggesting their involvement in the induction step of osteoblast to adipocyte conversion. For instance, after initial downregulation Dmp1 expression returned to the basal level, whereas Phex1 expression after initial upregulation was downregulated by 20-fold in differentiated adipocytes. Similarly, PPARγ2 affected expression of number of genes representing G-protein signaling family, which plays an important role in osteoblast differentiation [Hsiao et al., 2008; Peng et al., 2008; Teplyuk et al., 2008]. PPARγ2 affects early and temporally the expression of G-protein coupled receptor Gprc5a, Rho GTPase activating proteins Arhgap5 and 22, members of RAS oncogene family (Rab 5b and 30), and Rsg4 regulator.
Some of the genes, which are rapidly and temporally upregulated by Rosi may play a dual role as mediators of both, anti-osteoblastic and anti-diabetic activities of PPARγ. PIK3R1, the p85α regulatory subunit of phosphatidylinositol 3-kinase, plays an important role in insulin sensitivity and glucose homeostasis and mediation of growth factors signaling, including IGF-1 [Hallmann et al., 2003]. Socs2 negatively regulates growth hormone action and functions as a mediator of cross talk between PPARγ and growth hormone in vivo [Rieusset et al., 2004]. Animals deficient in Socs2 display an increased longitudinal skeletal growth but reduced bone mineral density, as a consequence of deregulated GH/IGF-1 signaling [Lorentzon et al., 2005]. Ppap2b, phosphatidic acid phosphatase type 2b, regulates phospholipids metabolism and might be involved in the production of PPARγ endogeneous ligands [Davies et al., 2001; Kanoh et al., 1999]. Ppap2b also interacts with Wnt pathway by inhibiting β-catenin-mediated TCF transcriptional activity [Escalante-Alcalde et al., 2003]. Thus, increased expression of Ppap2b may have two potential consequences: it may upregulate production of endogenous PPARγ ligand and/or it may contribute to downregulation of Wnt pathway pro-osteoblastic activity.
Genes with increased expression after 2h of Rosi treatment are candidates for direct regulation by PPARγ. Transcriptional activity of PPARγ involves binding to a specific regulatory sequence PPRE, which is present in upstream regions of genes regulated by this nuclear receptor. In silico analysis of 8 kb fragments located upstream of transcription start sites (TSS) for genes early regulated by Rosi-activated PPARγ2 revealed a high number of PPREs frequently organized in clusters, as compared to the occurrence of PPREs in randomly generated DNA sequences of the same length (Bonferroni-corrected t-test p=0.00004) and upstream regions of genes characterized by Rosi-suppressed expression (t-test p=0.004). Genes characterized by high and increasing expression over time show 6-fold higher concentration of PPREs within 1kb upstream of the TSS, as compared to random sequences and upstream regions of Rosi-suppressed genes (Fig. S4). In contrast, the frequency of PPREs in upstream regions of suppressed genes was not statistically different from the frequency recorded for random sequences (t-test p=0.1).
PPARγ2 transcriptional control involves upregulation and downregulation of gene expression. Locating multiple PPREs clustered in the regulatory regions of upregulated genes of the adipocytic pathway suggests direct transcriptional control by PPARγ2 and a modular organization of cis-regulatory regions with multiple transcription factor binding sites [Heinaniemi et al., 2007]. Similarly, the lack of PPREs in the regulatory regions of downregulated genes suggests a lack of PPRE involvement. This finding is consistent with reports that PPARγ has a suppressive effect on gene expression without direct transcriptional control [Ricote and Glass, 2007].
Osteoblast differentiation and function is regulated by a network of signaling pathways, phenotype-specific transcription factors, enzymes, and structural genes that permit the formation of a collagen based extracellular matrix and its subsequent mineralization [Lian et al., 2006]. Osteoblasts also produce factors, such as macrophage colony stimulating factor (M-CSF), receptor activator of the NFκB ligand (RANKL) and its decoy receptor osteoprotegerin (OPG), which support the development of bone resorbing cells, or osteoclasts [Boyce and Xing, 2008].
Rosi-activated PPARγ2 suppresses the expression of many genes linked to signaling, transcriptional regulation of osteoblast differentiation, and extracellular matrix formation and mineralization (Table 4, Fig. 6). Early changes in gene expression are seen among members of Wnt and TGFβ/BMP signaling pathways, followed by a suppression of increasing number of pathway-specific genes with time of treatment. The pattern of gene expression specific for these pathways at the 72h time point suggests that in cells converted to adipocytes Wnt and TGFβ/BMP pathways activities are substantially suppressed. Transcripts for essential osteoblast-specific transcription factors, such as Dlx5, Runx2, and Osterix, showed altered expression only after 24h treatment with Rosi and continued to be downregulated in cells converted to adipocytes. This pattern of expression suggests that the anti-osteoblastic PPARγ2 effect is mediated through early responders, but not osteoblast-specific transcriptional regulators.
The effect of Rosi on osteoblast specific regulators of transcription and function is followed by changes in the expression of genes encoding large number of proteins involved in extracellular matrix formation and osteoblast function. These include different types of collagen, fibronectin, integrins, alkaline phosphatase, osteocalcin, biglycan and bone sialoprotein. The large number of genes coding for G-protein family is affected over all time points of treatment. Through their roles in cytoskeleton organization and integration of cellular signaling, this group of proteins may represent potential regulators of osteoblast-to-adipocyte conversion.
Changes in the expression of other osteoblast-specific signaling pathways, including FGF, IGF-1, and PTH, either coincide with the effect on transcriptional regulators or occur later. This pattern indicates the relative distance of these signaling pathways from direct transcriptional control by PPARγ2.
The effect of Rosi on the expression of M-CSF, RANKL, VDR and OPG indicates that osteoblast support for osteoclastogenesis is under a direct control of PPARγ2. Even though PPARγ2 decreases the expression of pro-osteoclastic cytokines and the vitamin D receptor here, we have previously shown that this Rosi-induced anti-osteoclastic effect is completely abrogated in U-33/γ2 cells in the presence of 1,25-dihydroxyvitamin D [Lazarenko et al., 2007]. Additionally, Rosi has also been shown to induce RANKL expression in primary bone marrow cells, and Rosi administration to mice with attenuated bone formation increased bone resorption [Lazarenko et al., 2007]. Therefore, in some instances in vitro results may reflect regulatory associations between genes rather than defined transcriptional controls of their expression. In this light, our microarray results indicating control of osteoclastogenesis should be interpreted in light of the biological context.
A question-driven multivariate approach was used to probe the mechanistic consequences of Rosi-activation of PPARγ2 as cells convert from an osteoblast-like to a mature adipocyte phenotype. In contrast to previous studies, which examined a static effect of Rosi-activated PPARγ2 on marrow MSC gene expression [Lecka-Czernik et al., 2007; Lecka-Czernik et al., 1999; Lecka-Czernik et al., 2002; Shockley et al., 2007], the analysis presented here describes PPARγ2-directed temporal relationships in gene expression during progression toward adipocytic and suppression of osteoblastic phenotype. The following aspects of this process were considered: 1) the sequence of events which led from osteoblastic cells to cells of the adipocytic lineage, 2) whether PPARγ2-triggered pro-adipocytic and anti-osteoblastic mechanisms share common regulatory steps, and 3) whether PPARγ2 exerts transcriptional activities in the absence of exogenous activator.
A set of early responder genes provides information about initial triggers directing cells toward the adipocyte and away from the osteoblast lineage. Our analysis indicates that the pro-adipogenic mechanism includes PPARγ2 direct transcriptional control through PPREs located in the promoter region of other transcriptional regulators and genes involved in adipocyte function. While adipocytic-induction occurred very early and continued to increase through mature fat cells, most osteoblast-related genes were downregulated after initiation of adipocyte differentiation. The exceptions include several genes, which are good candidates for initial mediators of PPARγ2 anti-osteoblastic activity. These candidates consist of two transcripts encoding negative regulators of osteoblastogenesis (Tle3 and Tob1), and several transcripts which encode for known or potential positive regulators of osteoblastogenesis (e.g. Wisp1, BMP4, and Dmp1, Phex1). Interestingly, these genes function in the regulation of extracellular signaling. Thus, in contrast to the pro-adipocytic mechanism which involves direct effect on transcriptional regulators, the anti-osteoblastic mechanism is triggered by an effect on signaling pathways which, in turn, downregulates the expression osteoblast-specific transcriptional regulators. Without ruling out cross-talk between processes, we postulate that different mechanisms drive PPARγ2-mediated pro-adipocytic and anti-osteoblastic effects. This proposition is consistent with previous evidence that osteoblastic repression and adipocytic induction can be distinguished and are independently regulated [Lecka-Czernik and Suva, 2006].
The late effect of PPARγ2 on the activity of other pathways, such as FGF, IGF-1, and PTH, indicates their relative distance from the PPARγ2 direct control and suggests that their suppression results from either PPARγ2 negative effect on other signalings or osteoblast-specific transcriptional regulators. This pattern implies a hierarchical structure of PPARγ2-initiated negative regulation of osteoblast differentiation and suggests that PPARγ2 anti-osteoblastic activity involves pleiotropic interactions between signaling processes and transcriptional regulators over time.
Differences between U-33/γ2 and U-33/c cell types underscore the activity of PPARγ2 even without an activating ligand. Our results demonstrate that PPARγ2 alone modulates mesenchymal stem cell phenotype in a cell culture model. It suppresses the cell cycle response and affects the expression of genes specific for the maintenance of stem cell phenotype. We have previously showed that PPARγ2 suppresses the expression of “stemness” genes including Lif and Lif receptor, while increasing the expression of genes supporting hematopoiesis including Kitl, RANKL, and the ligands for CXC chemokines [Shockley et al., 2007]. In light of increased PPARγ2 expression in aged marrow MSCs, these findings suggest that PPARγ2 may be a factor responsible for age-related changes in MSCs differentiation potential toward adipocytes and osteoblasts [Moerman et al., 2004].
In conclusion, our study shows that PPARγ2 is a major regulator of bone marrow MSC differentiation and provides information about hierarchical interactions between different regulatory pathways involved in fat- and bone-cell development. Marrow-derived adipocytes produce an expression signature that is similar to known adipocytic cells with respect to fatty acid metabolism, carbohydrate metabolism, and adipokine production. Most importantly, we provide a database of triggers and temporal expression patterns that should be helpful in our search to identify regulatory mechanisms by which PPARγ2 controls osteoblast differentiation.
Table S7. “Early responders” to Rosi treatment (2 hr) of U-33/γ2 cells
Grant support To BLC - Contract grant sponsor: NIH/NIA/NIAMS; Contract grant numbers: AG17482 and AG028935; Contract grant sponsor: American Diabetes Association; Contract grant number: 1-03-RA-46;
To KRS - Contract grant sponsor: NIH/NHGRI Ruth L. Kirchstein Postdoctoral Fellowship; Contract grant number: HG003968.