Genome-wide expression profiles of subchondral bone in osteoarthritis

Introduction The aim of this study was to evaluate, for the first time, the differences in gene expression profiles of normal and osteoarthritic (OA) subchondral bone in human subjects. Methods Following histological assessment of the integrity of overlying cartilage and the severity of bone abnormality by micro-computed tomography, we isolated total RNA from regions of interest from human OA (n = 20) and non-OA (n = 5) knee lateral tibial (LT) and medial tibial (MT) plateaus. A whole-genome profiling study was performed on an Agilent microarray platform and analyzed using Agilent GeneSpring GX11.5. Confirmatory quantitative reverse-transcription polymerase chain reaction (qRT-PCR) analysis was performed on samples from 9 OA individuals to confirm differential expression of 85 genes identified by microarray. Ingenuity Pathway Analysis (IPA) was used to investigate canonical pathways and immunohistochemical staining was performed to validate protein expression levels in samples. Results A total of 972 differentially expressed genes were identified (fold change ≥ ± 2, P ≤0.05) between LT (minimal degeneration) and MT (significant degeneration) regions from OA samples; these data implicated 279 canonical pathways in IPA. The qRT-PCR data strongly confirmed the accuracy of microarray results (R2 = 0.58, P <0.0001). Novel pathways were identified in this study including Periostin (POSTN) and Leptin (LEP), which are implicated in bone remodeling by osteoblasts. Conclusions To the best of our knowledge, this study represents the most comprehensive direct assessment to date of gene expression profiling in OA subchondral bone. This study provides insights that could contribute to the development of new biomarkers and therapeutic strategies for OA.


Introduction
Osteoarthritis (OA) is the most prevalent musculoskeletal disorder and the most common form of arthritis among older individuals in most countries worldwide; OA constitutes a large economic burden due to the associated costs of medical care and lost wages [1,2]. Although degeneration of cartilage is the major characteristic of OA, the disease also involves the entire joint organ, including structural modifications of underlying subchondral bone, pathological changes of the meniscus and synovitis [3,4].
Maintenance of normal joint structure and function relies on load adaptation of the cartilage and bone. Disruption of the physiological relationship between these tissues contributes to the development of OA pathology [5].
Increasing evidence indicates that the subchondral bone, consisting of the subchondral plate and subchondral spongiosa, plays a major role in the initiation and progression of OA [6]. Magnetic resonance imagingbased visualization of the whole knee structure demonstrates that increased tibial subchondral bone volume is associated with severity of knee OA [7]. Kraus and colleagues demonstrated that subchondral bone texture can be used as a biomarker to predict progression of knee OA [8]. A study comparing two guinea pig strains indicates that an increasing rate of knee subchondral bone remodeling is associated with the process of cartilage deterioration [9]. Moreover, human bone cell culture studies have shown that factors released from bone cells might affect cartilage metabolism [10,11]. These studies provide insights into a temporal relationship between subchondral structural changes and alterations in articular cartilage during the development of OA; they also underscore the importance of elucidating the molecular changes in human subchondral bone to improve our understanding of the pathogenesis of OA.
Whole-genome microarray is a common technology for studying the behavior of many genes simultaneously. All of the gene expression microarray profiling studies in OA so far have been performed on human articular cartilage, meniscus or synovium; however, none have been performed on human subchondral bone tissue directly. Only one study has reported the gene expression profiles of subchondral bone in an early OA experimental mouse model [12]. A few studies have evaluated the gene expression profiles of distal trabecular bone from human OA [13,14], but this site is remote from subchondral bone rather than a reflection of alterations locally in the OA joint. This paucity of subchondral bone microarray studies is most probably due to the difficulties associated with isolation of high-quality RNA from subchondral bone. As described in our previous study [15], we have developed a method of precisely harvesting specific regions of subchondral bone tissue and for subsequently isolating high-quality RNA from these specimens. Our methodology has made it possible to perform microarray analyses of human subchondral bone samples.
Our goal in this study was to evaluate the association of subchondral bone gene expression with bone histomorphometry at sites of intact articular cartilage and osteoarthritic lesioned cartilage. To our knowledge, this is the first study to successfully perform microarray analyses of human knee subchondral bone in OA and non-OA samples, thereby providing clues to the pathogenic mechanisms of OA that could inform development of new diagnostic markers and therapeutic targets.

