Temporal Gene Expression Profiling during Rat Femoral Marrow Ablation-Induced Intramembranous Bone Regeneration

Enhanced understanding of differential gene expression and biological pathways associated with distinct phases of intramembranous bone regeneration following femoral marrow ablation surgery will improve future advancements regarding osseointegration of joint replacement implants, biomaterials design, and bone tissue engineering. A rat femoral marrow ablation model was performed and genome-wide microarray data were obtained from samples at 1, 3, 5, 7, 10, 14, 28, and 56 days post-ablation, with intact bones serving as controls at Day 0. Bayesian model-based clustering produced eight distinct groups amongst 9,062 significant gene probe sets based on similar temporal expression profiles, which were further categorized into three major temporal classes of increased, variable, and decreased expression. Osteoblastic- and osteoclastic-associated genes were found to be significantly expressed within the increased expression groups. Chondrogenesis was not detected histologically. Adipogenic marker genes were found within variable/decreased expression groups, emphasizing that adipogenesis was inhibited during osteogenesis. Differential biological processes and pathways associated with each major temporal group were identified, and significantly expressed genes involved were visually represented by heat maps. It was determined that the increased expression group exclusively contains genes involved in pathways for matrix metalloproteinases (MMPs), Wnt signaling, TGF-β signaling, and inflammatory pathways. Only the variable expression group contains genes associated with glycolysis and gluconeogenesis, the notch signaling pathway, natural killer cell mediated cytotoxicity, and the B cell receptor signaling pathway. The decreased group exclusively consists of genes involved in heme biosynthesis, the p53 signaling pathway, and the hematopoietic cell lineage. Significant biological pathways and transcription factors expressed at each time point post-ablation were also identified. These data present the first temporal gene expression profiling analysis of the rat genome during intramembranous bone regeneration induced by femoral marrow ablation.


determined th
t the increased expression group exclusively contains genes involved in pathways for matrix metalloproteinases (MMPs), Wnt signaling, TGF-b signaling, and inflammatory pathways.Only the variable expression group contains genes associated with glycolysis and gluconeogenesis, the notch signaling pathway, natural killer cell mediated cytotoxicity, and the B cell receptor signaling pathway.The decreased group exclusively consists of genes involved in heme biosynthesis, the p53 signaling pathway, and the hematopoietic cell lineage.Significant biological pathways and transcription factors expressed at each time point post-ablation were also identified.These data present the first temporal gene expression profiling analysis of the rat genome during intramembranous bone regeneration induced by femoral marrow ablation.

Introduction

Bone is a dynamic organ that undergoes continuous remodeling by contr lled cycles of bone resorption and bone formation, which are balanced to preserve bone mass.In the case of common metabolic bone disorders, such as osteoporosis, reduction in skeletal mass is caused by an imbalance between bone resorption and bone formation.Both types of bone regeneration, intramembranous and endochondral in response to fracture healing, are generally known to have parallels with developmental bone formation, and involve distinct yet interdependent healing phases consisting of biologically complex processes regulated by a very large number of transcriptional events [1][2][3].Bone marrow ablation in long bones induces intramembranous bone formation and subsequent bone resorption in order to regenerate normal bone marrow, and was originally established as an experimental model to study hematopoiesis [4,5].

