Transcriptional Profiling of Bone Regeneration

The healing of skeletal fractures is essentially a replay of bone development, involving the closely regulated, interdependent processes of chondrogenesis and osteogenesis. Using a rat femur model of bone healing to determine the degree of transcriptional complexity of these processes, suppressive subtractive hybridization (SSH) was performed between RNA isolated from intact bone to that of callus from post-fracture (PF) days 3, 5, 7, and 10 as a means of identifying up-regulated genes in the regenerative process. Analysis of 3,635 cDNA clones revealed 588 known genes (65.8%, 2392 clones) and 821 expressed sequence tags (ESTs) (31%, 1,127). The remaining 116 cDNAs (3.2%) yielded no homology and presumably represent novel genes. Microarrays were then constructed to confirm induction of expression and determine the temporal profile of all isolated cDNAs during fracture healing. These experiments confirmed that ∼90 and ∼80% of the subtracted known genes and ESTs are up-regulated (≥2.5-fold) during the repair process, respectively. Clustering analysis revealed subsets of genes, both known and unknown, that exhibited distinct expression patterns over 21 days (PF), indicating distinct roles in the healing process. Additionally, this transcriptional profiling of bone repair revealed a host of activated signaling molecules and even pathways (i.e. Wnt). In summary, the data demonstrate, for the fist time, that the healing process is exceedingly complex, involves thousands of activated genes, and indicates that groups of genes rather than individual molecules should be considered if the regeneration of bone is to be accelerated exogenously.

potential of both tissue and genetic engineering, it is anticipated that exogenous acceleration of fracture healing could increase the overall numbers of fractures that heal successfully, as well as reduce the number of patient days lost due to incapacity. To achieve this goal, however, it is essential that we expand our understanding of the interdependent stages of fracture healing (inflammation, chondrogenesis, ossification, remodeling), and ideally, how they contribute to a biomechanically functional bone.
A number of biochemical and biophysical interventions have been devised with the goal of accelerating the healing of skeletal fractures. The biological strategies include biomimetically inspired skeletal graft substitutes (hydroxyapatite, calcium carbonate), purified or recombinant molecules with chondrogenic and osteogenic attributes (i.e. growth factors), gene therapy, and stem cell reservoirs introduced by biodegradable matrices (2). The biophysical arsenal includes low intensity ultrasound (3,4), mechanical stimuli (5), and electromagnetic fields (6). While the basic science foundation supporting these modalities is strong, the clinical results have been inconclusive. Therefore, it becomes reasonable to conclude that the complexity of the healing process is being underestimated and that the healing process cannot ultimately be determined by a singular idealized molecule, material, or stimulus.
The great majority of studies that have examined the molecular basis of healing have focused on the expression of specific genes, with the bulk of these studies concentrating on the regulatory role of known extracellular matrix (ECM) 1 genes (7)(8)(9)(10)(11) and growth factor genes/proteins (12,13). The range of genes evaluated has recently expanded to include other protein families, including intracellular signaling molecules (14), transcription factors (15,16), cytokines (17,18), adhesion molecules (19), and enzymes (20,21). While this extensive body of excellent work has helped demonstrate the temporal and spatial roles of specific genes, it has also served to indicate that we have limited knowledge of how extensive the transcriptional control of the repair process may be or identify which specific processes are the critical regulators of successful bone healing. With the advent of more sophisticated techniques in gene expression analysis (i.e. microarrays), it becomes possible to simultaneously examine the transcriptional activity and interaction of thousands of genes.
Given the biological complexity of fracture healing, a process morphologically characterized by inflammation, chondrogenesis, and osteogenesis, it is reasonable to hypothesize that it is regulated by a very large number of transcriptional events (22,23). This hypothesis is supported by the marked similarities between the repair process of bone and embryonic development of the skeleton, marked by key cellular events (migration, adhesion, proliferation, and differentiation), all of which require the tightly orchestrated activity of thousands of proteins whose expression patterns rely on both extracellular and intracellular signals.
The work reported here combines SSH (24) and custom cDNA microarrays (25) to elucidate the transcriptional complexity of bone regeneration, as modeled in the rat femur (26). The data demonstrate that a very large number of genes are transcriptionally regulated during the repair process over the distinct, but interdependent, stages of healing. A large number of expressed sequence tags (ESTs) are also identified, which are then clustered with respect to functionally known gene families. Furthermore, we identify the Wnt pathway as an example of a signaling pathway involved in the repair process. Taken together, these results constitute the first demonstration of the transcriptional complexity underlying the early phases of the healing mammalian fracture callus and may help to identify critical components, as well as novel candidate genes, essential for successful bone regeneration.