Human knee joint tissues
Human osteoarthritic tibial plateaus with medial compartment dominant OA and macroscopically normal lateral compartments were obtained during total knee joint replacement surgery from knee OA patients. The diagnosis of OA was based on the criteria for knee OA of the American College of Rheumatology [16]. Normal tibial plateaus were obtained at autopsy or within 8 hours after amputation surgery from donors with nonarthritic knee joints.
To ensure consistency of sampling of prespecified regions of interest, the anatomic orientation was indicated on the freshly isolated specimens by marker pen and then all specimens were stored immediately in liquid nitrogen. Methods for precisely obtaining specimens of a particular type from specific regions of interest of the joint have been described previously [15]. The tissues were manipulated at all times under liquid nitrogen to prevent RNA degradation. Regions of interest were harvested for histological analysis of the articular cartilage, and subchondral bone structural parameters were determined by micro-computed tomography; regions immediately adjacent were harvested for RNA isolation and microarray analysis. For microarrays, RNA was isolated from subchondral bone from four regions: the outer and inner lateral tibial (LT) plateau and the inner and central medial tibial (MT) plateau (Additional file 1). The gene expression intensities changed little within the LT plateau (between inner and outer LT regions) and within the MT plateau (inner and central MT regions) (Additional file 2); we therefore based the analyses in the current study on a comparison of gene expression from subchondral bone of 20 outer LT plateaus and 20 MT plateaus from knee OA specimens. These 40 samples were derived from 30 individual knee OA specimens, yielding 10 paired LT and MT samples and 10 unpaired LT and MT samples. The 40 OA samples chosen for microarray analysis represented a subset from our previous study [15], and were chosen to have LT OARSI scores ≤8 and MT OARSI scores ≥18. In addition, five LT and five MT paired subchondral bone samples from non-OA donors were used as controls for the microarray analysis. The characteristics of the samples used for microarray analysis are summarized in Table 1 and Additional file 3.
The study was approved by the institutional review board of Academia Sinica, Tri-Service General Hospital and Taipei Medical University Hospital, and written informed consent was obtained from all of the participants.

Subchondral bone tissue harvest and RNA isolation
All tissues were sectioned to generate osteochondral slabs 2 to 2.5 mm thick from regions of interest for microcomputed tomography scanning and histological analysis ( Figure 1); regions immediately adjacent were isolated for subsequent grinding and RNA isolation using a custommade workstation as described previously [15]. Approximately 100 mg powdered subchondral bone from the tibial plateau was obtained for RNA extraction. Five milliliters of Trizol (Invitrogen, Carlsbad, CA, USA) were added to 100 mg (dry weight) ground tissue powder and mixed by vortexing until homogeneous. A total of 1 ml chloroform was added to the homogenized sample and centrifuged at 12,000 × g for 15 minutes at 4°C to achieve phase separation. The clear aqueous phase was transferred to a fresh tube for a second phase separation with phenol/ chloroform/isoamyl alcohol (25:24:1; Sigma-Aldrich, St. Louis, MO, USA), followed by RNA precipitation with pure isopropanol. After washing with 70% ethanol, the air-dried pellet was dissolved in RNase-free water and genomic DNA was removed by DNase digestion (Qiagen, Hilden, Germany) according to the manufacturer's instructions. To further improve RNA quality, one-half volume of 7.5 M LiCl (Invitrogen) was added and incubated at −20°C for at least 30 minutes followed by RNA precipitation by centrifugation at 12,000 × g for 15 minutes. After 70% ethanol washing, the RNA was resuspended in RNase-free water. The RNA concentration and purity were determined by Nano-Drop (NanoDrop Technologies, Wilmington, DE, USA). The RNA integrity number and the 28 s/18 s ratio were estimated using RNA 6000 Nano Assays on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer's instructions.

Microarray analysis
Four hundred nanograms of total RNA per sample were used for one round of cRNA synthesis and amplification. Cyanine 3-labeled cRNA was then purified and hybridized to Agilent whole human genome 44 k microarray chips (Agilent Technologies). All procedures were carried out according to the manufacturer's instructions (Agilent Technologies, Taipei, Taiwan). The array signal  intensities were further analyzed by the Agilent Gene-Spring GX software (version 11.5; Agilent Technologies).
Gene expression values of all datasets were normalized using quantile normalization; probes with low signal intensities were excluded by setting the filter above 32.
The normalized values were log transformed. Significant differentially expressed genes between lateral and medial samples in OA patients were identified using an adjusted Student's t test (P <0.05) corrected for multiple comparisons with the Benjamini-Hochberg test. The results were used to run principal component analysis to uncover the trends among OA and normal samples. Differentially expressed genes between medial and lateral regions among groups were identified by a threshold of ≥2 fold-change and P ≤0.05. The hierarchical clustering method, with Euclidean distance and centroid linkage, was further applied to cluster the differences in gene expression levels between samples. Lists of differentially expressed genes were imported into ingenuity pathway analysis (IPA; Ingenuity Systems, Redwood City, CA, USA) to identify functional annotations and biological interactions from the many known relationships between proteins, genes and diseases. The biological interaction scores were defined by the IPA statistical algorithm based on its knowledge base, and the adjusted P value was calculated by Fisher's exact test and corrected for multiple comparisons with the Benjamini-Hochberg test. All microarray raw data are available through the GEO database [GEO:GSE51588].