Several groups, including ours, have used the rat marrow ablation model for investigations related to implant fixation [6][7][8][9][10][11][12].The marrow ablation model has been utilized further as an e perimental model to study intramembranous bone regeneration with concentration on histological and biochemical approaches and focused gene expression analysis [2,[13][14][15][16]. From these investigations, it is understood that following marrow ablation there are three major and distinct, yet overlapping, phases of healing which can be generally described.The first phase primarily consists of clot formation and inflammation from days 1 to 5. The second major phase is that of repair from day 3 to 14, and involves neovascularization, erivascular maturation, mesenchymal stem cell migration, proliferation and condensation, osteoblastic differentiation, and woven bone formation.Lastly, there is a remodeling phase from approximately days 10 to 28, until restoration of fatty and hematopoietic marrow is achieved by 56 days.

Although the distinct phases of inflammation, bone repair, and bone remodeling following rat femoral ablation surgery can be well-defined by histological methods, enhanced understanding of temporal gene expression profiling and identification of significant biological pathways associated with the distinct phases of intramembranous bone regeneration will greatly improve future advancements of fixation and osteointegration of joint replacement implants [17], design and synthesis of novel biomaterials, and bone tissue engineering [18].

Therefore, the mechanical ablation of the rat femoral marrow cavity is an established and suitable model for gene expression profiling studies regarding intramembranous bone regeneration.Our group has previously used this rat femoral marrow ablation model to characterize coexpression patterns of 39 genes during repair phases of intramembranous bone regeneration up to 14 days [2] and to report modulation of expression of 21 osteogenic genes following local application of rhTGF-b2 [19].To date, there have been no published reports of genome-wide temporal transcriptional analysis of intramembranous bone regeneration that takes place in a marrow ablation model.A recent report utilized microarray data and temporal transcriptio

the mouse
ranscriptome of an endochondral bone formation process occurring during a 21 day period of fracture healing [20].Certain gene expression data analysis methods used for that study are similarly used in the present study, including a Bayesian modeling approach to cluster temporal gene expression profiles providing by the Cluster Analysis of Gene Expression Dynamics Program (CAGED), and specific analysis to identify significant biological processes and pathways using the DAVID bioinformatics resources.In this study, we identify significant pathways including matrix metalloproteinases (MMPs), Wnt signaling, TGF-b signaling, and notch signaling for major temporal expression groups, as well as differential pathways and transcription factors expressed at each time point up to 56 days.


Materials and Methods


Rat Model

In an Institutional Animal Care and Use Committee of the Rush University Medical Center (IACUC; Protocol #06-005) approved study, 45 adult male rats (Sprague Dawley, 400-425 g) were divided into nine groups: intact control (0 day), and 1, 3, 5, 7, 10, 14, 28, and 56 days post marrow ablation.Forty animals received unilateral femoral ablation, adapting the method described previously [2,14,15].Briefly, the rats were anesthetized by intraperitoneal injection of ketamine hydrochloride (100 mg/ kg) and xylazine (5 mg/kg) supplemented as necessary.Surgery was performed with adherence to aseptic technique.The hindlimb was shaved and scrubbed with ethanol and betadine solution.An incision (,1 cm) was made along the medial aspect of th

patella, and the pate
la along with the quadriceps tendon and patellar ligament were retracted to expose the distal condyles of the femur.A 2.0 mm hole was drilled through the patellar surface of the femur to gain access to the medullary canal.The contents of the marrow cavity were disrupted mechanically and reamed by hand with a 2.0 mm drill bit proximally up to the lesser trochanter.The canal was then irrigated with 10 ml of saline.The distal opening hole in the bone was sealed with bone wax, the patella was repositioned, and the deep fascia and skin sutured separately.Following surgery, the locomotion, grooming, and eating habits of the surgically treated rats were normal.At each time point, animals were killed in a carbon dioxide chamber, and whole femurs were harvested and denuded of soft tissue.


Total RNA Extraction

For three samples per time point, both the distal and proximal ends of the femurs were cut off in order to exclude the epiphyses and growth plate regions.Diaphyseal marrow t

sue and cells from the mid-shaft were flushed out
sing TrizolH.Following homogenization, 1 ml of solution was transferred to a 1.5 ml Eppendorf tube and centrifuged at 12,000 g for 10 minutes at 4uC to remove insoluble material.The supernatant containing RNA was collected, mixed with 0.2 ml of chloroform, and centrifuged at 12,000 g for 15 minutes at 4uC.After RNA in th

aqueous phase was transferred into a new tube,
the RNA was precipitated by mixing 0.5 ml of isopropyl alcohol and recovered by centrifuging the tube at 12,000 g for 10 minutes at 4uC.The RNA pellet was washed briefly in 1 ml of 75% ethanol and centrifuged at 7,500 g for 5 minutes at 4uC.Finally, the total RNA pellet was dissolved in diethylpyrocarbonate (DEPC) water, and its quantity was assessed by spectrophotometric analysis.


Experimental Design for Gene Expression Profiling

A flow chart displaying the major steps involved in the experimental design is provided in Figure 1.Gene expression profiling was performed on total RNA samples that were collected at day 0, which is the biological reference of intact bone, and at 1,

5, 7, 10, 14, 28, and 56 days post-ablati
n of rat femoral marrow.


Microarray Hybridization and Data Acquisition

For each of the three samples per time point, gene expression was analyzed with GeneChipH Rat Genome 230 2.0 Arrays (Affymetrix, Santa Clara, CA), which is comprised of over 31,000 probe sets representing approximately 28,700 well-characterized rat genes.For each gene, eleven pairs of oligonucleotide probes are synthesized in situ on the arrays.Total RNA from the samples was hybridized using optimal reagents and standardized protocols, a GeneChipH Hybridization Oven, a GeneChipH Fluidics Station, and a GeneChipH Scanner 3000 enabled for High-Resolution Scanning.Data acquisition from the microarrays also required the use of Affymetrix GeneChipH Command ConsoleH Software (AGCC).


Microarray Expression Data Pre-processing

The raw gene expression data from Affymetrix were available in the form of binary CEL (Affymetrix cell intensity) files.The opensource integrated software system BRB-Ar

yTools (Dr.Richard Simon and BRB-Arr
yTools Development Team) [21] was used and the CEL files were collated with RMA (robust multi-array average) [22] method and ''affy'' R/Bioconductor package to compute probeset summaries [23].This utilized a three-step approach of background correction on PM (Perfect Match) data, quantile normalization, and Tukey's median polish algorithm for summarization of probe level data.The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [24] and are accessible through GEO Series accession number GSE22321 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token = hbybhuqmksgeydy&acc = GSE22321).Specific genes are mentioned by the gene name and/or official gene symbol, along with the Entrez Gene ID number in square brackets.


Microarray Expression Data Analysis

A univariate F-test (with random variance model) with a significance threshold of p,1610 23 (an appropriately stringent significant level to reduce the chance of false positives) was performed with BRB-ArrayTools and used to determine differentially expressed gene probe sets over all time points.Gene expression profiles were clustered with Cluster Analysis of Gene Expression Dynamics (CAGED version 1.1) program, which utilizes a Bayesian model-based clustering method on temporal gene expression data [25] and uses an agglomerative procedure to identify the most probable set of clusters, where genes assigned to certain clusters have similar temporal expression profiles.In turn, genes clustered together in this way are likely to share similar physiological functions or regulation.Data were normalized as ratios to the expression values on Day 0. A model order of 0 was used, where d ta from each time point are assumed to be independent from each other.The prior precision and gamma value were set to 1 and 0, respectively, where the prior precision is the size of the sample upon which the prior distribution is built, while the gamma value is the rate to zero of the prior precision, with 0 representing the case of perfect ignorance.A Bayes Factor of 1 was used to impose this minimum limit for accepting the merging of two clusters if the Bayes Factor of their merging is at least the value of 1.The method required a similarity measure to guide the search procedure and a Euclidean distance measure between gene expression profiles was adopted.Goodness of fit of the resulting model was assessed by checking the normality of the standardized residuals of each cluster.

Each cluster identified by CAGED analysis was assigned to one of the three major temporal groups according to increased, variable, or decreased expression, which were then further analyzed with DAVID EASE (version 2.0) [26] software to identify significant gene ontology categories or biological pathways.All significant gene categories associated with the Biological Processes domain of the Gene Ontology (GO) Consortium [27] were determined, using a significance threshold of p,0.05.Overall major categories of biological processes were formed by manually combining specific subcategory terms having related or overlapping functions.DAVID EASE was also used to analyze the gene sets for each major temporal group (increased, variable, decreased expression) to find significant known pathways determined by KEGG [28] and GenMAPP [29] databases, using a significance threshold of p,0.05.Cluster 3.0 [30] was used to cluster the log base 2 expression values (for each day time point vs. day 0 time point) of genes associated with certain significant biological pathways of interest from the three major temporal expression groups (increased, variable, and decreased expression).For each significant pathway identified, the number of significant genes from the initial significant gene list known to be associated with a given pathway is noted as ''Gene List Hits''.Additionally, the total number of genes on the Affymetrix GeneChipH Rat Genome 230 2.0 Array known to be associated with a given pathway is noted as ''Gene Total Hits''.The log base 2 fold-change ratios were clustered using hierarchical clustering with a centered correlation distance/similarity metric and average linkage clustering method.The clustered data table file was viewed in TreeView [31] using the pixel setting contrast default of 3 and using blue and red to represent positive and negative fold-change expression values, respectively.A univariate two-sample T-test with significance threshold of p,1610 23 n BRB-ArrayTools was further used to determine significant gene lists of gene probe sets differentially expression on each time point (day 1, 3, 7, 10, 14, 28, and 56 vs. day 0 intact control) post-ablation.DAVID EASE was used to analyze gene lists for each time point post-ablation and find all kn

determined by KEGG and GenMAPP databases using a significance threshold of p,0.05.Additionall
the total number of probe sets known to be transcription factors and present on the GeneChipH Rat Genome 230 2.0 Array was determined to be 1,254, and significant gene lists for each time point (each day vs. day 0) were compared with the total number of transcription factor list, and the number of significant transcription factors expressed for each time point was found.


Histological methods for imaging of femoral marrow samples at each post-ablation time point

For two animals per time point, the whole femur was dissected and fixed in 4% paraformaldehyde in phosphate-buffered saline.Tissues were decalcified with 0.5 M ethylenediaminetetraacetate (EDTA) and embedded in paraffin.Four-mm-thick sections were cut in the sagittal plane through the femur and stained with hematoxylin and eosin.Images from selected sub-metaphyseal sites were recorded at 10x and 40x magnification with a NikonH600L photomicroscope.


Results


Representative histological images of femoral marrow samples at each time point post-ablation

Stages of bone regeneration in the ablated intramedullary space are illustrated in Figure 2.This regenerative sequence occurs in spatially and temporally complex domains within the metaphyseal and diaphyseal regions.The process can be divided into eight subphases.Immediately following the ablation, a blood clot fills the space producing a fibrin and platelet laden substrate for a succession of cellular infiltrates (clot consolidation phase, days 1-3).The perimeter of this clot shows polymorphonuclear inflammatory infiltration associated with damaged penetrating vessels (inflammatory phase, days 3-5).Subsequently, there is a rapid invasion of fibrovascular progenitors that actively generate a collagenous matrix within the clot (granulation phase, days 4-7).Vasculogenesis and angiogenesis become evident in this primitive matrix structure wherein new capillary networks become connected to feeder vessels to establish a circulation (neovascularization phase, days 3-7).Perivascular infiltrates, including pericytes, lymphocytes, mast, and mononuclear cells become prominent as this primitive vascular network matures and establis

s flow (perivascular maturation phase, da
s 5-10).Giant cells become conspicuous within this perivascular domain.The extravascular collagenous domains increase in density with sites of mineralization characteristic of woven-bone that assumes trabecular form (osteogenic phase, days 5-10).Arrays of osteoblastic cells begin to occupy the woven bone surfaces to produce lamellar bone on the woven bone cores (trabecular maturation phase, days 7-14).Variable sites of osteoclastic activity accompany this trabecular maturation phase, finally tipping the balance to the remodeling and resorption of the intramedullary trabecular structure (resorption phase, days 10-28) as the fatty and hematopoietic marrow is reconstituted (marrow reconstitution phase, days 28-56).The stages represented in the figure are overlapping as new trabecular bone rapidly forms to occupy the intramedullary space by day 7 to 14 and, following resorption of much of this bone by day 28, the fatty and hematopoietic marrow composition is restored.
emporal clusters of expression profiles

The graphical representations of the average expression profiles for each of the 8 unique temporal clusters determined by cluster analysis of significant gene probe sets are shown in Figure 3.These clusters can be differentiated according to the pattern and In-situ histology of intramedullary bone formation at each post-ablation time point (10X & 40X).Intact Day 0, submetaphyseal regions of the intact bone characteristically are populated by adipocytes and hematopoetic cell populations among scattered trabecular profiles.Day 1, the post-ablation clot containing scattered polymophonuclear cells fills the marrow space.Day 3, the clot is beginning to show cellular infiltration by cells of indeterminate origin.Day 5, organization of the clot with cell transformations in a fibrovascular structure with small vessels and an immature collagenous network.Day 7 and Day 10, increases in collagenous interstitial matrix.Day 14, consolidation and modeling of pretrabecular matrix structure with expansion of cellularity in perivascular space of the maturing vascular network.Day 28 and Day 56, maturation of trabecular structure and marrow shows reconstitution of the pre-ablation tissue architecture.Osteocytes and bone-lining osteoblastic cells are evident.doi:1 .1371/journal.pone.0012987.g002degree of differential expression.Clusters 1, 7, and 8 were assigned to the increased expression group (Figure 3A), clusters 2, 5, and 6 were assigned to the variable expression group (Figure 3B), and clusters 3 and 4 were assigned to the decreased expression group (Figure 3C).Clusters were assigned to the variable group if the expression profiles of all genes contained in that cluster exhibited patterns of both increased and decreased expression temporally (namely cluster 2 and 5) or if the log base 2 fold-change of expression values were less than 0.5 for all time points (cluster 6).Table S1 provides detailed information and expression data regarding the gene probe sets comprising each of the 8 clusters.


Identification of significant biological process ontologies and pathways associated with major temporal groups of clustered gene expression profiles

Pie graphs (Figure 4) present the percentage distribution of statistically significant (p,0.05)biological process ontologies identified for increased, decreased, and variable expression groups.All individual biological ontology results were manually examined for related or overlapping functions in order for consolidation into overall major categories of biological process and this detailed information is summarized and shown for each major temporal group in Table S2A-D.For all three of the major temporal groups and therefore during all time points of intramembranous bone regeneration following marrow abl tion, certain biological processes comprise large percentages out of the total biological functions identified.These aforementioned biological processes were specific aspects of metabolism, a large variety of developmental related processes, cell cycle, cellular si naling, as well as a large category referred to as miscellaneous.In contrast, there are unique differences in the percentage distributions of specific biological functions identified in each of the three temporal expression groups, and these will be described henceforth.

The major biological processes identified exclusively for increased gene expressions, which all peaked between days 5 and 10 as shown by their average expression profiles (Figure 4A), consisted of a large variety of processes related to neurogenesis, skeletogenesis, cell motility, cell adhesion, vasculogenesis, and Wnt signaling, as supported by detailed information in Tables S2A   Figure 3. Expression profiles of temporal clusters during intramembranous bone regeneration.Clustering of expression profiles of the 9,062 statistically significant gene probe sets by Bayesian modeling produced 8 temporal clusters, which were further grouped into three major temporal groups of (A) Increased Expression (cluster #1, 7, and 8; graphed as green), (B) Variable Expression (cluster #2 5, 6; graphed as blue), and (C) Decreased Expression (cluster #3 and 4; graphed as red).Data is presented as log base 2 fold-change values over the time points of day 1, 3, 5, 7, 10, 14, 28, and 56 days post-ablation, and is graphed as a solid line.The cluster number and the number of genes for each clustered gene expression profile graph are denoted in the title.Note that the scale for the vertical axis is not constant.doi:10.1371/journal.pone.0012987.g003and S2B.Biological processes identified for the variable expression group (Figure 4B) were dominated mostly by elements of metabolism but also processes related to DNA replication, apoptosis, cytoskeleton, catabolism, DNA transcription, and Tcell functions.

Likewise, he majority of the biological processes found for the decreased expression group (Figure 4C) were related to metabolism and DNA replication, however, other notable percentage distributions were exclusively identified for hematopoiesis and myelopoiesis.

Another more focused method of analysis was performed on the gene expression data of the three major temporal groups (increased, variable, and decreased expression).Significant (p,0.05)known biological pathways determined by KEGG and GenMAPP databases were identified and are presented in Table 1.This analysis revealed that certain significant pathways identified were unique to the gene expressions comprising only one of the three major temporal groups.Specifically, the increased expression group exclusively contains genes involved in pathways for matrix metalloproteinases (MMPs), Wnt signaling, axon guidance, TGF-b signaling, and inflammation.Contrastingly, only the variable expression group contains genes associated with glycolysis and gluconeogenesis, notch Signaling Pathway, natural killer cell mediated cytotoxicity, and the B cell receptor signaling pathway, among others.Furthermore, only the decreased expression group consists of genes involved in heme biosynthesis, the p53 signaling pathway, and the hematopoietic cell lineage.

Genes associated with significant pathways of interest from the three major temporal expression groups were hierarchically

lustered based on their expression profiles and are presented as heat maps (Figure 5) showing fold-change
values (vs.day 0 time point) of log base 2 expression values.Pathways associated with the increased expression group that are displayed as heat maps (Figure 5A) incl

e matrix me
alloproteinases (MMPs) (MMP), Wnt signaling pathway, TGF-b signaling pathway, and inflammatory pathway.Other pathways associated with the variable or decreased expression group that are displayed as heat maps are notch signaling pathway (Figure 5B) and hematopoietic cell lineage pathway (Figure 5C), respectively.


Identification of significant biological pathways associated with each time point post-ablation of marrow

As opposed to significant pathways identified to be associated with the three major temporal groups based on clustering of all significant genes, significant (p,0.05)biological pathways were also identified by analyzing the expression of significant genes at each time point following bone marrow ablation.Pathway results for significant gene expressions on each day vs. day 0 are shown in Table 2, Table 3, and Table 4.All days except Day 56 exhibit the general pathways of ECM-receptor interaction and focal adhesion as significant in their resulting lists, with the majority of time points exhibiting other general pathways related to cell communication or cell cycle.Collectively, the inflammatory response and blood clotting cascade pathways are shown to be significant at all time points.Other pathways of interest, such as those shown in Figure 5 are significantly expressed only at certain time points.For example, Wnt signaling and matrix metalloproteinases (MMPs) were identified to be differential pathways involved on Days 1 to Day 7, and Day 1 to Day 10, respectively, post-ablation of bone marrow.Interestingly, axon guidance and the p53 signaling pathway are found to be significantly expressed exclusively on Day 10, and the VEGF signaling pathway is identified as significant uniquely on Day 14.


Identification of significant transcription factors expressed at each time point post-ablation of marrow

The XY-scatter plot graph (Figure 6) displays the percentage distribution of significant transcription factors expressed at each time point, and Table S3 provides the list of Affymetrix GeneChipH Rat Genome 230 2.0 Array probe set ID and transcription factor names for each time point.The complete list of probe sets for transcription factors present on the Affymetrix GeneChip H Rat Genome 230 2.0 Array was found to be 1,254.Comparing the significant probe set lists for each time point postablation (each day vs. day 0) with the master list of known transcription factors, the number of transcription factors (TFs) expressed at each time point was determined to be: 42


Discussion

This study was the first to identify significant biological processes and pathways from genome-wide transcriptional analysis during all phases of intramembranous bone regeneration following marrow ablation in a rat model.We used an oligonucelotide microarray technology, specifically fabricated by Affymetrix for Rat Genome 230 2.0, and used univariate F-test to determine a list of 9,062 gene probe sets significantly differentially expressed over all time points, which is roughly 30% of the 31,042 gene probe sets present on the total microarray.The time points of days 1, 3, 5, 7, 10, 14, 28, and 56 days post-ablation of marrow were selected based on knowledge of the complete timeline of all general bone healing processes involving phases of inflammation, repair and new bone formation, and remodeling.This timeline is confirmed by the data presented, implicating that most expression profiles return to baseline at 56 days, as shown by the average expression profiles (Figure 3).A Bayesian model-based clustering method Figure 4. Percentage distribution of significant Gene Ontology/Biological Processes associated with each major group of clustered gene expression profiles.DAVID EASE (version 2.0) software was used to analyze the gene lists for each major temporal group (increased, variable, and decreased expression) to find all categories associated with the Biological Processes domain of the Gene Ontology (GO) consortium, using a significance threshold of p,0.05.Individual gene ontology results with related or overlapping functions were manually combined into overall major categories.Each biological process categ ry is assigned a distinct color, and consistent colors were used for comparisons of the same biological process category found for more than one of the three major groups of clustered gene expression profiles/each of the three pie graphs, which include Increased Expression Group shown in Figure 4A, Variable Expression Group shown in Figure 4B, and Decreased Expression Group shown in Figure 4C.For example, the major category of ''metabolism'' is consistently assigned the color red.The number of genes within each of the three major temporal groups (increased, variable, and decreased expression) th t underwent gene ontology/biological process analysis is denoted in the title of each pie graph.doi:10.1371/journal.pone.0012987.g004was used to determine major temporal expression groups based on gene expression profiles.The percentage distribution of significant biological processes (Figure 4 and Table S2A-D) and the differential biological/signaling pathways (Table 1) identified for increased, decreased, and variable expression groups serve as a genome-wide overview of the main biological functions that occur during marrow ablation-induced intramembranous bone regeneration in the rat femur.Additionally, identification of biological pathways (Table 2, Table 3, and Table 4) and transcription factors (Figure 6 and Table S3) significantly expressed at each time point post-ablation provides important insight into the regulation of intramembranous bone regeneration in a complete temporal context.Each of the phases of new bone formation identified in the histologic analysis generally progress in an ''outside-in'' fashion toward the core of the initial clot matrix.Thus, even though distinct days can be assigned to a given phase, because of the spatial and temporal heterogeneity inherent in the regenerative process, almost all stages can be identif ed in some, albeit limited, regions of the ablation site in samples recovered between days 7 and 10.This heterogeneity implies that peaks of activity occur concurrently in different locations and that critical events and phenotypic transformations reflected in gene expression at these transformative sites may occupy small footprints in the overall expression profiles.Alignment of gene expression patterns with the phases described above (see Results section) in their spatial context may eventually help identify genes that are critical to bone regeneration but appear at low levels within the more dominant spatial domains at any particular time.

Previous studies of gene expression using marrow ablation models selectively investigated individual or very limited number of osteogenesis-associated cytokines or growth factors which are described and summarized in previous publications by our group [2,19].Results determined by the present study complement conclusions highlighted in those previous publications, in that intramembranous bone regeneration generally relies upon the control of gene expression for inflammatory cytokines, osteogenic growth factors, transcription factors, matrix proteins and proteinases, and bone remodeling markers.

Focusing on the most differentially expressed genes in Clusters 8, 1, and 7 in Figure 3A (Increased Expression) and Table S1, se eral osteoblastic-and osteogenesis-associated genes of interest [2] are identified including osteopontin (SPP1 in Cluster 8 ).It is known that periostin has a functional role as a matricellular protein in a variety of tissue remodeling and wound repair situations [32], and it has been recently confirmed that periostin is a novel marker of intramembranous ossification, but not endochondral ossification, and exhibits upregulated expression in the fibrous component of fibrous dysplasia, due to its prominent influence on collagen fibrillogenesis [33].Many of the aforementioned osteoblastic-associated genes exhibit peak upregulation at day 5 or Day 10, which subsequently and, respectively, result in bone formation observed in the histology tissue by Days 10 and 14 (Figure 2).

Clusters 8, 1, and 7, in Figure 3A (Increased Expression) and Table S1, also include a number of osteoclastic-and bone matrix remodeling-associated genes of interest [2] ).The aforementioned osteoclastic-associated genes exhibit peak upregulation at day 10, which subsequently result in bone resorption and bone remodeling observed histologically in the tissue by Day 14 and 28 (Figure 2).

No evidence of chondrogenesis was observed histologically in the tissue at any time point, and two major chondrogenic marker genes collagen type II (Col2) and transcription factor Sox9 were not identified as significantly expressed in our model.Therefore, true chondrogenesis was not identified in this model.However, other chondrocytic-related marker genes were found within the significantly expressed genes in Table S1 and in Figure 3A and ), where they exhibit peak upregulation at day 5 followed by a rather steep decline to return to baseline by day 14.This is suggestive of a transient existence of chondrogenesis.This observation is similar to the one reported by Nah et al [34] on the ''covert chondrogenic phase'' in their study of intramembranous bone formation in chick embryonic frontal bone and sternum.

Adipocyte marker genes that are also found to be significantly expressed (Table S1, Figure 3B  Interestingly, all the expression profiles of Clusters 2, 4, and 5 which contain adipocytic markers display a downregulated expression profile during day 5 until day 10, which contrasts the expression profiles of Clusters 8, 1, and 7 which contain osteoblastic-associated marker genes and which show peak upregulation during day 5 to day 10.In other words, adipocyticassociated genes are downregulated while osteoblastic-associated genes are upregulated.Adipocytic-assocatied genes in Clusters 4 and 5 return to baseline by Day 28, leading to restoration of the fatty marrow, while adipocyte markers genes contained within Cluster 2 become upregulated at Day 14 and remain so until day 56.

One of the major pathways identified in our model to be significantly expressed in the increased expression group (Table 1), was the molecular function network of matrix metalloproteinases (MMPs), which was specifically determined to be significantly expressed at 1, 3, 5, 7, and 10 days post-ablation (Table 2 and Table 3).Eight MMP-specific genes with significant increased expression that were identified using the GenMAPP database were Bsg, Timp 1, Timp 2, MMP2 (gelatinase A or type IV collagenase), MMP12 (metalloelastase), MMP13 (collagenase), MMP14 (membrane-type 1 MMP), and MMP16 (membrane-type 3 MMP), and these are highlighted in the heat map for MMP genes in Figure 5A.It can be seen that the expression values of MMP2, MMP14, TIMP2, and MMP16 are hierarchically clustered ogether.This result is expected since MMP14 (also known as MT1-MMP) and TIMP-2 are specifically involved in the cell surface mechanism to activate proMMP-2 [35,36], that essentially contributes to the promoting the invasive capacity of tumor cells [35] and mesenchymal stem cells [37], and it is also known that MMP16 (MT3-MMP) activates MMP2 by cleavage of proMMP-2 [38].From the results of the temporal expression profile clustering, the most up-regulated genes were TIMP1, MMP2, and MMP12 (Cluster 8 in Figure 3A and Table S1).Other MMPs that exhibited significant increased expression were TIMP2, MMP13, MMP14, and MMP16 (Cluster 1 in Figure 3A and Table S1).MMP13, a known collagenase, exhibited sustained strong induction in our model for up to 14 days postablation, and is known to have important functions in bone formation and remodeling.It has previously been reported that MMP13 is specifically required to cleave collagen type II and aggrecan during the transition from cartilage to bone at the growth plate during long bone development, in addition to mediating initial and continual remodeling of trabeculae during ossification [39,40].This knowledge correlates with the relationship of MMP13 and Aggrecan gene expression observed in the data from our model.Aggrecan was found as a significantly upregulated gene in our model (Cluster 1 in Figure 3A and Table S1) with its peak expression at the day 5 and 7 time points.MMP13 reached peak expression on day 10 and 14, during which a significant decline towards baseline was seen for Aggrecan expression.It has also been recently reported, using non-stabilized or stabilized fracture models in mice, that MMP13 is required for both endochondral and intramembranous ossification during bone repair, likely for initial degradation of ECM prior to the invasion of blood vessels and osteoclasts [41].These interesting findings implicate that some extent of transient chondrogenesis is occurring in our model in addition to the predominant intramembranous bone regeneration induced by bone marrow ablation, even though cartilage were not histologically evident.Furthermore, our previous study demonstrated that a biomimetic polymer hydrogel with a semi-interpenetrating network that incorporates an MMP13 degradable crosslinker peptide and an integrin-binding peptide (containing the RGD domain) implanted in a rat femoral ablation model leads to significantly enhanced bone regeneration, likely due to enhanced osteoblast migration and proliferation in such constructs [18].Further studies involving tailored design of biomaterials could benefit from the knowledge presented here regarding MMP-specific genes to target particular MMP-related peptide sequences for certain types of tissue regeneration.

Another key pathway identified to be significantly expressed in the increased expression group (Table 1) was the cellular process pathway of Wnt signaling, which was found to be significantly expressed at the study timepoints of 1, 3, 5, and 7 days postablation.A study by Kim et al. used a model where bone regeneration in a 1.0-mm monocortical tibial defect occurs exclusively through intramembranous ossification, and they reported that Wnt signaling is upregulated at the injury site, thereby prompting marrow-derived osteoprogenitor cells to proliferate and mediating subsequent osteoblast differentiation [42].In our model, Wnt signaling related genes for Wnt5a  S1).In the previously mentioned study by Kim et al., Wnt5a was also identified as significantly expressed in all domains of the injury site in their primarily intramembranous ossification model, which included the osteocytes, endosteum, bone marrow, growth plate, and the periosteum, which was the only domain adjacent to Figure 5. Clustered heat maps of genes involved in significant pathways identified for each major group of clustered gene expression profiles.Cluster 3.0 was used to cluster the log base 2 fold-change expression values (vs.day 0 time point) of genes associated with certain significant biological pathways of interest from the three major temporal expression groupings (increased, decreased, and variable expression).Pathways of interest from the Increased Expression Group include those for matrix metalloproteinases (MMPs) (GenMAPP), Wnt signaling (GenMAPP), TGF-b signaling (KEGG and GenMAPP), and inflammatory response pathway (GenMAPP) shown in Figure 5A.A pathway of interest from the Variable Expression Group is notch Signaling Pathway (KEGG) shown in Figure 5B, and a pathway of interest from the Decreased Expression Group is hematopoietic cell lineage (KEGG) shown in Figure 5C.The log base 2 fold-change ratios were clustered using hierarchical clustering with a centered correlation distance/similarity metric and average linkage clustering method.The clustered data table file was viewed in TreeView using the pixel setting contrast default of 3 and using blue and red to represent positive and negative fold-change expression values, respectively.The resulting clustered heat maps display the fold-change (vs.day 0 time-point) expression values for each gene in the rows, and for each time point in the columns, and the order of the rows is based on the result of clustering genes of similar profiles.Each heat map displays the total list of genes on the Affymetrix GeneChipH Rat Genome 230 2.0 Array that are associated with the particular pathway of interes .Genes from the initial significant gene list of 9,062 (and noted as ''Gene List Hits'' in Table 1) associated with each pathway are highlighted with a bold green font for their gene symbols on the right side of the heat maps.doi:10.1371/journal.pone.0012987.g005

Table 2. Significant biological pathways identified at day 1, 3, and 5 (vs.day 0) post-ablation. the injury site where there was an endochondral ossification response [42].

The present results with marrow ablation-induced intramembranous bone regeneration share similarities with two other microarray studies using rat closed fracture models, which also identified genes for Wnt5a, Frizzled 2, and b-catenin [1,43].In the study by Hadjiargyrou et al., the identified temporal expression for Wnt5a was interestingly reported to increase at day 3 postfracture, then decline at days 5-7, then return to an increased expression at 10-14 days, and finally drop back down to baseline levels by 21 days [1].A similar but more early and transient temporal expression profile for Wnt5a was found in our model, with an increase at day 1 post-ablation, then a decrease on day 3, then a marked increase at day 5, then another decrease at days 7 and 10, followed by a further decrease towards baseline by 14 days.Wnt5a is known as a non-canonical Wnt ligand that signals independent of b-catenin through the Wnt/Ca 2+ pathway, which not only regulates Ca 2+ flux and Ca 2+ -sensitive protein kinases and transcription factors, but can inh bit the canonical Wnt pathway by promoting degradation of b-catenin [44].Specific mechanisms involved with the regulation of bone regeneration by noncanonical Wnt pathway are still relatively unknown, but several studies have shown evidence supporting the role of Wnt5a in the early stages of fracture repair, namely inflammation and chondrogenesis [1,44].Additionally, there is other evidence that Wnt5a induces osteoblastogenesis by reducing PPARc-induced adipogenesis in bone marrow mesenchymal stem cells [45].Considering this knowledge of Wnt5A, in addition to the information previously discussed regarding MMP13, it is possible that early and transient chondrogenesis may be occurring in our predominantly intramembranous regeneration model induced by marrow ablation.Altogether, these intriguing findings from our study and many other reports implicate that Wnt signaling in bone repair models, fracture or marrow ablation, is highly intricate and currently not well understood, evidently involving the activation of both canonical and non-canonical Wnt signaling pathways [44,46].

It is generally known during the early bone repair phase that mesenchymal stem cells (MSCs) migrate, proliferate, and undergo differentiation towards the osteoblastic lineage.It is of interest to further investigate signaling pathways and genes expressed by MSCs during this phase for overall enhancement and acceleration of bone regeneration following injury and surgery, and the temporal transcriptional profiling data results presented here are expected to aid in that research.Signaling pathways known to be in involved in the self-renewal and osteogenic differentiation of MSCs are Wnt signaling, BMP/TGF-b signaling, and notch signaling [47], which are pathways identified in this study to be significantly expressed within the increased or variable expression groups, respectively.There is evidence that MSCs express Wnt2, Wnt4, Wnt5a, Wnt11, Wnt16, Fzd2, Fzd3, Fzd4, Fzd5, Fzd6, and Dkk1 [48], and that canonical Wnt3a increases levels of b-catenin and the proliferation rate while noncanonical Wnt5a impedes the process of chondrogenesis [47].Indeed, there is cross-talk between canonical and noncanonical Wnt signaling in MSCs, as canonical Wnt3a suppresses Wnt5a to maintain MSC in an undifferentiated and self-renewing state, while noncanonical Wnt5a inhibits Wnt3a in order to mediate enhancement of osteogenic differentiation of MSCs [49,50].Specifically, Wnt5a suppresses b-catenin/TCF signaling to decrease the level of cyclin D1 and proliferation rate of MSCs [49,51].

The notch signaling pathway is thought to be involved in osteoblastogenesis and skeletogenesis but its particular role is not well understood [47].It has been reported that within the bone marrow, notch signaling suppresses osteoblastic differentiation and the Wnt/b-catenin signaling pathway generally acts to maintain a pool of proliferating mesenchymal progenitors [52,53].Although not specifically identified as significant by the KEGG database through DAVID EASE, genes involved in the notch Signaling Pathway and differentially expressed over all time points in our model, include Hes1 [29577], Jag1 [29146], and Numbl [292732].It was recently confirmed that Hes1 plays an important role in mediating the inhibition of osteoblastogenesis and the Wnt/bcatenin pathway by intracellular domains of notch signaling, likely by preventing the interaction of b-catenin with the transcriptional c -repressor Groucho/TLE and LEF-1 [53].The notch ligand gene Jag1 is known as an evolutionarily conserved target of the canonical Wnt signaling pathway and is a key molecule for induction of self-renewal and maintenance of homeostasis of stem and progenitor cells [54], while Numb and Numbl induce differentiation by inhibiting notch signaling in progenitor cells [55].

Regarding the TGF-b signaling pathway, many significantly expressed genes were identified by the KEGG or GenMAPP databases and are shown as three major clusters in the TGF-b signaling pathway heat map in Figure 5A.One notable cluster   [11,[57][58][59].Interestingly, in our data Thbs2, which is clustered with TGF-b2, is highly expressed from day 3 to day 10 post-ablation, and has been described as an autocrine inhibitor of proliferation secreted by MSCs [60] and an important regulator, along with osteonectin (also known as Sparc), of the osteoblast lineage and bone remodeling [61].There are a few recent reports regarding genome-wide transcriptional profiling of mouse and human MSCs and the identification of differentially expressed genes and pathways [62,63].Specifically it was reported that transcripts uniquely expressed by murine bone marrow-derived MSCs are mostly transcription factors and genes downstream of the Wnt signaling pathway and also MSCs exhibit a unique therapeutic effect on T cells by preventing their proliferation and supporting their survival [62].Another study performed microarray-based genome-wide differential gene expression analysis to identify the 20 most significantly expressed genes in human bone marrowderived MSCs, including several we have identified as differentially expressed in our model of intramembranous bone regeneration such as Periostin (Postn [361945]), Col3a1 [84032], Col6a3  S1) [63].

There are many noteworthy similarities and differences of the genome-wide transcriptional analysis results presented here for intramembranous bone regeneration induce by a marrow ablation model with results previously reported by Bais et al. regarding the endochondral bone formation process in fracture healing [20].In their paper, it was reported that a little over onehalf of the genes expressed in the mouse genome were differentially regulated during fracture healing [20].Our data indicates that roughly one-third of the genes expressed on the Affymetrix GeneChipH Rat Genome 230 2.0 Array were differentially expressed during intramembranous bone regeneration and healing following marrow ablation.Interestingly, many commonalities can be seen when comparing gene ontology/ biological process analysis of the three major temporal expression groups (up, variable, and down) from Bais et al. with the results presented here.Specifically, the many biological process that share similar percentage distributions for the up or increased expression group identified in both studies are cell adhesion, cell cycle, cytoskeleton/actin, development, immune response, ion transport, many aspects of metabolism, motility, neurogenesis, a variety of cellular signaling including Wnt (see Table S2A-D from this study and Table S1 from Bais et al.), and vasculo enesis [20].Marked differences include that a greater percentage distribution (7.7%) for skeletogenesis is found from our increased expression group data for intramembranous bone regeneration model compared to a much lower value (1%) for skeletogenesis shown in the increased expression group in Bais et al. which uses a fracture healing model [20].Other similarities and differences can be observed by comparing and contrasting the biological process percentage distribution for the variable (Figure 4B) and decreased (Figure 4C) expression groups from this study with the variable and down groups in Bais et al. [20].

More similarities exist with our results of significant biological pathways identified (from KEGG and GenMAPP databases) for each of the three major temporal groups (Table 1) when compared to pathways identified (by KEGG database) for the three major temporal groups in Bais et al. [20].The increased expression groups in both studies showed significant association with ECM-receptor interaction, focal adhesion, axon guidance, focal adhesion, TGF-b signaling, cell communication, adherens junction, basal cell carcinoma, and Wnt signaling.The same pathways identified for variable group and the variable group in Bais et al. include B cell receptor signaling pathway and natural killer cell mediated cytotoxicity, and the same pathways identified for the decreased expression groups in both studies include cell cycle (negative regulators), hematopoietic cell lineage, ABC transporters (general), and porphyrin and chlorophyll metabolism [20].Comparing these results from our intramembranous bone regeneration model and previously reported results from a fracture healing model implicate many similarities regarding significant biological processes and pathways during bone repair.

Given the result from this study that roughly 30% of the rat genome was significantly expressed during all time points of marrow ablation-induced intramembranous bone regeneration, it is of interest to discuss what fraction of the genome is involved in other forms of wound and tissue repair.Similar to bone, the liver has the ability to regenerate and repair with no resulting scar tissue.A genome-wide expression study using a rat liver regeneration model, in which a partial hepatectomy was performed recently reported that approximately 5.4% of the rat genome was differentially regulated during rat liver regeneration [64].Similar to the regenerati

of bone, liver regener
tion involves a coordinated cascades of biological events, although our data confirms that intramembranous bone regeneration involves a much greater number of significantly activated/expressed genes compared to liver regeneration.A very recent publication using a human genome-wide microarray analysis with a cutaneous wound healing model of thermally injured skin determined that approximately 4.4% of the human genome was significantly expressed during a 18-day healing period [65].It is known that unlike bone repair, skin repair results in scarring.Interestingly in their study, osteopontin (also known as SPP1), was highly upregulated in the first 17 days and was significantly expressed in our intramembranous bone regeneration model (see Cluster 8 in Figure 3 and Table S1 and the heatmap for TGF-beta signaling genes in Figure 5A).The significant expression of osteopontin in wound healing supports the knowledge that it is has a broad yet critical role in injury site extracellul r microenvironments containing highly proliferative cells and u dergoing rapid remodeling.Another study using a rat genome-wide microarray analysis with an injured spinal cord model found that during a 90 day post-injury period, approximately 15% of the rat genome was significantly expressed [66].It was also noted that the overall expression profiles of significantly increased tissue repair genes and the timing and sequence of post-spinal cord injury gene expression resemble the phas