EXPERIMENTAL PROCEDURES
Fracture Model-All methods and animal procedures were reviewed and approved by the University's Laboratory Animal Users Committee and met or exceeded all federal guidelines for the humane use of animals in research. The rat femur fracture model was described previously (26) and used extensively in our studies (22,23,27). A set of four animals was euthanized each at 3, 5, 7, 10, and 14 days post-fracture and a set of three animals at 21 days. Following euthanasia by CO 2 inhalation, the contralateral control femur from each animal, as well as the fracture calluses, were dissected free and processed for RNA extraction.
RNA Purification-Total RNA was isolated from each individual fracture callus, as well as from each intact bone, which included bone marrow and articular and normal growth plate cartilage, using the ToTALLY RNA kit (Ambion) based on the method of Chomczynski and Sacchi (28) and as described previously (22,23,27).
Suppressive Subtractive Hybridization-The samples used for these analyses were derived by pooling the RNA from two intact femurs and comparing it with RNA pooled from the fracture callus tissues of one animal harvested at each of PF days 3, 5, 7, and 10. Complimentary DNAs were synthesized from 1 g total RNA using the SMART PCR cDNA synthesis kit (CLONTECH), and subtractive hybridization was performed with a PCR-select cDNA subtraction kit (CLONTECH). cDNAs derived from the fracture callus material, considered the "tester" pool, and cDNAs from intact femurs, considered the "driver" pool, were digested in RsaI (New England Biolabs). To select transcripts up-regulated by the fracture repair process, PCR adaptors were ligated to the tester pool population (fracture callus). The tester cDNA pool was then hybridized with excess cDNAs (ϳ15-fold) from the driver pool (control bones). After hybridization, suppression PCR using primers specific for the tester PCR adaptors selectively amplified differentially expressed transcripts. Amplified cDNA sequences were then ligated into the T/A cloning vector pT-Adv (CLONTECH). Approximately, 4,500 cDNA clones were sequenced using an ABI 3700 DNA sequencer (Applied Biosystems). Finally, all sequences were checked for homology using the BLAST algorithm found at www.ncbi.nlm.nih.gov/blast/. For blastn, ϳ65% of the sequences had E values of 2.0 ϫ 10 Ϫ19 or better (parameters: BLOSUM62; word size 12). For blastp, ϳ40% of the sequences had E values of 1.0 ϫ 10 Ϫ6 or better (parameters: BLOSUM62; default NCBI Gap Costs).
cDNA Array Construction-Differential screening of the fracture callus-induced transcript library was performed through a series of steps. First, individual colonies were grown overnight at 37°C in 500 l of Luria broth medium containing 50 g/ml ampicillin. The next morning, 2 l of the culture was transferred to individual wells of a 96-well PCR plate containing 98 l of PCR master mix (10 l of 10ϫ PCR buffer, 2 l of 10 M pT-Adv nested primers 1 and 2, 1 l of 20 M dNTPs, 82 l of sterile H 2 O, and 1 l of Taq DNA polymerase (PerkinElmer Life Sciences). Using the GeneAmp PCR System 9600 (PerkinElmer Life Sciences), PCR amplification was performed by cycling for 2 min at 94°C followed by 35 cycles consisting of 40 s at 95°C and 3 min at 68°C. PCR amplification products were analyzed by electrophoresis on a 1% agarose gel, and ϳ90% of the clones were found to contain inserts, with an average insert size of 250 -500 bp.
To prepare nylon arrays of the PCR-positive clones, free primers were removed from the PCR reactions by filtration and the samples concentrated to ϳ50 ng/l. Each sample was then spotted on the nylon filters (0.25 l/spot). In addition to the 3,635 subtracted cDNA clones (representing all 588 known genes, the 821 ESTs and 116 novel sequences), 257 control samples were also spotted (mitochondrial DNA, ribosomal RNA, genomic DNA, plasmid, actin, tubulin, GAPDH, yeast DNA, buffer, bacterial DNA (Escherichia coli/Bacillus subtilis), as well as mammalian DNA (mouse, rat, human, monkey)). Membranes were denatured by soaking for 10 min in a solution of 0.5 mM NaOH, 1.5 mM NaCl, and neutralized by successive 5-min washes in 0.5 mM Tris-HCl (pH 8.0), 1.5 mM NaCl, and 2ϫ sodium chloride-sodium citrate (SSC) buffer. The membranes were then UV cross-linked using a Stratalinker (Stratagene).
Target Labeling and Microarray Hybridization-RNA derived from the calluses of three animals was pooled for each separate time point (PF days 3, 5, 7, 10, 14, and 21). The RNA sample representing intact bone was established by pooling specimens from intact femurs of three different animals. The protocol used for labeling target, hybridization, and washings was identical to that for GeneFilters from Research Genetics, with the exception that the last wash was carried out in a stringent wash in 0.2ϫ SSC, 0.1% SDS at 65°C for 30 min. Subsequently, all the blots were simultaneously exposed to a PhosphorImager screen (Amersham Biosciences) for 3 days prior to capturing the final image using a PhosphorImager (Amersham Biosciences).
Image Analysis-Each membrane image was analyzed using the GenePix Pro (version 3.0) microarray software package. Measurements of the optical intensity of the 244 pixels contained within each spot generated a median intensity value. The average background intensity was then calculated from the 1,060 blank spots scattered throughout each membrane. This average background intensity was then subtracted from the median intensity of each spot and normalized to the 18 S ribosomal RNA-positive control spots. Finally, replicate spots for each gene were averaged prior to further analysis.
Cluster Analysis-Clustering of the up-regulated genes, as a function of time, was performed using average linkage analysis following data normalization and filtering (includes background subtraction and normalization as stated above (Image analysis)). In addition, a block effect correction was also performed. We defined blocks in two directions. First, 38 rank blocks were generated based on the mean intensity value for each gene, and second, four (2 ϫ 2) adjacent blocks were combined. Each of these block effects was calculated using a generalized linear model and corrected from the intensity which led to the residual value from the generalized linear model fitting. Finally, the median of each rank block for each array was added to this residual. Since a single representative value was needed for each gene, and since some genes were represented by multiple cDNA clones (more than one spot on the membrane), some spots that generated inconsistent expression patterns as compared with others (of the same gene) were considered as outliers and were filtered out. Based on correlation (both Pearson's and Spearman's), 403 spots were excluded due to their inconsistent expression patterns (as described previously), and another 29 spots were excluded because of their low expression level (negative intensity values). With this normalized intensity value, we performed pairwise average linkage cluster analysis to establish clusters by grouping genes which shared the same/similar temporal expression patterns. Pearson's correlation coefficient, which statistically captures similarity in "shape" was used as similarity measure. The final number of clusters was determined to be 16 using pseudo F and pseudo t 2 statistics (29,30). To show that the final selection and number of clusters were reasonable, intracluster correlations of each of the 16 clusters, as well as the results from discriminant analysis, were examined. This data analysis was performed in logarithm base 2 space after standardization.

RESULTS
Transcriptional Activity in the Fracture Callus-To establish an in depth understanding of the transcriptional activity occurring during the early stages of a healing rat femur fracture, RNA was examined from four different time points (PF days 3, 5, 7, and 10) and compared with that of intact bone (containing cartilage and bone marrow). These initial time points were selected to represent specific physiological events of the healing callus, including inflammation, chondrogenesis, and ossifica-tion (Fig. 1). The spatial complexity and structural interdependence of these early events of the mammalian callus have been described previously (7,8,19,22,23).
Of the 4,500 cDNAs derived from SSH, ϳ3,635 were successfully sequenced (contained a readable insert). After BLAST searches, of 3,635 clones, 65.8% had homology to 588 known genes (represented as 382 singletons), 31% had homology to 821 ESTs, and the remaining 3.2% (116) had no homology match and presumably represent completely novel genes. The known genes reflected a variety of families with diverse functions in cell cycle regulation, cell matrix and cell adhesion, ECM construction, inflammation, general metabolism, signaling, transcriptional regulation, protein transport, etc. (Supplemental Table 1). The abundance of each gene within the library (represented as # of clones) is also indicated in this table.
Temporal Gene Expression Analysis via Custom-made cDNA Arrays-Given the large number of cDNA clones present in the subtracted library, it was deemed practical to analyze the expression of all cDNA clones simultaneously through the use of custom microarrays. These microarrays included the complete subtracted library (3,635 cDNAs), as well as 257 control spots (see "Experimental Procedures"). Initially, experiments utilizing the identical RNA samples (from single animal per time point) used in SSH were utilized to confirm the hybridization results, as well as the production and integrity of the microarrays. Results from these experiments yielded reproducible data between each run (data not shown). Once confirmed, RNA samples from multiple animals at each specific PF time point were pooled (n ϭ 3), and hybridizations of arrays were repeated to show how transcriptional activity was modulated as a function of time in the healing process as compared with RNA samples from intact bone (n ϭ 3). Hybridization was performed using a total of seven identical membranes, one for each of the six PF time points (days 3, 5, 7, 10, 14, and 21), as well as one for intact bone (control). We found that ϳ90 and ϳ80% of the subtracted known genes and ESTs are up-regulated (Ն2.5-fold) during the repair process (at any given PF day), respectively. The remaining cDNA clones probably represent false positives, most likely resulting from the SSH. Representative membranes following hybridization of RNA isolated from intact bone and PF day 10 callus are shown in Fig. 2 and clearly demonstrate differential gene expression.
Data resulting from imaging the seven filters was used to generate scatter plots, representing the ratio for each gene between each of the PF day callus (y axis) over that of intact bone (x axis, Fig. 3). The figure clearly shows the increased gene expression (shift of dots to the y axis) during the progression of the fracture callus through its early stages (PF days 3-10), thus confirming the initial subtraction (designed to identify up-regulated genes). Relative to intact bone, the highest level of expression is observed at PF day 14 (note the higher intensity and shift of dots), and by PF day 21, expression levels begin to decline toward the control (represented by diagonal line, Fig. 3). Based on these data, the exact temporal expression levels of all known genes are shown in Supplemental Table 2 and demonstrate that the expression patterns change extensively throughout the healing process and that these shifts are distinctly different from gene to gene. Finally, the fold change in expression observed in this study for many of these genes (i.e. collagens, osteopontin, osteonectin, etc.) are consistent with those previously determined using Northern analyses by our laboratory (22,23), as well as others (8 -10, 15).
Clustering Analysis of ESTs and Novel cDNAs-In an attempt to group known genes, ESTs and novel cDNAs based on their co-regulated expression patterns, cluster analysis was performed using data derived from the microarray experiments. Statistical analysis generated 16 distinct clusters ( Fig.  4 and Supplemental Table 3). A common feature of each cluster is the trend of the expression pattern, which invariably shifts upwards from intact to PF, indicating the expected up-regulated expression of genes relative to intact bone. Furthermore, the clusters show distinct patterns of expression, ranging from gene activity which peaked early and then declined (Fig. 4A,  Clusters 1, 10, 13, and 16), to clusters that continue to rise through the healing process (Fig. 4B, Clusters 4, 9, 11, and 12). A number of other clusters show a pattern of successive increased and decreased expression, indicating fluctuations in general metabolic activity (Fig. 4C, Clusters 5, 8, and 14). A final group of clusters display a pattern of a steady increase for the first 2 weeks and then a sharp decline by the 3rd week (Fig.  4D, Clusters 2, 3, 6, 7, and 15), indicating the end of various physiological processes (Fig. 1).
Clusters 2 and 3 contained the greatest number of genes that included the majority of the well known matrix genes (e.g. collagen types I, II, III, IV, V, VI, XI, and XII and bone sialoprotein, osteonectin, osteopontin, fibronectin, laminin, lumican, versican, tenascin, decorin, biglycan, and glypican), growth factors (IGF-I, TGF-␤, FGF-7), growth factor receptors (PDGF and NGF), transcription factors (hypoxia inducible factor 1␣, c-fos, and Sox9), as well as a very large number of genes representing other gene families (Supplemental Table 3). In  2. cDNA microarray analysis. Images of the same area from two microarray filters following hybridization with radioactively labeled RNA obtained from intact bone and PF day 10. Arrows indicate the varying signal intensity between the same spots on each membrane, indicating differential gene expression. addition, these two clusters also included over 600 functionally unknown genes (ϳ57%) represented either as novel or EST sequences or known genes with no assigned function.
Identification of the Wnt Signaling Pathway-Since there were many signaling molecules present in our subtracted cDNA library (Supplemental Table 1), we decided to examine the possibility of identifying active pathways that participate in bone regeneration. One of the more complete signaling pathways identified (based on the number of involved molecules) is the Wnt signaling pathway (31) (Fig. 5). More specifically, we identified Wnt-5A, Frizzled, casein kinase II, ␤-catenin, and phosphatase 2A, all of which were robustly up-regulated during the repair process (Table I). In addition, we also show the identity and expression levels of a number of genes that represent known transcriptional targets of the Wnt pathway (i.e. c-myc, fibronectin, retinoic acid receptor gamma, connexin 43, and OSF-2) (Table I). DISCUSSION The ability of bone to heal without leaving a scar makes it unique in the body. This is achieved through the intricate, interdependent stages of the healing process, which mirrors the tightly regulated development of the skeleton. Previously, it was shown that inflammation begins immediately following a fracture with the formation of a hematoma, which functions as a source of signaling molecules (e.g. growth factors, cyto-kines, etc.) that induce cascades of cellular events crucial to the repair process (32). Shortly thereafter, intramembranous ossification begins as a result of cellular changes within the periosteum resulting in the formation of new bone (hard callus) within a few millimeters of the fracture site (19). Simultaneously, the formation of cartilage begins at the fracture site (soft callus) and is followed by endochondral ossification that begins first with calcification of the cartilage and continues with its replacement by woven bone.
Despite this level of understanding of skeletal repair, there is surprisingly little known of the scale or complexity of the process at the molecular level. This study was designed to address this issue, by characterizing the global increases in gene expression that occur during the critical phases of bone fracture repair. As further evidence of how little is known of the healing process, this study demonstrated that a large number (1,243 or ϳ34%) of cDNAs identified through SSH represented functionally unknown genes. In addition, there were many known genes whose expression has never been described in bone. Inevitably, some of these genes will prove to hold central roles in other models of wound repair in the musculoskeletal system and possibly other tissues.
Clustering analysis revealed that the majority of these functionally unknown genes (ϳ663) were grouped within Clusters 2 and 3 that also contained ϳ389 known genes (Supplemental Table 3). The transcriptional profile pattern of these two clusters is marked by a sharp increase at PF day 3, which remains relatively high throughout the first 2 weeks, peaking at PF day 14, and then decreasing by week three. This high level of initial activity correlates well with the temporal sequence of biological events occurring during the early stages of fracture repair. As shown in Fig. 1, the highest level of metabolic activity predominantly occurs during the first 2 weeks following a fracture and is clearly associated with inflammation, chondrogenesis, and ossification, processes characterized by dramatic physical changes in cell type, number, size, shape, and connectivity, as well as changes in the ECM produced by the resident fibroblasts of the callus, chondroblasts/cytes, and osteoblasts/cytes. Furthermore, these data indicate that the genes grouped in Clusters 2 and 3 play a role in the growth and differentiation of the fracture callus, after which their role diminishes in the remodeling phase (PF day 21). Given the large number of known genes present in these two clusters, as well as the diverse functional families they represent, it would be premature and inaccurate to "assign" a possible function to the unknown genes. Clearly, further examination is required to help define the nature (i.e. structure and function) of some these novel genes.
In contrast, the large number of known genes identified in this study represent a multitude of families with functional involvement in the cell cycle, cell adhesion and communication, cytoskeleton, ECM, growth factors/cytokines, immune/inflammation, general metabolism (enzymes), muscle, protein/processing and degradation, RNA processing, signaling, transcriptional activation, and transport, etc., all indicating the complex, interdependent nature of the healing process. In addition, there were a large number of known genes without a defined function, which were labeled as miscellaneous (Supplemental Table 1).
Potentiating the possibility of identifying some of the key signal transduction pathways involved in the repair process, a large number of signaling molecules were found to be upregulated. The identification of these signaling molecules will enable us to assemble putative signal transduction pathways that are critical to the repair process. One clear example is the diverse Wnt signaling pathway (31) that involves many molecules, some that were identified as activated in this study (Wnt-5A, frizzled, casein kinase II, ␤-catenin, protein phosphatase 2A), as well as some of the transcriptional targets such as c-myc and fibronectin (46), connexin 43 (47), retinoic acid receptor ␥ (48), and OSF-2/periostin (49). Since it is already known that Wnt signaling is intimately involved in early embryonic development, organogenesis, and cell differentiation (31), and since fracture repair is essentially a replay of embryonic development, it is consistent then that members of the Wnt signaling pathway are activated. Interestingly, that the expression pattern of Wnt-5A increases sharply at PF day 3, declines (PF days 5-7), rises again at PF days 10 and 14, and declines back down to intact bone levels by PF day 21 suggests  that this signaling pathway may be activated during inflammation and chondrogenesis, but not during osteogenesis or remodeling. Evidence supporting the Wnt-5A role in chondrogenesis stems form previous work in skeletal development that showed that Wnt-5A is expressed in the limb-forming region during the earliest stages of development and that its misexpression delayed the maturation of chondrocytes and the onset of bone collar formation (50), as well as the truncation of long bones due to retarded chondrogenic differentiation (51). Supporting evidence of the role of Wnt-5A in the early stages of fracture healing (i.e. inflammation), in particular as a regulator of the diversification of hematopoietic cells, arises from a study that showed that the formation of macrophages was suppressed when bone marrow cells were exposed to Wnt-5A, while monocytes and red blood cells became the prevalent cell type (52). This finding is consistent with the formation of new blood vessels that occurs following a bone fracture (53), resulting in increased blood flow and thus increased numbers of red blood cells.
Finally, Wnt-5A has recently been shown to be a target of sonic hedgehog, suggesting a critical interaction of the Wnt and sonic hedgehog signaling pathways (54). Since it was previously shown that sonic hedgehog and Indian hedgehog were expressed in both chondrocytes and osteoblasts in the developing skeletal system, and thus hold important roles as regulators of skeletogenesis (55)(56)(57)(58), it is likely that these pathways are also critical to the fracture healing process. The hedgehog pathway was also shown to be linked to the well established BMP signaling pathway during skeletal development and repair (41, 56, 59 -60). Taken together, these data support the interaction of multiple signaling pathways (e.g. Wnt, hedgehog, BMP) during skeletal repair and suggest that these proteins may play a pivotal role in the successful healing of fractures.
In summary, this work demonstrates for the first time that the expression of thousands of genes is up-regulated during the repair of skeletal fractures in mammals, indicating that the processes of inflammation, chondrogenesis, and osteogenesis involved are as transcriptionally complex as they are spatially sophisticated. Although the majority of these genes are known, there remains a large number that are functionally novel. Considering the temporal activity of both known and unknown genes, a critical objective will be to molecularly and functionally characterize these genes, as a means of determining whether they are central to the repair process. Clearly, there is a great deal of information still to be grasped before the central physiologic stages are understood. Nevertheless, given the strong parallels of mammalian bone regeneration, skeletal development, and wound repair, this information will surely provide great insight toward some of the critical molecular mechanisms involved in a wide array of skeletal processes and dysfunctions.