Quantitative reverse-transcription polymerase chain reaction validation
To validate results from microarray analysis, we performed reverse-transcription polymerase chain reaction (qRT-PCR) for 85 genes (plus GAPDH) on nine additional OA paired LT and MT subchondral bone specimens; eight of these additional OA knee joint specimens have been described previously [15], as well as the qRT-PCR results for some of the genes combined with data yielded from one additional knee OA specimen. Forty-three of these genes were identified for verification on the basis of significant fold-change (≥2) comparing LT and MT microarray results. The remaining 42 genes were identified on the basis of their potential involvement in OA as described previously [15]. The expression of these 42 genes was previously analyzed for seven regions of interest from eight OA cartilage and subchondral bone specimens.
In the current study, we analyzed two regions of interest (LT and MT subchondral bone corresponding to outer LT and center MT of the prior work) in nine specimens (the eight former specimens and an additional OA specimen with available LT and MT regions). The qRT-PCR was performed using the Taqman high-density microfluidic cards (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. RNA from each region was converted into cDNA using Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA). The qRT-PCR reactions were performed on the ABI Prism 7900HT Sequence Detection system and the fluorescent signal intensity was analyzed by Sequence Detector v2.3 software (Applied Biosystems, Foster City, California). The levels of cDNA among samples were normalized to the expression of GAPDH and analyzed with the ΔΔCt relative quantification method to identify significant variation in gene expression comparing LT and MT regions in OA samples. Two-tailed unpaired Student t tests were performed to evaluate for statistically significant differences in gene expression levels between regions.

Immunohistochemistry staining
Approximately 2 mm diameter sections from the tibial plateau were fixed in 4% paraformaldehyde (Sigma-Aldrich) overnight and decalcified in 10% ethylenediamine-N, N, N′, N′-tetraacetic acid, disodium salt, dehydrate (Sigma-Aldrich) for 2 weeks. After decalcification, the sections were embedded in paraffin and 10 μm sections were prepared. The embedded tissue sections were then deparaffinized, hydrated, and the endogenous peroxidase activity was quenched in 0.3% H 2 O 2 -methanol for 30 minutes. After incubating in 0.5% trypsin/1% CaCl 2 at 37°C for 30 minutes and blocking with diluted goat serum at 25°C for 1 hour, the sections were incubated with anti-human PSOTN (dilution 1:300) antibodies (ab14041; Abcam, Cambridge, MA, USA), and LEP (dilution 1:1,000) antibodies (ab16627; Abcam) or nonimmune goat serum at 4°C overnight and reacted with biotinylated secondary antibody for 30 minutes. Signal amplification and staining were performed using Vectastain ABC kits (Vector Laboratory, Burlingame, CA, USA) according to the manufacturer's protocol and counterstained with hematoxylin solution Gill no. 3 (Sigma-Aldrich).

Analysis of differential gene expression
For the microarray analysis, 972 differentially expressed genes with ≥2 fold-changes between OA MT and OA LT regions were identified; 420 of these genes were upregulated and 552 were downregulated in OA MT compared with OA LT (Additional file 4). In total, 640 genes were differentially expressed with greater than twofold expression level changes (308 upregulated and 332 downregulated), 163 genes with changes greater than threefold (65 upregulated and 98 downregulated), 111 genes with changes greater than fourfold (28 upregulated and 83 downregulated), 35 genes with changes greater than sixfold (11 upregulated and 24 downregulated), and 23 genes with changes greater than eightfold (eight upregulated and 15 downregulated) (Figure 2A). A graphical depiction of microarray gene expression results across the entire extent of the tibial subchondral bone can be seen in Additional file 4.
Unsupervised hierarchical clustering was performed on the differentially expressed genes to identify the groups or clusters of genes based on similarities between the expression profiles of the samples. On the basis of gene expression profiles, OA LT and OA MT samples were grouped into two distinct clusters. However, the pattern for LT and MT samples from non-OA donors was mixed and did not cluster by sample site; the clustering of all non-OA samples was more akin to the OA LT samples than to the OA MT samples ( Figure 2B).

Quantitative reverse-transcription polymerase chain reaction validation
Nine joint tissue samples were available for validation of microarray results; qRT-PCR analysis of LT and MT subchondral bone regions of these samples was performed on 86 selected genes with GAPDH as endogenous control. Based on the results of histological evaluation and microcomputed tomography scanning, the subchondral bone samples selected from the OA MT regions displayed characteristics of severe stages of OA; whereas the OA LT regions displayed characteristics of relatively normal joints. The mean expression pattern of the LT region was compared with the mean expression pattern of MT regions. By qRT-PCR, 77 of 85 genes showed a statistically significant difference between regions (Additional file 5). Expression by qRT-PCR matched the pattern (upregulation or downregulation) of expression observed by microarray, and the two datasets were strongly correlated (R 2 = 0.5764, P <0001; Figure 2C). The most differentially upregulated genes in OA MT regions were IL11, POSTN, ASPN, COL6A3, TUBB3, COL3A1, OGN, DIO2, PPEF1, and TNFSF11. The most differentially downregulated genes in OA MT regions were LEP, APOB, ADH1B, CCL8, TUSC5, KLB, NPR3, PLIN1, PCK1, and MYO16. All of these genes were validated by qRT-PCR except for ADH1B and MYO16, because these two genes were not available on the microfluidic cards used for qRT-PCR (Table 2).

Functional and pathway analysis
IPA of gene expression profiles revealed 279 functional networks (canonical pathways) (Additional file 6) that were significantly associated with the differentially expressed genes. The predominant associated networks and biological functions were related to lipid metabolism and mineral metabolism (Additional files 7 and 8), connective tissue disorders, cellular growth and proliferation, connective tissue development and function (Additional file 9). Most of these pathways were associated with OA pathogenesis such as bone remodeling, adiposity, connective tissue, and cell activity and death (Table 3). Interestingly, on performing pathway analysis with the most upregulated and downregulated genes, a novel signaling network in OA was identified; this network was generally characterized by upregulation of bone mineralization or collagen-associated genes, such as COL3A1 (upregulated 7.33-fold, ranking sixth in the top upregulated gene list), COL6A3 (upregulated 8.27-fold, ranking fourth in the top upregulated gene list), BMP1, BMP7, POSTN (upregulated 15.65-fold, ranking second in the top upregulated gene list), WISP1, HTRA1, SOST, and ITGA11, and the downregulation of genes associated with cellular metabolism, proliferation or differentiation (NMB, LEP (downregulated 23.27-fold, ranking first in the top downregulated gene list), CHRDL2, GRB14, CIDEC, CIDEA, BTG2, PLAC8, and PRDM16). Additional differentially expressed genes in this network included ANGPTL1 (upregulated 2.16-fold), which is involved in angiogenesis. Both POSTN and LEP may affect the phosphorylation of AKT signaling ( Figure 3A).

Validation of protein expression in subchondral bone specimens
The regions of interest of subchondral bone samples ground for the gene expression study were about 5 mm deep under the calcified cartilage, which included the subchondral bone plate and subarticular spongiosa. These regions have a very heterogeneous complement of cells, including osteoblasts, osteoclasts, osteocytes and bone marrow cells. Although the specific cell types contributing to the changes in gene expression cannot easily be confirmed, all of the cell types in the subchondral bone would be expected to contribute to the subchondral bone gene expression profile. To obtain further independent validation of the microarray results, immunohistochemical staining was performed to investigate the protein expression of proteins Periostin and Leptin corresponding to genes POSTN and LEP in human OA and non-OA tissue sections from LT and MT regions. The highly expressed Periostin could be found in osteoblasts beneath the tidemark, in osteocytes, and in lining cells of the OA MT regions, but not OA LT regions ( Figure 3B,C). Protein expression of Leptin was low in OA MT regions, but strong in OA LT regions ( Figure 3D,E). By immunohistochemical analysis, the proteins Periostin and Leptin were not differentially expressed in the MT and LT regions of non-OA samples (Additional file 10).

Discussion
To the best of our knowledge, this represents the first comprehensive whole genome expression profiling of non-OA and OA subchondral bone. The validity of the microarray results was confirmed by qRT-PCR for a selected subset of genes. Based on analysis by IPA, we identified biological functions implicated in the pathogenesis of OA in bone. Many of the functional networks identified in this study were also associated with processes in connective tissue development and function.
In our previous study [15], we evaluated the expression of 61 OA genes and correlated the expression patterns with the bone morphometric measurement (bone volume). The expression of a total of 45 of 61 genes in the subchondral bone was significantly correlated with the alteration of bone structure including ADAMTS1, ASPN, BMP6, BMPER, CCL2, CCL8, COL5A1, COL6A3, COL7A1, COL16A1, FRZB, GDF10, MMP3, OGN, OMD, POSTN, PTGES, TNFSF11, WNT1 and others, including the ratio of OPG (also known as TNFRSF11B) to RANKL (also known as TNFSF11). RANKL is primarily produced from osteocytes, osteoblasts and/or marrow stromal cells and is a key factor for osteoclast differentiation and activation. RANKL binds to OPG, also known as osteoclastogenesis inhibitory factor, and interacts with the receptor (RANK) for RANKL, to modulate the level of osteoclast activity and regulate bone homeostasis in response to various endocrine, paracrine, cytokine and mechanical stresses [17,18]. As shown in Additional file 11, significant correlations were observed for OPG/RANKL and the bone parameters including percent bone volume, bone structure (structure model index), trabecular thickness and trabecular space. These data confirmed our hypothesis that investigations of subchondral bone gene expression changes could provide clues to the pathogenic mechanisms of OA and inform development of new diagnostic markers and therapeutic targets. We focused on the outer LT and central MT regions for the purposes of comprehensive microarray analyses in the current study. The characteristics of the central MT region generally included cartilage matrix loss, cyst formation within cartilage matrix, denuded sclerotic bone, and thickened subchondral bone plate, consistent with a bone sclerosis phenotype [19], whereas the outer LT region was relatively normal. Many of the most significantly regulated genes identified by this microarray study have documented roles in the pathogenesis of OA, arthritis or bone formation including seven of the top 10 upregulated genes (POSTN, ASPN, COL6A3, COL3A1, OGN, DIO2, TNFSF11) and two of the top 10 downregulated genes (LEP, APOB) [20][21][22][23][24][25][26][27]. Through IPA we could also identify the involvement of a lipid metabolism protein network and a mineral metabolism protein network that involved 29 and 22 differentially expressed genes respectively; these networks included PLINE, LIPE, DGAT2, ADRB1, NPY1R, HCAR3 and P2RRY14 that changed −9.6-fold, -8.7-fold, -7.8-fold, -6.2-fold, -7.4-fold, -5.1-fold, and −5.1-fold respectively comparing OA MT with OA LT (Additional files 7 and 8).
These functional networks indicated that the bone cells in the OA MT regions were characterized by a condition of low energy production that included a low rate of mineral metabolism. These results are compatible with a low rate of bone remodeling in MT regions. Increased bone remodeling is an important predictor of OA progression [8] and is characteristic of early stages of OA including resorption of bone and increased porosity in the subchondral bone region [28]. Given that OA MT regions of interest could be defined as representing late stages of OA, we were interested to determine whether OA LT regions of interest could be described as representing early stages of OA. The OARSI histopathological scores of all OA LT samples were <6, representing a relatively intact overlying cartilage matrix with minimal superficial fibrillation not too dissimilar from normal. However, in contrast to normal samples, the bone structural parameters of the OA LT regions had a lower bone volume, thinning of the subchondral bone plate, increased porosity in the subchondral region and reduced bone density (bone density data not shown) suggestive of early OA. Another indicator of early OA, as demonstrated by animal models of OA and in OA patients, is increased bone remodeling characterized by a decreased expression ratio of OPG to RANKL [29,30]. In contrast to non-OA tissues, a reduced OPG to RANKL ratio was identified in OA LT regions compared with non-OA LT regions (OPG lacked significant change, but RANKL was upregulated 1.63-fold, one-tailed unpaired Student t test P = 0.048), consistent with early OA. Thus, based on the bone parameters and OPG to RANKL ratios, comparing LT with non-OA we concluded that the bone gene expression patterns of OA LT regions were consistent with an early-stage OA phenotype.
According to a hypothetical model of OA pathogenesis proposed by Burr and Gallant [6], joint loading increases bone remodeling in the early stage of OA, and a bonecartilage crosstalk may occur via diffusion of small molecules or increased vascularization at the deep layers of cartilage that could interfere with the normal collagen network. These changes would lead to a loss of aggrecan, increasing inflammation, and result in bone sclerosis at late stages of OA. Several of the genes identified in this study fit nicely into this hypothetical model. For instance, upregulation of SOST was identified (2.51-fold); SOST is an inhibitor of Wnt signaling that responds to mechanical loading in the proximal rat tibia and is associated with bone mass changes [31]. In addition, several genes relevant to cell differentiation and activity of osteoblasts, osteoclasts and osteocytes were identified in both a mouse load-regulated gene expression model [32] and in this human study; examples included ASPN (regulator of osteoblast collagen mineralization), WNT5A (agonist of WNT signaling pathway in osteoblasts and osteocytes) and VCAN (regulating transforming growth factor beta expression in osteoclasts). These findings suggest that mechanical loading plays an important role in structural changes of human subchondral bone, initiates bone remodeling, and contributes to the pathogenesis of OA. Moreover, our results could provide a possible novel molecular mechanism to explain changes in bone remodeling during OA development based on an AKT-regulated network. AKT is a serine/threonine-specific protein kinase that plays a pivotal role in many cellular processes such as metabolism, apoptosis, cell proliferation, transcription, cell migration and extracellular matrix alternation [33]. In an early-stage of OA, overexpressed LEP may increase phosphorylation of AKT [34], which will trigger the downstream signaling pathway to increase bone remodeling. In a late stage of OA, upregulated POSTN may inhibit phosphorylation of AKT [35] and will decrease the cellular metabolism, cell proliferation, and differentiation; the rate of bone remodeling would thereby decline.
One of the limitations of this study was that the subchondral bone samples were obtained from OA joints at end-stage disease due to the difficulty of obtaining joints with early stages of OA from humans. However, comparing our results with the gene expression profiles of subchondral bone in an early OA experimental mouse model [12], their reported fold-changes in SB gene expression of ASPN, CCL2, COL3A1, COL5A1, POSTN and TNFSF11 at the initial stages of OA were similar to those of our study. This result strongly supports the credibility of this study. Another limitation was the difficulty distinguishing whether differentially expressed genes reflected OA progression in cartilage, differences in innate bone structure, or were driven by changes in mechanical loading. In our currently study, significant correlation between overlying cartilage integrity and underlying subchondral bone structural parameters could be identified (Additional file 12), suggesting that differential gene expression across the subchondral bone is due, at least in part, to differences in cartilage integrity of the LT and MT regions in OA. Microarray analysis of additional non-OA control samples could help to exclude differentially expressed genes driven by innate bone structure. However, to exclude genes driven by mechanical loading, microarray analysis of lateral compartment dominant OA is necessary.
The site of initiation of OA, at bone or cartilage, has been disputed for decades. Numerous clinical and experimental studies have confirmed that increased bone volume and changes in bone quality of the tibial subchondral bone of the knee are related to loss of cartilage integrity [7,8,[36][37][38]. Goldring and Goldring have pointed out that subchondral bone responds more rapidly to adverse loading and events than cartilage [5]. Moreover, the integrity of subchondral bone predicts later worsening of radiographic osteoarthritis progression [8]. As summarized by Martel-Pelletier and Pelletier, an abnormal rate of bone remodeling, incomplete bone mineralization and increased osteoid collagen matrix will result in hypomineralization of subchondral bone and deterioration of cartilage as manifested by knee OA progression [39]. These data would suggest that drugs designed to target both the bone and cartilage compartments to optimize bone remodeling rates and inhibit cartilage loss in the early stages of OA could inhibit progression of disease. The recent success of strontium ranelate as a disease-modifying agent for OA provides the first clear evidence that this hypothesis has been confirmed [40]. Some other bone-targeting agents have been assessed as OA therapeutics, including bisphosphonates [41,42], calcitonin [43] and vitamin D [44]. These bone-acting agents are unlikely to be effective in late stages of OA, when remodeling is already suppressed. The development of diagnostic biomarkers to detect progressive changes in early stages of OA is therefore important. The genes highly regulated in early stages of OA in our microarray data could provide potential diagnostic biomarkers or therapeutic targets for OA.

Conclusion
Our microarray results were obtained by examining early and late stages of OA and non-OA samples. Most of the differentially expressed genes in bone are involved in cartilage and bone development, OA pathogenesis and bone remodeling in the early and late stages of OA. This direct assessment of gene expression profiling in OA subchondral bone provides a comprehensive understanding of disease pathology that could further the development of new OA biomarkers and therapeutic strategies.